Revisiting the bi-dimensional elastic three-layer beam analytical solution: a novel solution for various loads followed by an XFEM numerical verification

Abstract

An analytical solution based on Airy’s stress function is presented for a three-layer sandwich beam in a cantilever configuration, subjected to transverse force, bending moment, and linearly distributed load. Each layer is assumed to be isotropic and homogeneous. Closed-form expressions for displacements and stresses are derived, with the displacement field exhibiting a fifth-degree polynomial dependence. Classical solutions for two-layer and homogeneous beams are recovered as particular cases. The proposed solution is implemented and made publicly available as an open-source library on GitHub (https://github.com/rodrgz/2D3LayerCantileverBeam), released under the MPL license. To assess its accuracy, the extended finite element method (XFEM)-well-suited for problems with material interfaces-is employed. Numerical results show excellent agreement with the analytical predictions, demonstrating the proposed solution’s potential as a reliable benchmark for validating numerical methods in elasticity problems involving layered media.

Introduction

Sandwich beams are essential elements of various engineering structures in civil, mechanical, and aerospace engineering, owing to their high specific strength, damping capacity, and low weight [25]. Failures in these beams can be attributed to interfacial stress concentrations and singularities arising from the disparity in stiffness properties among their constituent layers, particularly in critical operating conditions [1].

The historical motivation behind examining such structures, as elucidated by [12], was derived from the observation that thin-walled aluminum alloy structures, widely utilized in aircraft manufacturing during World War II, manifested surface irregularities under considerable loads. To address the buckling, engineers introduced a sandwich configuration comprising two thin layers of high-stiffness facing material with a thick layer of ultra-lightweight core in between. Nevertheless, conventional theories such as Euler–Bernoulli and Timoshenko beam theory [33] were unable to accurately calculate the instabilities of sandwich beams due to the pronounced shearing deformations enabled by the lightweight core.

Before continuing, it is important to recognize the extensive body of literature, as emphasized by [15] in the early 1970s. Consequently, in this review, the main emphasis will be placed on studies that attempt to derive a two-dimensional analytical solution for sandwich beams in accordance with the theory of elasticity, in which each layer consists of an isotropic material.

Moving forward, the analysis of plane problems of beams is particularly crucial when considering thick beams or assessing the accuracy of refined shear deformation theories or numerical solution methods [36]. The Airy’s stress function method is typically employed to determine beam stresses and displacements in this context. Early work in this area was conducted by [33], who examined a range of scenarios involving isotropic beams. The study was later extended by [18] to encompass anisotropic beams. Nevertheless, the availability of exact elasticity solutions for the bending of sandwich beams is significantly scarce in the literature [11].

In the investigation of bending phenomena in sandwich beams, [15] highlighted the noteworthy contributions of [12] and [35]. Despite the former being widely considered pioneering, [5] drew attention to the overlooked endeavors of [17], likely due to the impact of World War II and the subsequent Cold War.

The focus of [12] was primarily on computing deflections and buckling loads, assuming the Euler–Bernoulli beam hypothesis and supposing that the bending moment was solely resisted by the faces. In contrast, the work by [35] accounted for the effects of transverse-shear deformation and rotatory inertia in both the core and faces of the sandwich structure. The methodology employed was an extension of the one proposed by [19], and the obtained results converged to Mindlin’s equations for homogeneous plates under certain limiting conditions. Furthermore, the solution could be simplified, yielding results that coincided with those obtained by [12]. A shared feature among these studies was that the total deformation comprised a component caused by bending and an extra contribution arising from the shear of the core.

The work of [28] proposed an elasticity solution wherein the elastic constants are functions of the transverse direction. This solution can be utilized to determine the stress distribution in a beam subjected to an end load or distributed loading, though the article only presents the solution for the former case. It is worth noting that in this solution, no interfacial conditions between the layers are assumed.

Expanding on the findings of Schile, [30] leveraged those results to derive a solution for a cantilever sandwich beam with transverse properties that vary in a piece-wise constant fashion. Inversely, [10] extended the concept of Airy’s stress polynomial functions in homogeneous beams to composites, taking into consideration interfacial conditions. They also gained a closed-form elasticity solution for a sandwich cantilever beam composed of three isotropic layers and subjected to an end load.

Moreover, [18] presented an exact solution for a cantilever sandwich beam that was subjected to a transverse force and a moment using Airy’s stress polynomial functions. The analysis assumed each layer to be orthotropic in the plane of bending. Lekhnitskii further specialized the solution for a two-layer unsymmetrical laminated cantilever beam. Moreover, [23, 24] formulated exact solutions for the elasticity of cylindrical bending of laminated composite plates, considering unidirectional laminates as well as two- and three-layer cross-ply laminates subjected to sinusoidal loads.

Additionally, [16] proposed a theory derived from the variation approach that is suitable for analyzing sandwich beams with arbitrary boundary conditions and possessing two outer layers of equal thickness. The established theory presumes that Bernoulli’s hypothesis is applicable to each layer individually but not to the entirety of the cross-section. It is shown that the generalized displacement can be chosen in a manner that renders the set of equations governing motions in which the beam remains straight independent of the set of equations describing bending and shear motions.

One noteworthy three-dimensional elasticity investigation was conducted by [20] for heterogeneous bars in flexion. In this analysis, it was postulated that all layers had the same Poisson’s ratio, thus simplifying the calculations.

In their work, [25] presented a plane elasticity framework to develop solutions for unsymmetrical sandwich beams, also employing Airy’s stress function. Their methodology was formulated to accommodate n layers with different dimensions and material properties, thus providing a general framework, though only a closed-form solution for a three-layer beam subjected to a uniformly distributed load was presented. Vertical forces were applied at the beam ends to enforce equilibrium. Additionally, they proposed a straightforward procedure for selecting an appropriate Airy’s stress function to address sandwich beams subjected to polynomially distributed loads.

The study conducted by [6] presents a three-dimensional elasticity solution for symmetric cantilever sandwich beams, where each layer of the sandwich beam is assumed to possess homogeneity, transverse isotropy, and unidirectional reinforcement, with the end section subject to a transverse force and a torque.

It is worth noting that in recent years there has been a surge in interest in functionally graded materials (FGM). A great number of articles have been published, primarily focusing on anisotropic beams whose material properties fluctuate arbitrarily in the direction of the thickness [13], as well as functionally graded cantilever beams [37]. A noteworthy contribution is the work of [34], in which they provided an analytical solution for a sandwich beam with a graded intermediate layer subjected to a uniform load on the upper surface. The authors also accounted for variable mechanical properties, with the Young’s modulus of the graded intermediate layer assumed to be an arbitrary function of the thickness coordinate while maintaining a constant Poisson’s ratio. See also the works of [2], where the variational asymptotic method was used to analyze the deformation behavior of bidirectional FGM beams under transverse loading and of [32], where the effect of axial force on the free vibration and buckling of FGM sandwich beams is studied.

Recent developments include the application of a fifth-order shear and normal deformation theory by [21] for laminated and sandwich plates under transverse loads using Navier’s method. Pagani et al. [22] proposed a solution for simply supported cross-ply laminated and sandwich beams using a layer-wise framework based on Lagrange polynomials, where the three-dimensional displacement field is represented by arbitrary-order approximations of the displacement variables within each layer. Sidhardh and Ray [29] presented a closed-form exact solution for the size-dependent elastic response of simply supported laminated beams, based on strain gradient elasticity and incorporating both isotropic and orthotropic laminae. In [3], explicit analytical estimates are provided for the through-the-thickness normal stress component (${\sigma }_{22}$) of sandwich beams that have identical skins and are governed by the Krajcinovic theory, which relies on zigzag warping under flexure.

For further insights into the bending, buckling, and free vibration solutions in laminated composites and sandwich structures, [11] and [27] extensively review various theories, including single-layer, layer-wise, zigzag, and exact elasticity solutions, as well as the utilization of finite element modeling.

The proposed solution in this work can be seen as an extension of [25], considering a cantilever configuration, a distributed linear load, and incorporating other loadings such as a bending moment and a transverse load at the end section. With the advent of computing, which has dramatically improved approximation methods and the systems of equations they use, more complex analytical solutions can also be obtained by means of symbolic mathematics and provided as programming code.

Therefore, the implementation of the analytical solution derived here is made available as a programming library via our GitHub repository at https://github.com/rodrgz/2D3LayerCantileverBeam. This repository hosts implementations in both Julia and FORTRAN, ensuring accessibility across different knowledge bases in programming languages. More information is in Sect. 5.

Although several analytical solutions for layered and sandwich structures have been proposed in the literature, the set of two-dimensional problems commonly employed to assess the efficacy of the extended finite element method (XFEM)—as well as other methods for modeling material interfaces—remains rather limited. Typical examples include the bi-material plane bar problem with a vanishing Poisson’s ratio, exhibiting linear displacements, and the bi-material plate with a circular inclusion, whose solution in polar coordinates involves rational terms. Nevertheless, a problem with cubic polynomials in the primal variable recently addressed in the literature, specifically discussed in [8] is the sandwich beam. Accordingly, this work also aims to extend the solution presented in the aforementioned paper and develop an analytical solution for a two-dimensional sandwich beam subjected to more intricate loadings, including a distributed linear load, bending moment, and transverse load in the end section, resulting in quintic polynomials in the primal variable. This solution will be represented by Airy’s stress functions with a quadratic polynomial degree for each layer.

Analytic model

This section presents the derivation of a two-dimensional solution for a three-layer sandwich beam with a combination of transverse force, bending moment at the free end, and a linearly distributed load applied to the top, as illustrated in Fig. 1. This type of structure is commonly employed in engineering applications, such as laminated or coated composites. The solution is determined based on the stress formulation for elasticity theory, utilizing the stress function concept. For this two-dimensional problem, the in-plane stresses are determined using a single function—the Airy’s stress function—which satisfies the equilibrium and compatibility relations through a single partial differential (bi-harmonic) equation. Moreover, the layers are completely bound to each other.

In the context of Airy’s stress solutions for a sandwich beam with three layers, closed solutions for this case are provided by two works: [10] for a cantilever beam subjected to an end load and [25] for a beam with constant distributed loading and, at the end, the resulting transversal forces to auto-equilibrate the beam. In relation to the latter, it would still be necessary to make considerations using Saint-Venant’s principle to address the cantilever configuration. Regarding the loads, the proposed solution presents, in addition to the transverse load, a bending moment at the end section, and a distributed linear load on the upper part of the beam. Thus, the proposed solution can be seen as an extension of these works, for which the principle of superposition is valid.

[See PDF for image]

Fig. 1

Geometrical model and boundary conditions of a three-material, two interfaces, cantilever beam. Consider a narrow rectangular cross-section with a width of one unit

Fundamental linear elasticity equations

Consider a linear elastic body under small deformations with no body forces. The strong form of the boundary value problem (BVP) is as follows:

Problem 1. Given a body force $\overline{\mathit{b}}:\mathrm{\Omega }\to {\mathbb{R}}^2$, Dirichlet condition $\overline{\mathit{u}}:{\Gamma }_u\to {\mathbb{R}}^2$ and Neumann condition $\overline{\mathit{t}}:{\Gamma }_t\to {\mathbb{R}}^2$, find $\mathit{u}:\overline{\mathrm{\Omega }}\to {\mathbb{R}}^2$, such that

1a

$$\begin{array}{cc}\hfill \mathbf{\nabla }·\mathit{\sigma }+\overline{\mathit{b}}=\mathbf{0},& \forall \mathit{x}\in \mathrm{\Omega },\hfill \end{array}$$

1b

$$\begin{array}{cc}\hfill \mathit{u}=\overline{\mathit{u}},& \forall \mathit{x}\in {\mathrm{\Gamma }}_u,\hfill \end{array}$$

1c

$$\begin{array}{cc}\hfill \mathit{\sigma }·\mathbf{n}=\overline{\mathit{t}},& \forall \mathit{x}\in {\mathrm{\Gamma }}_t,\hfill \end{array}$$

1d

$$\begin{array}{cc}\hfill 〚\mathit{\sigma }〛·{\mathbf{n}}_{\mathrm{\Gamma }\mathit{ij}}=\mathbf{0},& \forall \mathit{x}\in {\mathrm{\Gamma }}_{\mathit{ij}},\hfill \end{array}$$ where $〚\bullet 〛={\bullet }^{+}-{\bullet }^{-}$ is the jump operator describing the interface, and $\mathit{\sigma }$ is the Cauchy stress tensor defined by $\mathit{\sigma }=\mathbb{E}\mathit{\epsilon }$, with $\mathbb{E}$ being the fourth-order constitutive elastic tensor and $\mathit{\epsilon }$ the infinitesimal elastic strain tensor defined by

2

$$\begin{array}{c}\hfill \mathit{\epsilon }=\frac{1}{2}\left[\mathrm{\nabla },\mathit{u},+,{\left(\mathrm{\nabla },\mathit{u}\right)}^{\text{T}}\right].\end{array}$$For isotropic materials, Hooke’s law gives the stress–strain relations

3

$$\begin{array}{c}\hfill \mathit{\epsilon }=\frac{1+\nu }{E}\mathit{\sigma }-\frac{\nu }{E}\text{tr}\left(\mathit{\sigma }\right)\mathbf{I},\end{array}$$where: $\mathit{\epsilon }$ is the infinitesimal strain tensor; E the Young’s modulus; $\nu$ the Poisson’s ratio; $\text{tr}\left(\mathbf{·}\right)$ the trace of the tensor; and $\mathbf{I}$ the second-rank identity tensor.

The strain tensor is related to the displacements $\mathit{u}$ through the symmetric part of its gradient, as expressed in Eq. 2. Consequently, for a strain tensor to be integrable and an associated displacement field to exist, strain compatibility equations must be hold for simply connected regions, which are given by

4

$$\begin{array}{c}\hfill \mathbf{\nabla }\times {\left(\mathbf{\nabla },\times ,\mathit{\epsilon }\right)}^T=\mathbf{0}.\end{array}$$These are known as Saint-Venant’s compatibility equations. In three dimensions, they form a system of 81 relations, of which only six are independent. In two dimensions, only one equation is nontrivial, corresponding to $i=j=1$ and $k=l=2$, which is

5

$$\begin{array}{c}\hfill \begin{array}{c}\frac{{\partial }^2{\epsilon }_{11}}{\partial y^2}+\frac{{\partial }^2{\epsilon }_{22}}{\partial x^2}-2\frac{{\partial }^2{\epsilon }_{12}}{\partial x\partial y}=0\end{array}.\end{array}$$The strains can be eliminated from Eq. 4 using Eq. 3 and by manipulating Eq. 1a, resulting in the compatibility equations being expressed only in terms of stress as

6

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2\mathit{\sigma }+\frac{1}{1+\nu }\mathrm{\nabla }\mathrm{\nabla }\text{tr}\left(\mathit{\sigma }\right)=\mathbf{0},\end{array}$$also known as Beltrami–Michell compatibility equations.

Airy’s stress function

By neglecting all components of Maxwell’s stress function tensor except for the scalar component ${\mathrm{\Phi }}_{33}$, see Appendix 5, the system of equations associated with the compatibility statement, Eq. 46, reduces to a single partial differential equation. This equation, which also satisfies the equilibrium, is the bi-harmonic equation

7

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^4\varphi =0.\end{array}$$where the unknown scalar field $\varphi$, represents Airy’s stress function.

Then, the stress equations in terms of Airy’s stress function for plane stress and plane strain are

8

$$\begin{array}{cc}\hfill {\sigma }_{11}& =\frac{{\partial }^2\varphi }{\partial y^2},{\sigma }_{12}=-\frac{{\partial }^2\varphi }{\partial x\partial y},\hfill \end{array}$$

9

$$\begin{array}{cc}\hfill {\sigma }_{22}& =\frac{{\partial }^2\varphi }{\partial x^2},{\sigma }_{23}={\sigma }_{31}=0.\hfill \end{array}$$The analysis presented hereafter primarily focuses on the plane stress condition, under which ${\sigma }_{33}=0$.

Following Silverman’s (1964) proposal of a general way to obtain the Airy’s stress functions, when applied to cantilever beams subjected to the same loadings as here, the Airy’s stress function $\varphi (x,y)$ is assumed to have the same polynomial behavior as the tractions that describe the boundary loads, becoming

10

$$\begin{array}{c}\hfill \varphi (x,y)=\sum _{k=0}^{\mathcal{N}}{x}^k{\chi }_k^{⟨\alpha ⟩}(y),\end{array}$$where ${\chi }_k^{⟨\alpha ⟩}(y)$ are functions to be determined for each layer indicated by the label superscript1$⟨\alpha ⟩\in [1,2,3]$.

Silverman [31] stipulated for normal loads that $\mathcal{N}=\mathcal{M}+2$ and for shear loads $\mathcal{N}=\mathcal{M}+1$, where $\mathcal{M}$ is the polynomial order of the loading, taking into account that a moment will be applied at an end and a linear distributed load at the top, so $\mathcal{N}=3$, thus

11

$$\begin{array}{c}\hfill \varphi (x,y)=\sum _{k=0}^{\mathcal{N}=3}{x}^k{\chi }_k(y)=x^3{\chi }_3(y)+x^2{\chi }_2(y)+x{\chi }_1(y)+{\chi }_0(y),\end{array}$$and some of the above equations will be replaced to build the equations needed to get the beam solution.

Equation 11 will be replaced, over and over, in the elasticity equations: first into the stress-stress function relation given by Eq. 8, turning into

12

$$\begin{array}{c}\hfill {\sigma }_{11}(x,y)=x^3\frac{{\partial }^2{\chi }_3}{\partial y^2}+x^2\frac{{\partial }^2{\chi }_2}{\partial y^2}+x\frac{{\partial }^2{\chi }_1}{\partial y^2}+\frac{{\partial }^2{\chi }_0}{\partial y^2},\end{array}$$

13

$$\begin{array}{c}\hfill {\sigma }_{22}(x,y)=2\left(3,x,{\chi }_3,+,{\chi }_2\right),\end{array}$$and

14

$$\begin{array}{c}\hfill {\sigma }_{12}(x,y)=-3x^2\frac{\partial {\chi }_3}{\partial y}-2x\frac{\partial {\chi }_2}{\partial y}-\frac{\partial {\chi }_1}{\partial y};\end{array}$$second, into the strains given by the stress–strain relations, Eq. 3, yielding

15

$$\begin{array}{c}\hfill {\epsilon }_{11}=\frac{1}{E}\left[x^3,\frac{{\partial }^2{\chi }_3}{\partial y^2},+,x^2,\frac{{\partial }^2{\chi }_2}{\partial y^2},+,x,\frac{{\partial }^2{\chi }_1}{\partial y^2},+,\frac{{\partial }^2{\chi }_0}{\partial y^2},-,2,\nu ,,\left(3,x,{\chi }_3,+,{\chi }_2\right)\right],\end{array}$$

16

$$\begin{array}{c}\hfill {\epsilon }_{22}=-\frac{1}{E}\left[\nu ,,\left(x^3,\frac{{\partial }^2{\chi }_3}{\partial y^2},+,x^2,\frac{{\partial }^2{\chi }_2}{\partial y^2},+,x,\frac{{\partial }^2{\chi }_1}{\partial y^2},+,\frac{{\partial }^2{\chi }_0}{\partial y^2}\right),-,2,\left(3,x,{\chi }_3,+,{\chi }_2\right)\right],\end{array}$$and

17

$$\begin{array}{c}\hfill {\epsilon }_{12}=-\frac{1+\nu }{E}\left[3,x^2,\frac{\partial {\chi }_3}{\partial y},+,2,x,\frac{\partial {\chi }_2}{\partial y},+,\frac{\partial {\chi }_1}{\partial y}\right];\end{array}$$and last, in the sense of strain compatibility given by Eq. 5, leading to

18

$$\begin{array}{c}\hfill x^3\frac{{\partial }^4{\chi }_3}{\partial y^4}+x^2\frac{{\partial }^4{\chi }_2}{\partial y^4}+x\left(12,\frac{{\partial }^2{\chi }_3}{\partial y^2},+,\frac{{\partial }^4{\chi }_1}{\partial y^4}\right)+4\frac{{\partial }^2{\chi }_2}{\partial y^2}+\frac{{\partial }^4{\chi }_0}{\partial y^4}=0.\end{array}$$Since the above equation must be valid for all x, collecting the coefficients related to this variable and the unrelated one and then integrating each coefficient until obtaining the derivative-free form of ${\chi }_k^{⟨\alpha ⟩}(y)$, is possible to get

19a

$$\begin{array}{cc}& {\chi }_3^{⟨\alpha ⟩}\left(y\right)=A_1^{⟨\alpha ⟩}{y}^3+A_2^{⟨\alpha ⟩}{y}^2+A_3^{⟨\alpha ⟩}y+A_4^{⟨\alpha ⟩},\hfill \end{array}$$

19b

$$\begin{array}{cc}& {\chi }_2^{⟨\alpha ⟩}\left(y\right)=A_5^{⟨\alpha ⟩}{y}^3+A_6^{⟨\alpha ⟩}{y}^2+A_7^{⟨\alpha ⟩}y+A_8^{⟨\alpha ⟩},\hfill \end{array}$$

19c

$$\begin{array}{cc}& {\chi }_1^{⟨\alpha ⟩}\left(y\right)=-\frac{3}{5}{A}_1^{⟨\alpha ⟩}{y}^5-A_2^{⟨\alpha ⟩}{y}^4+A_9^{⟨\alpha ⟩}{y}^3+A_{10}^{⟨\alpha ⟩}{y}^2+A_{11}^{⟨\alpha ⟩}y+A_{12}^{⟨\alpha ⟩},\hfill \end{array}$$

19d

$$\begin{array}{cc}& {\chi }_0^{⟨\alpha ⟩}\left(y\right)=-\frac{1}{5}{A}_5^{⟨\alpha ⟩}{y}^5-\frac{1}{3}{A}_6^{⟨\alpha ⟩}{y}^4+A_{13}^{⟨\alpha ⟩}{y}^3+A_{14}^{⟨\alpha ⟩}{y}^2+A_{15}^{⟨\alpha ⟩}y+A_{16}^{⟨\alpha ⟩},\hfill \end{array}$$ where $A_j^{⟨\alpha ⟩}$ is the integral constants to be determined; until now, there are 16 constants for each layer.

Returning to $\varphi (x,y)$, given by Eq. 11, and introducing Eq. 19, it becomes

20

$$\begin{array}{cc}& \varphi (x,y)=-\frac{1}{5}\left(3,x,A_1,+,A_5\right)y^5-\left(x,A_2,+,\frac{1}{3},A_6\right)y^4+\left(x^3,A_1,+,x^2,A_5,+,x,A_9,+,A_{13}\right)y^3+\left(x^3,A_2,+,x^2,A_6,+,x,A_{10},+,A_{14}\right)y^2\hfill \\ \hfill & +\left(x^3,A_3,+,x^2,A_7,+,x,A_{11},+,A_{15}\right)y+x^3{A}_4+x^2{A}_8+x{A}_{12}+A_{16}.\hfill \end{array}$$The fulfilment of the bi-harmonic equation, given by Eq. 7, can be easily proved and thus omitted here.

Replacing Eq. 20 in Eqs. 12–14, the stress equations develop into

21

$$\begin{array}{cc}& {\sigma }_{11}^{⟨\alpha ⟩}=2\left(x^3,\left(3,y,A_1^{⟨\alpha ⟩},+,A_2^{⟨\alpha ⟩}\right),+,x^2,\left(3,y,A_5^{⟨\alpha ⟩},+,A_6^{⟨\alpha ⟩}\right)\right)\hfill \\ \hfill & \left(+,x,\left\{3,y,\left[A_9^{⟨\alpha ⟩},-,2,y,\left(y,A_1^{⟨\alpha ⟩},+,A_2^{⟨\alpha ⟩}\right)\right],+,A_{10}^{⟨\alpha ⟩}\right\}\right)\hfill \\ \hfill & \left(-,2,y^2,\left(y,A_5^{⟨\alpha ⟩},+,A_6^{⟨\alpha ⟩}\right),+,3,y,A_{13}^{⟨\alpha ⟩},+,A_{14}^{⟨\alpha ⟩}\right),\hfill \end{array}$$

22

$$\begin{array}{cc}& {\sigma }_{22}^{⟨\alpha ⟩}=2\left(3,x,y^3,A_1^{⟨\alpha ⟩},+,3,x,y^2,A_2^{⟨\alpha ⟩},+,3,x,y,A_3^{⟨\alpha ⟩},+,3,x,A_4^{⟨\alpha ⟩},+,y^3,A_5^{⟨\alpha ⟩}\right)\hfill \\ \hfill & \left(+,y^2,A_6^{⟨\alpha ⟩},+,y,A_7^{⟨\alpha ⟩},+,A_8^{⟨\alpha ⟩}\right),\hfill \end{array}$$and

23

$$\begin{array}{cc}& {\sigma }_{12}^{⟨\alpha ⟩}=-3y^2\left(3,x^2,A_1^{⟨\alpha ⟩},+,2,x,A_5^{⟨\alpha ⟩},+,A_9^{⟨\alpha ⟩}\right)-2y\left(3,x^2,A_2^{⟨\alpha ⟩},+,2,x,A_6^{⟨\alpha ⟩}\right)\hfill \\ \hfill & \left(+,A_{10}^{⟨\alpha ⟩}\right)-3x^2{A}_3^{⟨\alpha ⟩}-2x{A}_7^{⟨\alpha ⟩}+3y^4{A}_1^{⟨\alpha ⟩}+4y^3{A}_2^{⟨\alpha ⟩}-A_{11}^{⟨\alpha ⟩},\hfill \end{array}$$being the stress equations’ final form with the coefficients $A_j^{⟨\alpha ⟩}$ to be determined.

The displacement components are derived by integrating strain–displacement relations. However, due to space constraints, the detailed derivation is provided in 6. It’s worth noting that $A_j^{⟨\alpha ⟩}$ for all $j\in \left[12,,,15,,,16\right]$ does not appear in either the displacement or stress equations. Equation 20 shows that these integral constants are associated with the linear and constant terms. However, the procedure shown in 6, adds three more constants: $A_j^{⟨\alpha ⟩}$ for $j\in \left[17,,,18,,,19\right]$, bringing the total of non-null constants per layer to 16; for three layers, this yields 48 constants to be determined.

So far, no considerations beyond the degree of Airy’s stress function have been made, thus it can be employed for other configurations aside from the cantilever.

Determination of the integral constants for a three-layer cantilever beam

The boundary and continuity conditions between the layers must be applied to determine the integral constants, following the three-layer cantilever beam shown in Fig. 1, considering a narrow rectangular cross-section of unit width.

The upper surface is subject to the linear distributed load, thus

24

$$\begin{array}{c}\hfill {\sigma }_{22}^{⟨1⟩}=p_1+p_{2}x,{\sigma }_{12}^{⟨1⟩}=0,\text{at}y=-h^{⟨1⟩},\end{array}$$where $p_1$ and $p_2$ are the constant and linear parameters, respectively;

and on the lower one, it is free, so

25

$$\begin{array}{c}\hfill {\sigma }_{22}^{⟨3⟩}=0,{\sigma }_{12}^{⟨3⟩}=0,\text{at}y=h^{⟨2⟩}+h^{⟨3⟩}.\end{array}$$The continuity conditions between the layers will be applied on the interface between the layer $⟨1⟩$ and $⟨2⟩$, therefore

26

$$\begin{array}{c}\hfill \begin{array}{c}{\sigma }_{22}^{⟨1⟩}-{\sigma }_{22}^{⟨2⟩}=0,\\ {\sigma }_{12}^{⟨1⟩}-{\sigma }_{12}^{⟨2⟩}=0,\\ u_1^{⟨1⟩}-u_1^{⟨2⟩}=0,u_2^{⟨1⟩}-u_2^{⟨2⟩}=0,\end{array}\text{at}y=0;\end{array}$$and on the interface between the layer $⟨2⟩$ and $⟨3⟩$,

27

$$\begin{array}{c}\hfill \begin{array}{c}{\sigma }_{22}^{⟨2⟩}-{\sigma }_{22}^{⟨3⟩}=0,\\ {\sigma }_{12}^{⟨2⟩}-{\sigma }_{12}^{⟨3⟩}=0,\\ u_1^{⟨2⟩}-u_1^{⟨3⟩}=0,u_2^{⟨2⟩}-u_2^{⟨3⟩}=0,\end{array}\text{at}y=h^{⟨2⟩},\end{array}$$The boundary condition at the end $x=0$ are applied in a weakly way using

28

$$\begin{array}{cc}& {\int }_{-h^{⟨1⟩}}^0{\sigma }_{11}^{⟨1⟩}\text{d}y+{\int }_0^{h^{⟨2⟩}}{\sigma }_{11}^{⟨2⟩}\text{d}y+{\int }_{h^{⟨2⟩}}^{h^{⟨2⟩}+h^{⟨3⟩}}{\sigma }_{11}^{⟨3⟩}\text{d}y=0,\hfill \end{array}$$

29

$$\begin{array}{cc}& {\int }_{-h^{⟨1⟩}}^0{\sigma }_{12}^{⟨1⟩}\text{d}y+{\int }_0^{h^{⟨2⟩}}{\sigma }_{12}^{⟨2⟩}\text{d}y+{\int }_{h^{⟨2⟩}}^{h^{⟨2⟩}+h^{⟨3⟩}}{\sigma }_{12}^{⟨3⟩}\text{d}y-Q=0,\hfill \end{array}$$

30

$$\begin{array}{cc}& {\int }_{-h^{⟨1⟩}}^0{\sigma }_{11}^{⟨1⟩}y\text{d}y+{\int }_0^{h^{⟨2⟩}}{\sigma }_{11}^{⟨2⟩}y\text{d}y+{\int }_{h^{⟨2⟩}}^{h^{⟨2⟩}+h^{⟨3⟩}}{\sigma }_{11}^{⟨3⟩}y\text{d}y-M=0.\hfill \end{array}$$Displacement conditions at the fixed end are applied at a point within the lower interface, thus

31

$$\begin{array}{cc}& u^{⟨2⟩}=u^{⟨3⟩}=0,v^{⟨2⟩}=v^{⟨3⟩}=0,\frac{\partial v^{⟨2⟩}}{\partial x}=\frac{\partial v^{⟨3⟩}}{\partial x}=0,\hfill \\ \hfill & \text{at}x=l\text{and}y=h^{⟨2⟩}.\hfill \end{array}$$The equations which, following the substitutions, remain functions of x and/or y are dismembered, as they must be valid for any value of x and/or y. Thus, an equation of the type $\mathcal{A}{x}^2+\mathcal{B}x+\mathcal{C}y+\mathcal{D}=0$ is divided into four equations: $\mathcal{A}=0$, $\mathcal{B}=0$, $\mathcal{C}=0$ and $\mathcal{D}=0$, for all real x and y.

The constants for each layer can be determined by solving a system of 48 equations with 48 unknowns. These constants are incorporated into the code as the matrix C, which can be found in the computeC function within the files SandwichSolution.jl and sandwich_solution.f90 for Julia and FORTRAN, respectively. The matrix C has dimensions [3, 16], where 3 represents the number of layers and 16 represents the number of constants per layer.

Numerical XFEM verification

To demonstrate a numerical application of the analytical solution obtained in this work, the XFEM was employed to model the body and the material interface. Therefore, it is assumed that the reader already has a background in the fundamentals of the method. More details on the specific implementation of the XFEM used in this section can be found in [8].

Extended finite element method

The XFEM is an extension of the FEM, developed to deal more effectively with problems involving discontinuities, such as fractures, cracks, and material interfaces. The main idea of the method is to utilize the concept of partition of unity (PU) to augment the approximation space with non-polynomial terms without necessitating a conformal mesh, as in FEM. Moreover, since non-smooth effects are usually localized phenomena, the enrichment is solely applied in the vicinity of the non-smoothness. For instance, this allows the geometric transfer of a material interface, previously represented by the edges of an element in the FEM, to the approximation space as an enrichment.

A consequence of localized enrichment is the appearance of partially enriched elements in the intermediate region between fully enriched and standard FEM elements, referred to as blending elements. The presence of a discontinuous gradient in the primal field can lead to pathological oscillations due to the lack of reproducibility of the approximation space in these particular elements since the PU property is not met, resulting in an ill-conditioned stiffness matrix and influencing the convergence rate [14]. Therefore, a proposal to improve the limitation is put forth by [7], which employs hierarchical shape functions in the blending elements to compensate for the loss of reproducibility of the approximation space.

Therefore, the approximated XFEM displacement field can be written as

32

$$\begin{array}{cc}& {\mathit{u}}^h(\mathit{x})=\underset{\text{Standard FEM}}{\underbrace{{\sum }_{i\in I}{N}_i^{\mathrm{std}}(\mathit{x}){\mathit{u}}_i}}+\underset{\text{Enrichment}}{\underbrace{{\sum }_{i\in I^{\star }}{N}_i^{\star }(\mathit{x})\left(|,{\gamma }^h,(\mathit{x})|-|\gamma ,({\mathit{x}}_i),|\right){\mathit{a}}_i}}\hfill \\ \hfill & +\underset{\text{Hierarchical Blending}}{\underbrace{\sum _{i\in J}{\widehat{N}}_i(\mathit{x}){\mathit{b}}_i}},\hfill \end{array}$$with $N_i^{\text{std}}(\mathit{x})$ and ${\mathit{u}}_i$ represent the shape functions and the degree of freedom (DOF) associated with the FEM approximation, respectively. The set $I$ includes all nodes. $N_i^{\text{enr}}(\mathit{x})$ and ${\mathit{a}}_i$ represent the shape functions and the DOF associated with the enrichment, respectively. The subset $I^{\star }\subset I$ consists of nodes that are enriched. The set $J$ includes hierarchical nodes associated with the sides of the elements, connecting an enriched node to a non-enriched node. Notably, there are no changes in the shape functions in the first two terms due to the incorporation of the hierarchical corrector term.

The PU employed in this work was assumed to be the Q4 shape functions, which are low-order when compared to the fifth-order displacements of the analytical solution. Therefore, the hierarchical quadratic functions corresponding to the mid-side nodes, which could be added between the nodes 1–2, 2–3, 3–4, and 4–1, are respectively represented by

33a

$$\begin{array}{cc}& {\widehat{N}}_5\left(\xi ,,,\eta \right)=\frac{1}{2}\left(1,-,{\xi }^2\right)\left(1,-,\eta \right),\hfill \end{array}$$

33b

$$\begin{array}{cc}& {\widehat{N}}_6\left(\xi ,,,\eta \right)=\frac{1}{2}\left(1,-,{\eta }^2\right)\left(1,+,\xi \right),\hfill \end{array}$$

33c

$$\begin{array}{cc}& {\widehat{N}}_7\left(\xi ,,,\eta \right)=\frac{1}{2}\left(1,-,{\xi }^2\right)\left(1,+,\eta \right),\hfill \end{array}$$

33d

$$\begin{array}{cc}& {\widehat{N}}_8\left(\xi ,,,\eta \right)=\frac{1}{2}\left(1,-,{\eta }^2\right)\left(1,-,\xi \right).\hfill \end{array}$$

Since the focus of this work is on the analytical solution, the XFEM was numerically integrated using the Gauss–Legendre technique, as in FEM, to provide the slice plots to analyze the fields. Moreover, in the convergence analysis, sub-cell integration in the enriched elements was also employed. Therefore, unless otherwise stated, sixteen $(4\times 4)$ Gauss-Legendre integration points are used for $S_{\mathit{GL}}^{\star }$. For the sub-cell strategy $S_{\mathit{ubC}}^{\star }$, the region on either side of the interface and the element are partitioned into four triangles, each having three integration points, resulting in a total of 24 integration points per element.

Before proceeding, the $L^2$-displacement and energy norms are defined prior to the assessment of numerical results. The $L^2$ norm for the displacement is calculated using

34

$$\begin{array}{c}\hfill {∥\mathit{u}∥}_{L^2}^2={\displaystyle {\int }_{\mathrm{\Omega }}}\mathit{u}·\mathit{u}\text{d}\mathrm{\Omega },\end{array}$$and the respective error is

35

$$\begin{array}{cc}& {∥\mathrm{\Delta }∥}_{L^2}^2={∥\mathit{u},-,{\mathit{u}}^h∥}_{L^2}^2\hfill \\ \hfill & =\underset{\mathrm{\Omega }}{\int }\left(\mathit{u},-,{\mathit{u}}^h\right)·\left(\mathit{u},-,{\mathit{u}}^h\right)d\mathrm{\Omega },\hfill \end{array}$$where $\mathit{u}$ is the analytical solution displacement and ${\mathit{u}}^h$ is the approximated displacement using the different XFEM strategies evaluated in this work. Furthermore, the relative $L^2$ norm is

36

$$\begin{array}{c}\hfill {\eta }_{L^2}=\frac{{∥\mathit{u},-,{\mathit{u}}^h∥}_{L^2}^2}{{∥\mathit{u}∥}_{L^2}^2}.\end{array}$$The energy norm is calculated as follows

37

$$\begin{array}{c}\hfill {∥\mathit{u}∥}_E^2={\int }_{\mathrm{\Omega }}{\mathit{\sigma }}^T\mathbf{E}\mathit{\sigma }d\mathrm{\Omega },\end{array}$$and the respective error is

38

$$\begin{array}{cc}& {∥\mathrm{\Delta }∥}_E^2={∥\mathit{u},-,{\mathit{u}}^h∥}_E^2\hfill \\ \hfill & ={\displaystyle {\int }_{\mathrm{\Omega }}}{\left(\mathit{\sigma },-,{\mathit{\sigma }}^h\right)}^T\mathbf{E}\left(\mathit{\sigma },-,{\mathit{\sigma }}^h\right)\text{d}\mathrm{\Omega },\hfill \end{array}$$where $\mathit{\sigma }$ is the analytical stress and ${\mathit{\sigma }}^h$ is the approximated stress solution obtained through the evaluation of various XFEM strategies. The term ${\mathit{\sigma }}^h$ is gathered from integration points and the expressions are numerically integrated. Moreover, the relative energy norm is

39

$$\begin{array}{c}\hfill {\eta }_E=\frac{{∥\mathit{u},-,{\mathit{u}}^h∥}_E^2}{{∥\mathit{u}∥}_E^2}.\end{array}$$

Numerical cases

Two scenarios with reversed combinations of material properties are analyzed for this section, along with a particular case of the general solution consisting of only two layers, outlined in 7. As a result, the following sub-numbering is utilized to differentiate between these scenarios:

• Benchmark 3.1: Representing a classic sandwich beam featuring two rigid external skins with $E^{⟨1⟩}=E^{⟨3⟩}=500$, ${\nu }^{⟨1⟩}={\nu }^{⟨3⟩}=0.3$, and a compliant core with $E^{⟨2⟩}=100$ and ${\nu }^{⟨2⟩}=0.25$;

• Benchmark 3.2: Illustrating a coating structure comprising two outer films with $E^{⟨1⟩}=E^{⟨3⟩}=100$, ${\nu }^{⟨1⟩}={\nu }^{⟨3⟩}=0.25$, and a rigid core with $E^{⟨2⟩}=500$ and ${\nu }^{⟨2⟩}=0.3$;

• Benchmark 3.3: Depicting a structure with two distinct layers: an upper layer with $E^{⟨1⟩}=100$, ${\nu }^{⟨1⟩}=0.25$, and a bottom layer with $E^{⟨2⟩}=500$ and ${\nu }^{⟨2⟩}=0.3$.

For Benchmarks 3.1 and 3.2, a transverse load of $Q=-1·{10}^{-2}$ and a bending moment of $M=-5·{10}^{-2}$ were applied at the end section. The upper part of the beam had $p_1=3$ and $p_2=2/5$, resulting in $p_1+p_{2}x=3$ when $x=0$ and $p_1+p_{2}x=5$ when $x=l$. The dimensions and elastic properties of layers $⟨1⟩$ and $⟨3⟩$ are the same. The cantilever has a length of $l=5$ and a height of one unit. Layers $⟨1⟩$ and $⟨3⟩$ have a height of $h^{⟨1⟩}=h^{⟨3⟩}=0.25$, while layer $⟨2⟩$ has a height of $h^{⟨2⟩}=0.5$. Furthermore, in Benchmark 3.3, the load values were kept, but the distributed loading was removed. The beam maintained its dimensions, but the layers heights were equated: $h^{⟨1⟩}=h^{⟨2⟩}$.

The Cartesian meshes are discretized and refined linearly, with sizes of $a\times b$ with a varying in the set $\{30,50,70,90,110,270\}$ and b varying in the set $\{06,10,14,18,22,54\}$. The interfaces are positioned at the center of the elements. Furthermore, the coarser $30\times 06$ and finer $270\times 54$ meshes, including the interfaces and integration sub-cells, are depicted in Fig. 2.

[See PDF for image]

Fig. 2

Benchmarks 3.1 and 3.2—Meshes, integration sub-cells and material interfaces: a$30\times 06$ for Gauss-Legendre integration; b$30\times 06$ for sub-cell integration; c$270\times 54$ for Gauss–Legendre integration; d$270\times 54$ for sub-cell integration

According to [9], the incorporation of Neumann boundary conditions is a simple process in XFEM and is carried out in the same manner as in FEM. However, the same authors also note that in order to impose smooth Dirichlet boundary conditions in XFEM, it is only necessary to prescribe the $u_i$ and set $a_i$ to zero. However, for non-smooth behavior, such as those found at the fixed end of the cantilever beam, the values for the enriched DOFs are often unknown. The same authors highlight that it is possible to solve a local interpolation problem to determine the values of $u_i$ and $a_i$, or to impose the Dirichlet conditions weakly, for example, through the use of penalties.

In this work, as the analytical solution of the problem was previously known, the boundary displacement components $u_x$ and $u_y$ are respectively illustrated in Fig. 3a and b for the free end, and Fig. 3c and d for the fixed end. The boundary condition at the fixed end displays kinks2 in $u_x$, while both components behave smoothly at the free end. In light of this, a viable approach could be to apply displacements at the free end and traction at the fixed end, thereby eliminating the need for any special treatment within these boundaries.

[See PDF for image]

Fig. 3

Boundary displacement components: a$u_x$, b$u_y$ at free end; and c$u_x$, d$u_y$ at fixed end

Therefore, in terms of numerical modeling, the boundary conditions are set as follows: equivalent displacement was applied at the free end, $x=0$, corresponding to the loading conditions at this boundary, originally given by the integral Eqs. 28–30; traction was applied at the upper part of the beam, at $y=-h^{⟨1⟩}$. Similarly, an equivalent traction is applied at the fixed end, $x=l$, where the support is located.

Benchmark 3.1: classic three-layer cantilever beam

Convergence analysis

A convergence analysis of the classic sandwich beam is conducted by considering two strategies in terms of the relative error norms. The results for the relative errors in the displacement and energy norms, respectively, as h decreases in the elements, can be seen in Fig. 4.

[See PDF for image]

Fig. 4

Benchmark 3.1: Convergence behavior in terms of the relative error norms a$L^2$ norm and b energy norm

Overall, both norms exhibited linear convergence on a logarithmic scale, demonstrating minimal sensitivity to the type of integration employed. Regarding convergence rates, the type of integration made no significant difference. The relative $L^2$ norm had an average convergence rate of approximately $m=1.9$, while the relative energy norm displayed an average convergence rate of $m=0.93$.

Analysis of displacement and stress fields, along with $L^2$ and energy norms per element

Continuing the investigation, a qualitative analysis will be conducted on the fields of interest. Therefore, the $70\times 14$ mesh will be examined and visualized. This particular mesh allows for the discernment and analysis of the fields in each element type. Thus, the results3 for the stress field ${\sigma }_{\mathit{xx}}$ are presented in Fig. 5, demonstrating higher magnitudes of this stress in the outer layers, particularly near the fixed support and the upper and lower surfaces of the beam, as expected for this particular problem.

[See PDF for image]

Fig. 5

Benchmark 3.1: ${\sigma }_{\mathit{xx}}$ stress field for a $70\times 14$ mesh

In the stress ${\sigma }_{\mathit{yy}}$ presented in Fig. 6, negative stresses are generally observed on the upper surface due to the distributed loading on top of the beam. As the beam approaches the fixed end, the stress field becomes locally positive. Furthermore, a variation in this stress field is observed near the fixed end between the blending and standard elements. This can be attributed to the imposition of boundary conditions at this interface. It is worth noting that the behavior of the stress ${\sigma }_{\mathit{yy}}$ will be further addressed in 3.6, where cross-sections of the beam will be examined in more detail.

[See PDF for image]

Fig. 6

Benchmark 3.1: ${\sigma }_{\mathit{yy}}$ stress field for a $70\times 14$ mesh

In addition, the stress field ${\sigma }_{\mathit{xy}}$ shown in Fig. 7 demonstrates a parabolic pattern along the y-axis in the domain. Furthermore, there is a noticeable variation in the behavior of this stress near the fixed end. This can be explained by Saint-Venant’s principle, which states that stress fields near the boundaries are influenced by the local traction distribution and should decrease with distance from them. However, the exact solution presented here does not take these effects into account, which could result in a more accurate representation of the field in this region.

[See PDF for image]

Fig. 7

Benchmark 3.1: ${\sigma }_{\mathit{xy}}$ stress field for a $70\times 14$ mesh

For the relative $L^2$ element-wise norm, the largest error was most pronounced near the fixed end and did not exhibit significant errors near the material interfaces. Additionally, this error was symmetric with respect to the neutral axis of the beam. In contrast, the relative energy element-wise norm showed greater errors in regions where the stress field was less smooth, corresponding to areas with the highest stress amplitude (Figs. 8, 9).

[See PDF for image]

Fig. 8

Benchmark 3.1: Relative $L^2$ norm errors per element $70\times 14$ mesh

[See PDF for image]

Fig. 9

Benchmark 3.1: Relative energy norm errors per element for a $70\times 14$ mesh

The stress behavior for a refined mesh of size $270\times 54$ is depicted in Fig. 10, along with the corresponding relative norms displayed in Fig. 11. These figures closely approximate the exact solution, featuring smooth contours similar to it. However, the stress ${\sigma }_{\mathit{yy}}$ exhibits artificial abrupt changes near both the interfaces and the fixed support due to the approximate nature of the XFEM implementation and the low-order shape functions employed. In contrast, these abrupt changes do not occur in the exact solution, which is smooth. Additionally, the relative $L^2$ norm consistently shows higher errors in elements closer to the fixed support, while the relative energy norm indicates significantly greater errors concentrated in the enriched elements within the fixed support region.

[See PDF for image]

Fig. 10

Benchmark 3.1: a${\sigma }_{\mathit{xx}}$, b${\sigma }_{\mathit{yy}}$ and (c) ${\sigma }_{\mathit{xy}}$ stress fields for a $270\times 54$ mesh

[See PDF for image]

Fig. 11

Benchmark 3.1: Relative error in the a$L^2$ norm and b energy norm for a $270\times 54$ mesh

Stress fields at the integration points

Furthermore, an analysis of the evolution of approximated fields is done in y-direction where data was extracted from the integration points from a slice nearly to the fixed end and nearly to the center in a mesh of size $110\times 22$; these results are illustrated in Figs. 12 and 13, respectively. To facilitate the analysis, vertical gray lines have been incorporated to signify the positions of the nodes linked with the blending and enriched elements.

Examining the stress profiles near the fixed end, at the integration points located at $x=4.98243$. The behavior of the longitudinal displacement $u_x$, transversal displacement $u_y$ and stress field ${\sigma }_{\mathit{xx}}$, respectively, in Fig. 12a–c, closely matches the exact solutions.

Regarding the stress field ${\sigma }_{\mathit{yy}}$, depicted in Figure 12(d), the overall behavior aligns with the analytical solution. However, there are minor, in comparison to the magnitudes of the other component stresses, jagged fluctuations, which can be attributed to the use of Q4 shape functions, as these are low-order interpolators for this problem. Also, it is possible to observe that the hierarchical blending elements exhibit lower absolute values of fluctuations and a better approximation compared to the rest of the enriched elements. Moreover, in Fig. 12e the stress ${\sigma }_{\mathit{xy}}$ shows lower fluctuations near the interface due to the hierarchical blending elements.

Additionally, a slice was taken near the beam’s center, and the results were depicted in Fig. 13. Concerning the displacement $u_x$, a uniform behavior consistent with the expected response for this beam and its loading can be observed, unlike the region near the fixed end, which exhibits a more complex behavior due to the prescribed boundary condition.

[See PDF for image]

Fig. 12

Benchmark 3.1: Profile of the displacements $u_x$, $u_y$, stresses ${\sigma }_{\mathit{xx}}$, ${\sigma }_{\mathit{yy}}$, ${\sigma }_{\mathit{xy}}$, nearly to the fixed end, respectively in (a), (b), (c), (d), and (e)

[See PDF for image]

Fig. 13

Benchmark 3.1: Profile of the displacements $u_x$, $u_y$, stresses ${\sigma }_{\mathit{xx}}$, ${\sigma }_{\mathit{yy}}$, ${\sigma }_{\mathit{xy}}$, nearly to the center, respectively in (a), (b), (c), (d), and (e)

Benchmark 3.2: three-layer coating beam

Convergence analysis

In accordance with the methodology outlined in 3.3, a convergence analysis of the coating beam is carried out by considering the two strategies in terms of the relative error norms. Figure 14 presents the results, showcasing the variations in relative errors in $L^2$ and energy norms as the element size h decreases.

[See PDF for image]

Fig. 14

Benchmark 3.2: Convergence behavior in terms of the relative error in the a$L^2$ norm and in b energy norm

In this benchmark, the behavior of the relative norms was the same as in Benchmark 3.1, also with the type of integration having no influences. The relative $L^2$ norm, shown in Fig. 14a, had a mean convergence rate of $m=1.99$, while the relative energy norm, depicted in Fig. 14b, exhibited a mean convergence rate of $m=0.95$, showing small differences from the previous benchmark.

Analysis of displacement and stress fields, along with $L^2$ and energy norms per element

Similar to the qualitative analysis conducted for Benchmark 3.1, an analysis will now be performed for the coating case.

The result for the stress field ${\sigma }_{\mathit{xx}}$ are shown in Fig. 15, revealing stress concentrations at the material interface, near the fixed support. This concentration is expected to occur in the core, which has a higher stiffness.

[See PDF for image]

Fig. 15

Benchmark 3.2: ${\sigma }_{\mathit{xx}}$ stress field for a $70\times 14$ mesh

In the stress field ${\sigma }_{\mathit{yy}}$ shown in Fig. 16, negative stresses are observed on the upper surface due to the distributed loading on top of the beam. As the beam approaches the fixed end, the stress field becomes locally positive, as demonstrated in 3.5. Additionally, a variation in this stress field can be observed near the fixed end, between the blending and standard elements, as well as in the more rigid layer. This can be attributed to the boundary condition being imposed at this end and its interaction with enriched DOFs (Figs. 17, 18, 19).

[See PDF for image]

Fig. 16

Benchmark 3.2: ${\sigma }_{\mathit{yy}}$ stress field for a $70\times 14$ mesh

The same qualitative observations regarding the stress field ${\sigma }_{\mathit{xy}}$, $L^2$, and energy norms as those in Benchmark 3.1 can be applied here.

[See PDF for image]

Fig. 17

Benchmark 3.2: ${\sigma }_{\mathit{xy}}$ stress field for a $70\times 14$ mesh

[See PDF for image]

Fig. 18

Benchmark 3.2: Relative $L^2$ norm errors per element for a $70\times 14$ mesh

[See PDF for image]

Fig. 19

Benchmark 3.2: Relative energy norm errors per element for a $70\times 14$ mesh

In addition, the stresses and relative norms for a refined mesh of $270\times 54$ size are illustrated in Figs. 20 and 21. The observations made in Benchmark 3.1 regarding the results for the finer mesh apply here as well. Again, these figures closely approximate the exact solution, featuring smooth contours similar to it. However, the stress ${\sigma }_{\mathit{yy}}$ exhibits artificial abrupt changes near both the interfaces and the fixed support due to the approximate nature of the XFEM implementation and the low-order shape function employed. On the other hand, this does occur in the exact solution, which is smooth.

[See PDF for image]

Fig. 20

Benchmark 3.2: a${\sigma }_{\mathit{xx}}$, b${\sigma }_{\mathit{yy}}$ and ${\sigma }_{\mathit{xy}}$ stress fields for a $270\times 54$ mesh

[See PDF for image]

Fig. 21

Benchmark 3.2: Relative error in the a$L^2$ norm and b energy norm for a $270\times 54$ mesh

Displacement and stress fields at the integration points

The examination of the cross-section will begin by scrutinizing the stress profiles in the vicinity of the fixed end, specifically focusing on the integration points located at $x=4.98243$.

The behavior of the longitudinal displacement $u_x$, as depicted in Fig. 22a, was consistent with the exact solution. Additionally, the transversal displacement $u_y$, shown in Fig. 22b, exhibited absolute values that were slightly higher than those of the exact solution.

Furthermore, the behavior of ${\sigma }_{\mathit{xx}}$, illustrated in Fig. 22c, was found to be identical to the exact solution. The stress field ${\sigma }_{\mathit{yy}}$, shown in Fig. 22d, generally followed the analytical solution. However, small jagged fluctuations around the exact solution were observed due to the use of low-order approximation functions.

Regarding the stress ${\sigma }_{\mathit{xy}}$, as seen in Fig. 22e, the stress field did not replicate the steep behavior of the core, as observed in the exact solution. Instead, it appeared to be shifted toward lower values compared to the exact solution. However, lower fluctuations near the interface were also evident in the hierarchical elements.

Additionally, a slice was taken at the beam’s center, and the profiles were redrawn in Fig. 23. The behavior in $u_x$ was linear, indicating the absence of support imposition effects observed in the previous slice. Conversely, the behavior in $u_y$ was more aligned with the exact solution. The improvement in the approximation in this region of the beam is also evident in ${\sigma }_{\mathit{xy}}$, in Fig. 23e, where stress was more accurately interpolated than in the previous cross-section, clearly highlighting the benefits of hierarchical blending elements.

[See PDF for image]

Fig. 22

Benchmark 3.2: Profile of the displacements $u_x$, $u_y$, stresses ${\sigma }_{\mathit{xx}}$, ${\sigma }_{\mathit{yy}}$, ${\sigma }_{\mathit{xy}}$, nearly to the fixed end, respectively in (a), (b), (c), (d), and (e)

[See PDF for image]

Fig. 23

Benchmark 3.2: Profile of the displacements $u_x$, $u_y$, stresses ${\sigma }_{\mathit{xx}}$, ${\sigma }_{\mathit{yy}}$, ${\sigma }_{\mathit{xy}}$, nearly to the center, respectively in (a), (b), (c), (d), and (e)

Benchmark 3.3: Two-layer cantilever beam

Stress fields at the integration points

To conclude the section, the analysis now turns to a simplified scenario with only two layers, utilizing the solution presented in 7. The resulting stress fields, as depicted in Figure 24, show that the stress field ${\sigma }_{\mathit{xx}}$ agrees with the exact response. However, ${\sigma }_{\mathit{yy}}$ still exhibits small jagged oscillations, albeit around the zero value, since there is no distributed loading. This indicates that while the approximate solutions are close to the exact ones, they still need refinement. Nevertheless, as emphasized throughout this work, there is a significant reduction in oscillations in the blending elements.

[See PDF for image]

Fig. 24

Benchmark 3.3: ${\sigma }_{\mathit{xx}}$, ${\sigma }_{\mathit{yy}}$ and ${\sigma }_{\mathit{xy}}$ profiles, nearly to the center, respectively in (a), (b), and (c)

Conclusions

This study presented a new two-dimensional analytical solution for a sandwich beam, which was employed to verify an XFEM implementation. The proposed solution, involving fifth-degree polynomials and material interfaces, demonstrated significant value for both engineering applications and the verification of numerical methods. The XFEM, applied with low-order shape functions to capture the nontrivial solution field, proved effective, although the results indicate that higher-order approximations would further enhance stress field continuity across elements. The analytical solution has been implemented and made publicly available as an open-source library in Julia and FORTRAN, providing a valuable resource for benchmarking and verification. The analytical solution is derived under the assumption that each layer consists of an isotropic and homogenous material. Expanding the current formulation to include sandwich beams with orthotropic or transversely isotropic material layers is recognized as a promising direction for future investigation.

Acknowledgements

The author Rodrigo Rossi wishes to acknowledge the support of CNPq, National Council for Scientific and Technological Development of Brazil, grant number 309430/2021-6.

Author Contributions

The authors contributed equally to this work and approved it for publication.

Data Availability

The implementation of the analytical solution derived in this study is open-source and available in our GitHub repository at github.com/rodrgz/2D3LayerCantileverBeam. It is provided in both the Julia and FORTRAN languages. The code is licensed under the Mozilla Public License, version 2.0. Additionally, a README file is included to guide users through the process of coupling the modules in each language and running the provided examples.

Declarations

Conflict of interest

The authors declare that there is no Conflict of interest.

1

The label superscript $⟨\alpha ⟩$, which is not an exponent, will be propagated in several quantities in the subsequent developments and, to simplify the notation, will be omitted, being shown only in the first definition and in the last; thus, the quantities already defined will be referred to for each layer as: ${\mathit{\sigma }}^{⟨\alpha ⟩}=\mathit{\sigma }$; ${\mathit{\epsilon }}^{⟨\alpha ⟩}=\mathit{\epsilon }$; $E_{\alpha }=E$; ${\nu }_{\alpha }=\nu$. The superscript also appears when there are distinctions in the equations for each layer.

2

This could be attributed to the inherent behavior of analytical solutions in proximity to displacement constraints, where in the analytical solution the absence of displacement is imposed at a single point in the material interface, through Eq. 31. However, it is important to note that this may not perfectly replicate the expected real-world behavior of the beam support.

3

The post-processing method treats values within elements intersected by the interface as conventional finite elements, thus not visually depicting the XFEM field. Stress values are extrapolated from integration points to nodes using least squares with bi-linear shape functions. Nodal values from each element sharing the node are averaged to obtain the plotted values. A dedicated subsection will examine magnitudes at the integration points for a more comprehensive analysis.

Santosh Kapuria

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Abrate S (1998) Impact on composite structures, 1st edn. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511574504

2. Amandeep, S; Padhee, SS. Analytic solution of timoshenko-like deformation in bidirectional functionally graded beams. J Eng Mech; 2024; 150, 2 04023118. [DOI: https://dx.doi.org/10.1061/JENMDT.EMENG-7317]

3. Bardella, L. Explicit analytical solutions for the full plane-stress field in sandwich beams under flexure governed by zigzag warping. Compos Struct; 2024; 329, [DOI: https://dx.doi.org/10.1016/j.compstruct.2023.117754] 117754.

4. Beltrami, E. Osservazioni sulla nota precedente. Rendiconti Acc dei Lincei; 1892; 1, 5 pp. 141-143.

5. Carrera, E. Historical review of zig-zag theories for multilayered plates and shells. Appl Mech Rev; 2003; 56, 3 pp. 287-308. [DOI: https://dx.doi.org/10.1115/1.1557614]

6. Cheng, S; Wei, X; Jiang, T. Stress distribution and deformation of adhesive-bonded laminated composite beams. J Eng Mech; 1989; 115, 6 pp. 1150-1162. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-9399(1989)115:6(1150)]

7. Chessa, J; Wang, H; Belytschko, T. On the construction of blending elements for local partition of unity enriched finite elements. Int J Numer Meth Eng; 2003; 57, 7 pp. 1015-1038. [DOI: https://dx.doi.org/10.1002/nme.777]

8. Da Rosa, RE; Rossi, R. Assessment of eqp in xfem for weak discontinuities. J Braz Soc Mech Sci Eng; 2023; 45, 6 312. [DOI: https://dx.doi.org/10.1007/s40430-023-04211-z]

9. Fries, TP; Belytschko, T. The extended/generalized finite element method: an overview of the method and its applications: The gefm/xfem: An overview of the method. Int J Numer Meth Eng; 2010; 84, 3 pp. 253-304. [DOI: https://dx.doi.org/10.1002/nme.2914]

10. Gerstner, RW. Stresses in a composite cantilever. J Compos Mater; 1968; 2, 4 pp. 498-501. [DOI: https://dx.doi.org/10.1177/002199836800200410]

11. Ghugal, YM; Shimpi, RP. A review of refined shear deformation theories for isotropic and anisotropic laminated beams. J Reinf Plast Compos; 2001; 20, 3 pp. 255-272. [DOI: https://dx.doi.org/10.1177/073168401772678283]

12. Hoff, NJ; Mautner, SE. Bending and buckling of sandwich beams. J Aeron Sci; 1948; 15, 12 pp. 707-720. [DOI: https://dx.doi.org/10.2514/8.11699]

13. Dj, H; Hj, D; Wq, C. Analytical solution for functionally graded anisotropic cantilever beam under thermal and uniformly distributed load. J Zhejiang Univ-Sci A; 2007; 8, 9 pp. 1351-1355. [DOI: https://dx.doi.org/10.1631/jzus.2007.A1351]

14. Khoei, AR. Extended finite element method: theory and applications; 2014; 1 Wiley: [DOI: https://dx.doi.org/10.1002/9781118869673]

15. Krajcinovic, D. Sandwich beam analysis. J Appl Mech; 1972; 39, 3 pp. 773-778. [DOI: https://dx.doi.org/10.1115/1.3422787]

16. Krajcinovic, D. Sandwich beams with arbitrary boundary conditions. J Eng Ind; 1975; 97, 3 pp. 873-880. [DOI: https://dx.doi.org/10.1115/1.3438695]

17. Lekhnitskii, SG. Strength calculation of composite beams. Vestnik Inzhen I Teknikov; 1935; 9, pp. 137-148.

18. Lekhnitskii, SG. Anisotropic plates; 1968; New York, Gordon and Breach:

19. Mindlin, RD. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. J Appl Mech; 1951; 18, 1 pp. 31-38. [DOI: https://dx.doi.org/10.1115/1.4010217]

20. Muskhelishvili NI (1977) Some basic problems of the mathematical theory of elasticity. Springer, Dordrecht. https://doi.org/10.1007/978-94-017-3034-1

21. Naik, NS; Sayyad, AS. 2d analysis of laminated composite and sandwich plates using a new fifth-order plate theory. Latin Am J Solids Struct; 2018; [DOI: https://dx.doi.org/10.1590/1679-78254834]

22. Pagani, A; Yan, Y; Carrera, E. Exact solutions for static analysis of laminated, box and sandwich beams by refined layer-wise theory. Compos B Eng; 2017; 131, pp. 62-75. [DOI: https://dx.doi.org/10.1016/j.compositesb.2017.08.001]

23. Pagano, N. Exact solutions for composite laminates in cylindrical bending. J Compos Mater; 1969; 3, 3 pp. 398-411. [DOI: https://dx.doi.org/10.1177/002199836900300304]

24. Pagano, N. Influence of shear coupling in cylindrical. Bending of anisotropic laminates. J Compos Mater; 1970; 4, 3 pp. 330-343. [DOI: https://dx.doi.org/10.1177/002199837000400305]

25. Rao, KM; Ghosh, BG. Exact analysis of unsymmetric laminated beam. J Struct Div; 1979; 105, 11 pp. 2313-2325. [DOI: https://dx.doi.org/10.1061/JSDEAG.0005279]

26. Rieder G (1960) Topologische fragen in der theorie der spannungsfuktionen. Braunschweigischen Wissenschaftlichen Gesellschaft 12:4–65. https://doi.org/10.24355/DBBS.084-201301101242-0

27. Sayyad, AS; Ghugal, YM. Bending, buckling and free vibration of laminated composite and sandwich beams: a critical review of literature. Compos Struct; 2017; 171, pp. 486-504. [DOI: https://dx.doi.org/10.1016/j.compstruct.2017.03.053]

28. Schile, RD. A nonhomogeneous beam in plane stress. J Appl Mech; 1962; 29, 3 pp. 582-583. [DOI: https://dx.doi.org/10.1115/1.3640609]

29. Sidhardh, S; Ray, MC. Exact solution for size-dependent elastic response in laminated beams considering generalized first strain gradient elasticity. Compos Struct; 2018; 204, pp. 31-42. [DOI: https://dx.doi.org/10.1016/j.compstruct.2018.07.030]

30. Sierakowski, R; Ebcioglu, I. On interlaminar shear stresses in composites. J Compos Mater; 1970; 4, 1 pp. 144-149. [DOI: https://dx.doi.org/10.1177/002199837000400117]

31. Silverman, IK. Orthotropic beams under polynomial loads. J Eng Mech Div; 1964; 90, 5 pp. 293-319. [DOI: https://dx.doi.org/10.1061/JMCEA3.0000540]

32. Sun, SL; Li, XF. Higher-order effect of axial force on free vibration and buckling of functionally graded sandwich beams. ZAMM - J Appl Math Mech; 2024; 104, 74777604 [DOI: https://dx.doi.org/10.1002/zamm.202301054] e202301054.

33. Timoshenko S, Goodier JN (1952) Theory of elasticity. Engineering societies monographs, McGraw-Hill Series in Mechanical Engineering. https://www.cambridge.org/core/product/identifier/S036839310012471X/type/journal_article

34. Wang, M; Liu, Y. Analytical solution for bi-material beam with graded intermediate layer. Compos Struct; 2010; 92, 10 pp. 2358-2368. [DOI: https://dx.doi.org/10.1016/j.compstruct.2010.03.013]

35. Yu, YY. A new theory of elastic sandwich plates-one-dimensional case. J Appl Mech; 1959; 26, 3 pp. 415-421.111250 [DOI: https://dx.doi.org/10.1115/1.4012054]

36. Cx, Z; Liu, Y. Plane elasticity solutions for beams with fixed ends. J Zhejiang Univ-Sci A; 2015; 16, 10 pp. 805-819. [DOI: https://dx.doi.org/10.1631/jzus.A1500043]

37. Zhong, Z; Yu, T. Analytical solution of a cantilever functionally graded beam. Compos Sci Technol; 2007; 67, 3–4 pp. 481-488. [DOI: https://dx.doi.org/10.1016/j.compscitech.2006.08.023]