Mesh-adapted stress analysis of multilayered plates using a layerwise model

This paper proposes a new finite-element modelling of a recent layerwise model for multilayered plates. This layerwise model is built from a specific 3D stress-field expansion along the thickness direction and involves, in particular, interlaminar transverse shear and out-of-plane stresses as generalized stresses. Its main feature is that 3D equilibrium equations and free-edge boundary conditions are directly taken into account into the stress-based construction of the model. A dual displacement-based finite-element discretization is implemented using the FEniCS software package and a remeshing strategy is proposed based on a novel error indicator. The error indicator is built based on the 3D stress field directly deduced from the layerwise generalized stresses and compared to a reconstructed stress field based on the model generalized displacements. The proposed error indicator is shown to identify the most critical parts of a laminate structure associated with complex 3D stress fields such as boundaries or stress concentration/singularity regions (near free-edges or delamination fronts). Through the combination of thickness discretization and in-plane mesh refinement in regions of interest, the proposed framework therefore offers an attractive alternative to 3D solid finite elements for an accurate prediction of stress states in composite laminates.


Introduction
Multilayered plates have very interesting mechanical properties that make them widely used in aerospace, automotive, telecommunication structures and civil engineering.A multilayered plate is represented as a pile of homogenized anisotropic plies made of fiberreinforced composites.However, the highly anisotropic and heterogeneous nature of such laminates, the prediction of their overall properties is a challenging task.One of the major issues in design and analysis of such plates is related to free-edge effects.It has been proved that the differences in the elastic properties of adjacent layers generally result in a highly concentrated interlaminar stresses near free edges [1][2][3][4][5].Many models were derived to accurately capture these free-edge effects.Highly detailed three-dimensional (3D) finite-element models are computationally expensive and will only result in accurate stress predictions for sufficiently refined meshes since they rely on displacement interpolations.Two-dimensional plate models have therefore been introduced in order to simplify these computations while trying to keep a sufficiently accurate description of local 3D stress fields.
Equivalent single layer (ESL) models represent the laminate as an equivalent homogeneous plate.Many ESL models based on higher order theories have been proposed in the literature [6][7][8][9][10][11][12][13][14] and are usually derived using two main approaches: asymptotic approaches and axiomatic approaches.The first class derives the plate model from the full 3D formulation of the problem, assuming the thickness of the plate goes to zero and using asymptotic expansion in which the leading order leads to Kirchhoff-Love plate theory [15].The second approach is based on assuming a priori 3D fields, and the plate theory is derived by integration through the thickness and variational tools [16][17][18].Although ESL models can provide acceptable results for the laminate global response, they may lead to very inaccurate estimations of local response especially near free-edges.
Layerwise models, in which each layer is considered as an independent plate, have therefore been proposed to improve the local stress representation [19][20][21][22][23][24].Layerwise models have been proved to a very good alternative to 3D models since interpolation choices along the z-direction take into account the specificities of the laminate.The interested reader can refer to [25,26] for a general overview of such models and to [27,28] for a recent comparison between ESL, zig-zag and layerwise models.
Following the ideas of Pagano's model [29], a layerwise model named LS1 was developed in [30][31][32][33][34][35][36][37][38][39][40].In this model, the laminate is considered as a superposition of Reissner-Mindlin plates linked together by interfacial stresses which are considered as additional generalized stresses.The main difference between LS1 models and other models is that LS1 is a stressbased approach, while other models are either displacement or mixed stress/displacement approaches.
However, the LS1 model presents some conceptual drawbacks since, for instance, 3D stress-free boundary conditions cannot be exactly fulfilled.Second, the model is derived by means of the Hellinger-Reissner mixed variational principle, so that there is no theoretical guarantee of the convergence to the 3D model as the number of mathematical layers per physical layer increases.Generalizing upon the same ideas, a layerwise model called statically compatible (SCLS1), was introduced in [41] in which the divergence of the interlaminar transverse shear stresses is introduced as an additional generalized stress.Doing so, the SCLS1 model produces 3D stress field satisfying the local 3D balance equations and boundary conditions provided that their 2D plate counterparts are satisfied.The model can therefore be derived by means of the complementary potential energy minimum principle ensuring the convergence of its refined version to the exact 3D model as the number of mathematical layers per physical layer increases.Aiming at providing an operational tool for stress analysis in multilayered plates, this paper is concerned with the development of a mesh adaptation strategy based on an error indicator built from the local 3D stress field and a reconstructed 3D displacement field.
The paper is organized as follows: the SCLS1 model equations and finite-element implementation are discussed in section .Section is dedicated the reconstruction of the 3D displacement and to error indicator computation used in the mesh adaptation.Finally, section illustrate the method efficiency in capturing regions of interest in various configurations.

The SCLS1 model
In this section, the equations of the SCLS1 model for elastic multilayered plates are recalled.This model is derived from the 3D continuum equations by considering Statically Compatible Layerwise Stresses with first-order membrane stress approximations per layer in the thickness direction.The generalized stresses of the proposed model are actually those of a Reissner-Mindlin plate per layer in addition to the inter-laminar shear and normal stresses at the interfaces between layers and the divergences of these inter-laminar shear stresses.The plate kinematics is then obtained through duality arguments.

Problem description and notations
We consider a linear elastic multilayered plate composed of n monoclinic elastic layers.The plates occupies the 3D domain where ω ⊂ R 2 is the middle surface of the plate.In the following, x and y are the in-plane coordinates and z is the out-of-plane coordinate.The following notations are introduced: • The superscript i and j, j + 1 indicates layer i and the interface between the layer j and j + 1 with 1 ≤ i ≤ n and 1 ≤ j ≤ n − 1, respectively.By extension, the superscript 0,1 refers to the lower face ω − = ω × {h − 1 } and the superscript n, n + 1 refers to the upper face ω + = ω × {h + n }. • (S i = (S i klmn )) is the fourth-order 3D compliance tensor of layer i with the minor and major symmetries: S i klmn = S i lkmn = S i klnm = S i mnkl and it is positive definite.Its inverse is the 3D elasticity stiffness tensor and is denoted by (C i klmn ) for layer i.The tensor (C i klmn ) possesses the same symmetries as (S i klmn ) and it is also positive definite.• S i is monoclinic in direction z : S i αβγ 3 = S i α333 = 0 • σ αβ (x, y, z) are the in-plane stress components, σ α3 (x, y, z) are the transverse shear stresses and σ 33 (x, y, z) is the normal stress.• αβ (x, y, z) are the in-plane strain components, α3 (x, y, z) are the transverse strain stresses and 33 (x, y, z) is the normal strain.• u α (x, y, z) are the in-plane 3D displacement components, u 3 (x, y, z) is the normal 3D displacement component.
The plate is loaded on its upper face ω + and lower face ω − with the distributed surface forces T + = (T + k ) and T − = (T − k ), respectively.The lateral boundary is decomposed into two complementary parts: a free part ) is set to zero, and a restrained part where the displacement u = (u k ) is set to zero.Here, the subsets ∂ω T and ∂ω u are the partition of ∂ω, and n = (n k ) is the outer normal to ∂ω T .

Equations of the 3D model
The 3D elastic problem is to find in a statically compatible stress field σ = (σ kl ), a kinematically strain field = ( kl ) which comply with the constitutive equation: where a stress field σ is said to be statically compatible if it complies with the equilibrium equations: and the stress conditions on the lower and the upper faces: and on the lateral boundary: A strain field is kinematically compatible if there exists a displacement field u = (u k ) complying with the displacement conditions on the lateral boundary: and such that : The static of SCLS1 model The SCLS1 model assumes the following form of the 3D stresses in layer i: where P i k , k = 0, 1, 2, 3, are the orthogonal Legendre-like polynomial basis defined on layer i by: for and where N i αβ , M i αβ and Q i α are, respectively the classical membrane forces, bending moments and shear forces in layer i, and τ j,j+1 α and ν j,j+1 are the transverse shear and out-of-plane normal stresses at the interface j, j + 1. π j,j+1 is an additional variable whose interpretation will appear later.

Equilibrium equations
The σ 3D stress field will comply with 3D equilibrium equations ( 2), if and only if, the following equations hold true for all (x, y) in ω and for all i = 1, . . ., n and j = 0, . . ., n : The last equation gives the interpretation of π j,j+1 which is equal to the divergence of the interlaminar shear stress vector τ j,j+1 = τ j,j+1 α .Now stress boundary conditions also have to be enforced in addition to Eq. (11).The lateral boundary conditions σ 3D ij n j = 0 on ∂ T are equivalent to the following equations for i = 1, . . ., n and j = 0, . . ., n : And the boundary conditions (3) on the upper an the lower faces write, respectively, and It should be noticed that the boundary conditions ( 12) and ( 13) cannot be simultaneously verified unless T ± α n α = 0 on ∂ω T , which will be assumed in the sequel.Moreover, from the last equations of (11) for j = 0 and j = n, we see that: Finally, the stress field σ 3D is statically compatible when it complies with the generalized equilibrium equations on ω: (11) for i = 1, . . ., n and j = 1, . . ., n − 1, ( 13) and ( 14), and with the generalized stress free boundary conditions on ∂ω T : (12) for i = 1, . . ., n and j = 1, . . ., n − 1.

Generalized displacements and strains
The SCLS1 generalized displacements are U i α (x, y), U i 3 (x, y), i α (x, y) and V j,j+1 (x, y), respectively the two in-plane displacements, the vertical displacement, the two bending rotations, and V j,j+1 is a kinematical variable having the dimension of an area.The following expressions for i = 1, . . ., n and j = 1, . . ., n − 1 give the relation between the generalized displacements and the 3D displacement field. and, The generalized strains dual of the generalized stresses , ν j,j+1 , π j,j+1 for i = 1, . . ., n and j = 1, . . ., n − 1 are respectively expressed in terms of the generalized displacements as:

The SCLS1 model constitutive equations
The constitutive equations of the SCLS1 model are derived using the stress energy associated to σ 3D .They are given by : for 1 ≤ i ≤ n and for 1 • Bending constitutive equations of layer i: • Transverse shear constitutive equation of layer i: • Shear constitutive equation of interface j, j + 1:

Finite element discretization
A finite-element discretization of the SCLS1 model has been proposed in [41] using the MPFEAP in-house software described in [36].In our numerical study, the SCLS1 multilayered plate model has been implemented in the open-source finite element FEniCS package [42,43].The FEniCs Project is a collection of free and open-source software components with the common goal to enable automated solutions of differential equations.The components provide scientific computing tools for working with computational meshes, finite element variational formulations of ordinary and partial differential equations, and numerical linear algebra.We therefore benefit from FEniCS high-level domain-specific language for implementing the variational formulation associated with the SCLS1 model.Building upon the FEniCS implementation of a Reissner-Mindlin plate model [44], we define a generalized function space for the SCLS1 generalized displacement degrees of freedom.More precisely, the retained discretization is based on a mesh of triangular elements with quadratic interpolation for all kinematical variables.As is the case for classical FE discretization of Reissner-Mindlin plate models, FE discretization of the SCLS1 model leads to shear-locking in the thin plate limit.Selective reduced integration is then used on the shear part of the strain [44].

Comparison with other layerwise models
In [41], the SCLS1 model has been compared with the LS1 model and reference 3D computations.This work showed that the SCLS1 model is as accurate as the LS1 model and is even closer to refined 3D solutions near free edges since it can correctly satisfy stressfree boundary conditions.Besides, an intensive comparison between the LS1 model and other layerwise models derived of the Carrera Unified Formulation (CUF) family has been performed in [45].The main conclusion of this work was that the LS1 model exhibits a similar accuracy to LM4 (mixed fourth-order) and LD3 (displacement third-order) layerwise models.This conclusion therefore also holds for the SCLS1 model considered here.Moreover, in contrast to these models, LS1 and SCLS1 exhibit much fewer degrees of freedom per node.For instance for a laminate with n = 4, LS1 has 20 (5n) dofs/node, SCLS1 has 23 (6n − 1) whereas LD3 has 39 and LM4 has 102.One important feature of LS1 and SCLS1 is that no assumption is made on the displacement variations through the thickness but rather on the stress.Therefore, obtaining a complete 3D displacement field must be performed by a post-processing procedure which we will now describe.

Mesh adaptivity based on field reconstructions
Although being much cheaper than LM4 or LD3 CUF models, the SCLS1 model is still quite expensive due to its high number of degrees of freedom per node.It can be seen as specific, mechanically-based, discretization in the z direction and can therefore be compared to a 3D discretization with a more accurate representation of the stress fields in the z direction.It becomes therefore beneficial to optimize the in-plane mesh for improved computational efficiency.The purpose of this section is to fulfil this goal by building an error indicator for mesh adaptation.
We propose to define this indicator as follows: from the finite-element computed generalized displacements (U i α , U i 3 , i α , V j,j+1 ) fields in the (x, y)-plane, we first aim at reconstructing a 3D displacement field u(x, y, z).We then derive the associated 3D strain and stresses using the local constitutive equation.The so-obtained reconstructed stress field σ is then compared to the initial 3D stress σ 3D obtained from the generalized stresses , ν j,j+1 , π j,j+1 ) via Eqs.( 7)- (9).See the illustration of the scheme in Fig. 1.

Field reconstructions
In this subsection, we propose to reconstruct u by considering a continuous piecewise linear variation of its components u i along the z direction.This interpolation will have to be as close as possible to satisfying Eqs. ( 15)- (19).Let us mention that we tried other interpolations (in particular of higher-order) or reconstruction strategies but the latter gave the most satisfying results.
Let us first consider the in-plane displacement field.We first build an auxiliary in-plane displacement u d α , with α = 1, 2, as follows: Fig. 1 The reconstruction scheme u d α is piecewise-linear and complies with Eqs. ( 15), ( 16) but is not continuous at the interfaces.To achieve our initial goal, u α is obtained by performing an L2-projection of u d α over piecewise linear continuous functions over z.Now, we aim to find the reconstructed out-of-plane displacement u 3 as a continuous piecewise linear function of z which is compatible with the generalized displacements U i 3 and V j,j+1 in the sense of the following equations: Introducing a continuous piecewise linear interpolation for u 3 (x, y, z) of the following form: where φ j,j+1 are linear shape functions and q j,j+1 are the corresponding nodal values, the above equations become: where t [q] = q 0,1 , . . ., q n,n+1 , [B] is a matrix of dimension (2n − 1, n + 1) and [F 3 ] a vector of dimension (2n − 1).The solution to the above problem is computed in the least-squares sense and gives a direct characterization of the degrees of freedom [q] as a function of the generalized displacements U i 3 and V i,i+1 .Finally, from the previously reconstructed 3D displacement field u i , the strain field is computed using the 3D compatibility equations and then the reconstructed stress tensor σ is computed using the 3D constitutive equations.

Error indicator and mesh adaptation
The error indicator which will be used for mesh adaptation is then computed based on the difference between σ 3D and σ in terms of elastic energy.This error is computed for each triangular element: where e(σ) = 1 2 σ : S : σ and e denotes a given element e.Each mesh element is then ordered in a decreasing fashion based on its error indicator value: where N is the total number of elements.We then tag the first K elements which contribute to at least a fraction η of the total error The tagged elements are then automatically refined by FEniCS mesh adaptation procedures.

Illustrative applications
In this section, we investigate different illustrative applications assessing the quality of the stress field approximation, error indicator and mesh refinement strategy.The last examples consider more practical situations arising when designing composite laminates, namely stress concentrations near holes with associated free-edge singularities and interfacial stress singularities in the presence of interface delamination.

Homogeneous laminate
This first example considers a homogeneous square plate of length l = 1 and thickness h = 0.2.The constitutive material is assumed to be isotropic with E = 10 GPa and ν = 0.3.
The plate is fully clamped on its boundary and subject to a uniform loading of intensity q = 8.Calculations are performed considering a uniform discretization of n = 3 and n = 5 layers across the thickness and have been compared to finite-element computations using 3D solid elements on a very fine mesh.The initial mesh was a structured mesh with two triangular elements on each side of the square plate.First the case with n = 3 layers is considered.As expected, the multilayered plate solution is of very good quality near the plate center after only one refinement step as shown in Fig. 2.
We therefore investigate the quality of the computed stress field at a point of coordinate (x = 0.01, y = 0.5) near the left edge.In Fig. 3a, we compare the multilayered stress field σ 3D with its reconstruction as described in section at the same point near the edge.It can be observed that the reconstruction does not agree with σ 3D for the initial coarse mesh, indicating that mesh size should be refined in this region.Figure 3b, c illustrate the evolution of σ 3D and σ near the border when refining the mesh.It can be seen that mesh refinement provides a much better agreement between both stress fields.The error indicator therefore correctly identifies regions located near the clamped boundaries as the most critical regions as evidenced by the final mesh layout of Fig. 4a obtained after 6 refinement steps.
Performing the same comparison in the case when the plate thickness is discretized in n = 5 layers shows the same behaviour (Fig. 5).Although σ 3D and σ are a little closer for the initial coarse mesh, the deviation is still significant indicating that in-plane mesh resolution is not fine enough.The mesh refinement procedure yields a similar final mesh layout, with fine cells concentrated along the borders (see Fig. 4b), and better agreement between σ 3D and σ at the final stage.On both Figs.3c and 5c, the reference solid FE solution is also represented, showing a good agreement with the multilayered stress field after mesh refinement.The effect of mesh refinement is further illustrated when plotting the evolution of the total relative error indicator in Fig. 6, defined as: It can be seen that the relative errors decrease when refining the mesh and tend to stabilize after a few iterations only.Besides, relative errors are larger for n = 3 than n = 5 which may indicate that the mesh reconstruction if of higher quality for n = 5 layers than n = 3 layers.Let us point out that the value obtained for such errors cannot be considered neither as a guaranteed level of error with respect to an exact solution nor as an upper bound to the true error.It is however an error indicator, as showed by the previous results, which can be used qualitatively to assess the solution accuracy.

Triple laminate
The second example considers a heterogeneous square plate of length l = 1 and total thickness h = 0.2 made of a triple laminate consisting of a central core of thickness e 2 = 0.12 and two symmetric skins of thickness e 1 = 0.04 each.The constitutive materials are assumed to be isotropic with E = 50 GPa and ν = 0.2 for the skins and E = 10 GPa and ν = 0.3 for the core.Loading and conditions are the same as for the homogeneous plate.Calculations are performed considering a discretization consisting of one mathematical layer in both skins and in the core (total of n = 3 layers) and a discretization consisting of one mathematical layer per skin and 3 layers for uniformly discretizing the core thickness (total of n = 5 layers).Again the multilayered plate model computations have been compared to reference 3D solid finite-element computations on a fine mesh.First, we considered the case with n = 3 layers.As before, the solution is of lesser quality near the supports and stress fields are therefore compared at the same (x = 0.01, y = 0.5) location as before.Figure 7 shows the comparison of the stress field σ 3D with its reconstruction across the plate thickness for various mesh refinement steps.It can be observed that σ 3D and its reconstruction do not match for the initial coarse mesh, indicating the mesh size should be refined in this region.Mesh adaptation improves the quality of the solution in such regions as evidenced by the good agreement with the reference 3D FE solution.Performing the same comparison using a more refined discretization with n = 5 layers in the thickness exhibits a similar behaviour (Fig. 8).Although σ 3D and σ are a little closer for the initial coarse mesh, the deviation is still significant indicating that in-plane mesh resolution is not fine enough.A similar refined mesh layout is obtained with fine cells concentrated along the borders (see Fig. 9), and better agreement between σ 3D and σ at the final stage.
Finally, the evolution of the total relative error indicator as a function of mesh refinement steps in Fig. 10 exhibits a similar behaviour as for the homogeneous plate.

Laminate with a circular hole
The third example considers a rectangular multilayered plate of length l = 6, width w = 1 and total thickness h = 0.01.The plate is perforated by a circular hole of radius R = 0.15 in its center (Fig. 11a).The laminate is made of a transversely isotropic material of elastic properties E T = 14.48 GPa,E L = 137.9GPa,ν T = 0.21,ν L = 0.21,μ T = 5.86 GPa and μ L = 5.86 GPa with L (resp.T ) denoting the fiber longitudinal direction (resp.the perpendicular transverse direction).The laminate consists of 6 plies (one layer per ply) with fibers oriented at [0 • , 90 • , 45 • , −45 • , 90 • , 0 • ] with respect to the horizontal direction.A tensile loading is applied to the plate through an imposed horizontal displacement U i = ±Ue x for all plies i = 1, . . ., 6.
Applying the proposed reconstruction and error estimation on this problem yields a globally more refined mesh with finer regions located near the top and bottom boundaries of the circular hole, see Fig. 11b obtained after 4 refinement steps.More insight can also be gained at visualizing the individual layer contributions to the total error.For instance, Fig. 12 plots the contribution of the 45 • (layer 3) and −45 • (layer 4) layers to the total error.These two contributions are the most dominant one as regards stress concentrations near the hole.The effect of the material anisotropy on these two contributions can also be clearly observed.

Double-Cantilever Beam with delaminated interface
The final example we consider is that of a rectangular multilayered plate of the same dimensions as before (without the circular hole) and the same lamination properties.We model a portion of a delaminated interface located in the middle interface ((i, i+1) = (3, 4)) in the region x ≤ 1 by forcing the interface stresses ν 3,4 and τ 3,4  α to be zero on this region.This results in an appropriate modification of the constitutive equations of the SCLS1 model and the corresponding finite-element implementation.
The plate is clamped on its right boundary, and positive (resp.negative) vertical displacements U i 3 = +U (resp.U i 3 = −U) are enforced on the left part for the top layers i = 4, 5, 6 (resp.bottom layers i = 1, 2, 3), simulating a Double-Cantilever Beam test (see Fig. 13).
As expected, the mesh adaptation procedure mainly concentrates the finer cells near the delamination front at which interface stresses are the most singular, see Fig. 14.The proposed procedure can therefore be considered to be coupled with a delamination Fig. 12 Error indicator maps in layers 3 and 4 (top and middle) and total error for all layers (bottom) on the initial mesh Fig. 13 The initial mesh for the DCB problem propagation model for which stresses driving the delamination front propagation will be well resolved.

Conclusions and perspectives
In this paper, a statically compatible layerwise stress model for laminated plates (SCLS1) is considered for an accurate representation of 3D elastic fields.A mesh adaptation strategy is then developed which relies on the reconstruction of 3D displacement fields from the model generalized displacements, the error indicator being obtained by a constitutive error between both fields.The obtained results indicate that: • the error indicator is able to refine the mesh in regions with complex 3D stress fields • these critical regions indeed correspond to plate edges, notches or delamination fronts The proposed methodology can be further improved by pointing out that refined layerwise models such as the one considered here is appropriate in critical regions near boundaries, free-edges, delaminated interfaces, etc.This point is indeed properly identified by the proposed remeshing procedure.In the bulk region away from these critical zones, it would we sufficient to adopt an equivalent single-layer plate model based on a Love-Kirchhoff kinematics for instance.Although the remeshing procedure favours coarse cells in such regions, mitigating the number of unnecessary degrees of freedom, an additional gain could then be obtained by mixing a layerwise model for critical regions with an equivalent single-layer model for the remaining part.
A second potential line of work is concerned with the fact that, although the layerwise model is built at the continuous level from a stress-based perspective complying with the balance equations, its numerical resolution is performed through a displacement-based approximation for the in-plane variations.As a consequence, the resulting generalized stress fields, and therefore, the associated 3D stress field, do not satisfy strongly the balance equations.In order to maintain the initial philosophy of a stress-based statically compatible construction, developing a stress-based finite-element discretization of the model would be an interesting approach, potentially paving the way to obtaining more rigorous error estimators than the one considered here.

••
Normal constitutive equation of interface j, j + 1: Constitutive equation for the π generalized stress at interface j, j + 1:

Fig. 2 Fig. 3
Fig.2Energy densities across the plate thickness computed for σ 3D and σ at the plate center (n = 3)

Fig. 4 b c 5
Fig.4 Final refined meshes for the homogeneous plate for different thickness discretization levels

Fig. 6
Fig.6 Total relative error evolution for 3 and 5 layers discretizations

Fig. 7
Fig.7 Energy densities across the plate thickness computed for σ 3D and σ at the plate edge for n = 3

Fig. 8 Fig. 9
Fig. 8 Energy densities across the plate thickness computed for σ 3D , σ and σ ref at the plate edge for n = 5