On the coupling of local 3D solutions and global 2D shell theory in structural mechanics

Most of mechanical systems and complex structures exhibit plate and shell components. Therefore, 2D simulation, based on plate and shell theory, appears as an appealing choice in structural analysis as it allows reducing the computational complexity. Nevertheless, this 2D framework fails for capturing rich physics compromising the usual hypotheses considered when deriving standard plate and shell theories. To circumvent, or at least alleviate this issue, authors proposed in their former works an in-plane-out-of-plane separated representation able to capture rich 3D behaviors while keeping the computational complexity of 2D simulations. However, that procedure it was revealed to be too intrusive for being introduced into existing commercial softwares. Moreover, experience indicated that such enriched descriptions are only compulsory locally, in some regions or structure components. In the present paper we propose an enrichment procedure able to address 3D local behaviors, preserving the direct minimally-invasive coupling with existing plate and shell discretizations. The proposed strategy will be extended to inelastic behaviors and structural dynamics.


Introduction
Many mechanical systems and complex structures involve plate and shell parts or components whose main particularity is having a characteristic dimension (the one related to the thickness) much lower than the other ones (in-plane dimensions). In order to analyse such structures, beam, plate and shell theories have been developed in solid mechanics [1,2] and extended later to many other physics, like flows in narrow gaps, thermal or electromagnetic problems in laminates, among many others. In these theories the introduction of appropriate kinematic and mechanic hypotheses on the evolution of the solution through the thickness of the plate (shell) allows the reduction of the general 3D mechanical problem to a 2D one involving the in-plane coordinates.
However, in many cases, when addressing complex coupled physics, inelastic behaviors or any other exhibiting localization, the validity of hypotheses able to reduce models from 3D to 2D becomes doubtful and consequently in order to ensure accurate results 3D discretizations seem compulsory. However mesh-based solutions of models defined in such degenerated domains is a challenging issue because the resulting meshes usually involve too many degrees of freedom, where the mesh size is almost determined by the domain thickness and the material and/or solution details to be represented. In order to alleviate the associated computational complexity in [3] authors proposed computing the fully 3D solution employing an in-plane-out-of-plane separated representation whose computational complexity remains the one characteristic of 2D plate or shell simulations, using the proper generalized decomposition-PGD-method [4].
However the proposed in-plane-out-of-plane separated representation appears to be too intrusive to be implemented in structural mechanics commercial software that generally propose different plate and shell finite elements, even in the case of multilayered composites plates or shells. For this reason in this work we propose a minimally-intrusive method which allows integrating fully 3D local descriptions in plate or shell models implemented in any software, without affecting its computational complexity that remains the one related to standard 2D analyses. For that purpose, two different enrichment routes will be considered, the first based on the use of that separated representation and the second on a simple condensation. The former methodology allows for very fine 3D enriched descriptions while the last is particularly well adapted to address inelastic and dynamical behaviors.

Elastostatic problem definition
We consider the linear elastostatic problem defined in the plate domain depicted in Fig. 1 where C is the Hooke's fourth order tensor. The relation between strain ε and displacement u (with components u = (u, v, w)) writes where G = ∇ s • = 1 2 (∇ • +∇ T •) is the symmetric gradient operator. Considering an homogeneous and isotropic material and using the Voigt notation, the Hooke's tensor can be written as In absence of volumetric body forces, the displacement field evolution u(x) for x ∈ is described by the linear momentum balance equation The domain boundary ∂ is partitioned into Dirichlet, D , and Neumann, N , boundaries, where displacement u g and tractions T are enforced respectively.
The problem weak form associated to the strong form (4) lies in looking for the displacement field u verifying the Dirichlet boundary conditions such that the weak form fulfills for any test function u * , with the trial and test fields defined in appropriate functional spaces.
In this type of domains, plate theory is usually used in order to reduce the general 3D mechanical problem to a 2D one involving the in-plane coordinates only. Two kinds of theories exist: the thin plate theory proposed by Kirchoff [5] which establishes that the normal remains straight and orthogonal to the middle plane after deformation and the thick plate theory proposed by Reissner [6] and Mindlin [7] which assumes that normals remain straight, but not necessarily orthogonal to the middle plane after deformation.
In both theories the middle plane is taken as the reference plane (z = 0) for deriving the plane kinematic equations. In this work, we consider the Reissner-Mindlin theory whose fundamental hypotheses are the following: (i) on the middle plane (z = 0) the in-plane displacements vanish, i.e. u(x, y, z = 0) = v(x, y, z = 0) = 0 that implies that points located in the middle-plane only moves vertically; (iii) the plate thickness remains unchanged; (iv) the plane stress assumption remains valid, i.e. σ zz = 0 and (v) a straight line normal to the undeformed middle plane remains straight but not necessarily orthogonal to the middle plane after deformation.
From these assumptions the displacement field can be written as: where w is the vertical displacement (deflection) of the points on the middle plane and the rotations θ x and θ y coincide with the angles followed by the normal vectors contained in the planes xz and yz respectively in their motions. We define the generalized displacement vectorû defined at any point on the middle plane.
Injecting plate theory assumptions into the 3D elastostatic problem weak form, Eq. (5) reduces to the following 2D formulation [1] whose standard finite element discretization leads to where for notational simplicity the hat symbol (•) is omitted. In the previous expression (9), K is the stiffness matrix and U and F are the vector of the generalized displacements and forces, the former containing nodal rotations and deflections and the last the dual quantities: the nodal moments and vertical nodal forces. The 3D displacement field can be then recovered by using the relations (6). In many cases, the complexity of the solution makes impossible the introduction of pertinent hypotheses for reducing the dimensionality of the model from 3D to 2D. In that case a fully 3D descriptions seems compulsory, and the in-plane-out-of-plane separated representations become particularly suitable.

In-plane-out-of-plane separated representation
The in-plane-out-of-plane separated representation was proposed in [3] and applied to efficiently solve 3D elastic problems in plate geometries. The elastic problem was defined in a plate domain = xy × z with (x, y) ∈ xy , xy ⊂ R 2 and z ∈ z , z ⊂ R. The separated representation of the displacement field u = (u, v, w) consists in a finite sum decomposition on N terms, each one of them consisting in the product of two unknown functions, one depending on the in-plane coordinates (x, y) and one on the out-of-plane coordinate z, i.e.: Expression (10) can be written in a more compact form by using the Hadamard (component-to-component) product: The separated representation of the displacement field Eq. (10) leads to a separated representation of its derivatives and consequently of the strain tensor, as proved in [3]. Then, using the constitutive equation and the separated representation of both the trial and test displacements and the corresponding strains, the separated representation constructor proceeds as described in [3]. The procedure was easily generalized to the case of elastic problems defined in shell geometries [8] and also extended to dynamics in [9].  However the proposed in-plane-out-of-plane separated representation remains too intrusive to be implemented in structural mechanics commercial softwares that instead, propose a variety of plate and shell elements.
For this reason in the next section we propose a method which allows enriching standard plate and shell discretizations while limiting intrusivity.

PGD-based enriched elements
We assume the 2D mesh defined on the middle plane of a plate geometry , depicted in Fig.  2. Our goal here is to address one element of the mesh, for example the one whose boundary is highlighted in red and that is noted by e , by using a fully 3D description in order to extract its homogenized 9 × 9 element stiffness matrix corresponding to e , and therefore fully compatible with the plate kinematics enforced on the element boundary ∂ e .
For that purpose, e is 3D resolved using the PGD-based in plane-out-of-plane separated representation while a kinematics compatible with the plate kinematics (6) is enforced on the element boundary. All the other elements in depicted in Fig. 2 are described using the standard plate theory, the only 3D resolved is just the e whose 3D mesh is depicted in Fig. 3.
The separated representation of the 3D displacement field in e is denoted by u e = (u e , v e , w e ), whose first mode consists of the standard plate kinematics where N 1 (x, y), N 2 (x, y) and N 3 (x, y) are the shape functions related to the 2D linear triangular element e , and are the nine degrees of freedom associated to its three vertices, depicted in Fig. 4. Once the first mode has been imposed the following modes are constructed as described in section while ensuring that they vanish on the element boundary ∂ e . The last condition can be enforced by using the so-called bubble function that results from the product of the three triangle shape functions B(x, y) = N 1 (x, y)N 2 (x, y)N 3 (x, y) (even if other alternatives exist for that purpose [10] Thus, the solution on the triangle boundary ∂ e is given by the first mode, that is the standard plate kinematics, whereas inside e the solution becomes enriched with models that can describe richer kinematics than the one related to standard plate kinematics. In order to extract the effective enriched 9 × 9 element stiffness K e (related to element e ) to be assembled into the global stiffness matrix involved in the algebraic system (9), we consider the element average from which the element elastic energy U e becomes Now, with the strain defined from (2) we assume the existence of a localization tensor L such that that using Eq. (2) results where U e , as previously indicated, represents the plate generalized displacement degrees of freedom defined in (13). Thus, the components of the localization tensor result from the elastic problem solution in e by prescribing the canonical boundary displacements, i.e. U e 1 = (1, 0, 0, 0, 0, 0, 0, 0, 0), U e 2 = (0, 1, 0, 0, 0, 0, 0, 0, 0), and so on, with nine associated 3D elastic problems solved by using the PGD-based in-plane-out-of-plane separated representation.
Then, by defining M = GL we obtain which allows defining the element stiffness matrix K e as In order to check the procedure, the stiffness matrix K e obtained using the first mode of the displacement expansion is compared to that one corresponding to the plate kinematics. In that case, the resulting stiffness matrix was, as expected, the one related to standard plate theory.
Of course, when considering more modes in the separated representation of the 3D displacement field u e , the resulting stiffness matrix K e differs from the one associated to standard plate theory by including eventual 3D effects ignored in standard plate kinematics. However, in that case the expression of the effective enriched stiffness matrix K e depends on the element geometry. To avoid the necessity of calculating online stiffness matrices for each triangle, the PGD rationale allows including different parameters into the separated representation as proved in our former works. The interested reader can refer to [10], [11] for an intrusive and a non-intrusive approach, respectively. In the case here addressed, the natural choice consists of considering the triangle geometry as parameters, and then as extra-coordinates within the PGD rationale.

Static condensation based enrichment
In this section, we consider an alternative procedure for computing the homogenized element stiffness matrix using static condensation [12].
We consider again the problem defined in the domain e depicted in Fig. 3. The nodal degrees of freedom used for discretizing the displacement field u e , U e , can be decomposed in the ones related to internal nodes U e i and the ones related to nodes located on ∂ e , U e b . Thus, the standard 3D finite element formulation with can be expressed as Now, by developing the first row of the previous system, the internal degrees of freedom can be expressed from the ones on the element border, that inserted in the second leads to with Now we enforce on the element border a kinematic compatible with the one of the plate theory, that is, the border nodal displacements U e b are expressed form U e [(the plate theory degrees of freedom defined at the triangle vertices as expressed in Eq. (13)], according to Eq. (6), relation that is expressed in matrix form from U e b = BU e . Thus, Eq. (25) can be rewritten as that premultiplying by B T (looking for the Galerkin based discrete system) results that allows extracting the expression of the effective plate element stiffness matrixK e and the effective plate nodal forcesF e Our numerical simulation allowed proving that the effective plate element stiffness matrixK e obtained using the just described rationale based on the static condensation coincides with the one previously obtained by using the PGD-based separated representation K e .
An additional advantage of this second route is the fact of deriving an expression for the effective plate nodal forces that will be advantageously considered when addressing inelastic and dynamic behaviors.

Numerical validation
We solve the elastic problem in the plate domain depicted in Fig. 5 and compare the solutions obtained using 3D FEM, the standard plate theory and the enriched formulations considered in the previous section. The domain geometry Table 1 whereas Table 2 specifies the considered mesh.
On the domain right face (blue area in Fig. 5) a vertical traction is applied, T = (0, 0, 10000)N/m 2 whereas the displacement is prevented on the opposite face. Figure 6 compares the different computed solutions. In that figure "3D solution" refers to the solution obtained using 3D elements according to the mesh specified in Table 2. As reference solution we consider a 3D finite element solution using a much finer mesh (the one defined in Table 3).
The number of elements with enriched behavior (according to the procedures described in the previous section) can be arbitrarily increased, and it was noted, that as expected, by increasing it the solution approaches the reference one. It is important to note that even when all the element are enriched, the solution exhibits a gap with respect to the reference solution. This gap can be explained by the fact that despite of the valuable enrichments introduced at the elements level, at the elements boundaries standard plate kinematics is enforced with its consequent impact on the resulting kinematic that continue to be too constrained with respect to a fully 3D kinematics.
Thus one could conclude that the proposed technique could become advantageous when used for enriching the mechanical description inside an element or patch (as discussed Holed 3D triangular element and its associated mesh later) but whose effects remain confined in the interior of that element (or patch) and when approaching to the element border the plate kinematic is accurate enough.

Table 3 Parameters values used in the simulation for the reference solution
In the case here discussed plate theory and coarse 3D finite element solutions are closer one to the other, whereas the enriched one is the closest to the reference one (refined 3D finite element solution).
We consider the same problem as in the previous example but now we suppose that in the 3D element there is an hole as depicted in Fig. 7 (in the other elements we consider the enrichment presented before in absence of holes).
Here, there is clearly a double advantage in using the enrichment methodology. First, if the effect of the hole remains confined inside the element, vanishing when approaching the element boundary, the element enriched suffices and the plate meshing can be alleviated because the hole only exist for the element that will be accordingly enriched, but remains invisible at the plate level. On the other hand, its local 3D effects will be described very accurately.
Once the holed element stiffness matrix K e is computed using the procedure presented in the previous section, the problem is solved and the generalized displacement field U e can be extracted at the vertices of the holed triangle. Then the 3D solution inside the element can be easily computed by using (17). Figure 8 depicts different solutions. Even if e accounts the hole presence (the PGD and static condensation based enrichment proceed on a mesh that describe the hole presence), the plate mesh becomes unaltered.
As discussed, after solving the problem at the plate level, using the enriched stiffness matrices for the different elements, 3D fields can be reconstructed. Thus, Figs. 9, 10 and 11 represent the out-of-plane stress tensor components, including the component σ zz that is neglected when using standard plate theory. Moreover, a parabolic evolution of the components σ xz and σ yz along the domain thickness can be noticed in Fig. 11.

Extension of the method to patches
In this section we extend the method previously presented to the construction of an enriched patch (super-element). In fact a single triangle element is too small for capturing rich events which can occur on a large area of the domain. Our goal is to select the  Fig. 9 Out-of-plane components of the stress tensor around the hole patch delimited by the red line in Fig. 12, p , treat it using a 3D description by employing one of the two procedures previously presented (the PGD in-plane-out-of-plane separated representation or the static condensation procedure) and then construct the (3N b )×(3N b ) patch stiffness matrix (being N b the number of nodes in the border of the enriched patch). In the case here addressed numerically, the stiffness matrix size will be 84 × 84. In the construction procedure it is assumed that on the patch border the kinematic is the standard plate kinematic applied in the surrounding elements.
In the present case the generalized patch displacement reads with N b = 28. Figures 12, 13 and 14 depict respectively the patch location, its 3D representations and its connection with the surrounding plate elements.
In the numerical validation we consider again the problem depicted in Fig. 5, where the geometrical and mechanical properties are defined in Table 4. On the right face (the blue area in Fig. 5) a vertical traction is applied, T = (0, 0, 1100000)N/m 2 whereas on the opposite face the displacement is prevented. Figure 15 depicts the different solutions, where 3D FEM coarse and refined (reference) solutions are related to the meshes specified respectively in Tables 5 and 6. The solution referred as "with all patch elements" refers to the solution computed replacing all the plate elements with the 9 patches (super-elements) depicted in Fig. 16. The 3D reference The solution involving the nine patches is the closest to the reference solution. However, because at the border of the patches the usual plate kinematics is enforced, a gap between the nine-patches solution and the refined 3D FEM solution (reference solution) persists. In any case, as the patches solution involved less constraints (plate kinematics) than the  Since it is not pretended that plates should be discretized using enriched super-elements, our main interest is adopting accurate descriptions of rich behaviors that could be assumed confined inside a region (our patch). This is the route retained for replacing 3D FEM discretization in regions exhibiting rich behaviors by enriched super-elements keeping its 2D computational complexity.

Extension of the method to plasticity
In this section we extend the method to problems in which plastic behavior can occur. In fact localized plasticity phenomena can occur in several situation, as in spot-welds during crash simulations.
In its general form the infinitesimal relation between the stress increment dσ and the elastic strain d e increment reads [13] where C is the Hooke tensor, d is total strain and d p is the plastic strain.  In plasticity yielding can occur only if the stresses satisfy a general yield criterion; in the considered examples, for the sake of simplicity, we use the Von Mises criterion [14], assuming perfect plasticity. Ignoring volumetric body forces, the weak form of the problem, in its finite incremental form, associated to the strong form (4) lies in looking for the displacement field increment u verifying the Dirichlet boundary conditions, verifying for any test function u * in an appropriate functional space. In order to solve the resulting nonlinear problem, various computational procedures have been proposed and extensively used, among them [15][16][17][18][19][20][21][22].
In this work we consider a simple implicit approach that at the n load step solves that allows updating the stress and from it computing the plastic strain increment 1 n p using well known procedures.
Then, the stress is updated according to Then, 2 n is calculated from and from it computing the plastic strain increment 2 n p using well known procedures.
Then, the stress is updated according to At iteration j the problem to be solved reads where it can be noticed that the problem structure remains unaltered, with the left-hand member encountered when addressing elastic behaviors, the inelastic behavior appearing as a body force. Thus, the enrichment procedure based on the static condensation seems specially suitable. For validating the proposed strategy we consider the same problem as in the previous section, sketched in Fig. 5, with geometrical and mechanical properties specified in Table  4 and the considered mesh defined in Table 5. On the right face of the domain (blue area in Fig. 5) a vertical traction is applied, T = (0, 0, 1100000)N/m 2 , and on the opposite face displacement is prevented. The considered uniaxial stress yield is σ 0 = 250 · 10 6 N/m 2 . Again, Fig. 17 depicts the usual different solutions. Again "all patch elements" refers to the case in which the domain is composed by the nine super-elements shown in Fig. 16. The 3D reference solution has been computed using a mesh as fine as with the one used in the patch description (Table 6). Again, the nine-patches solution is the most accurate (when compared with the reference solution).

Extension to structural dynamics
In this section we address structural dynamic, again ignoring body forces without loss of generality. Thus, the displacement field evolution u(x, t) in the domain and time interval t ∈ I = (0, T ] is described by the linear momentum balance equation where ρ is the density (kg/m 3 ). The boundary ∂ is decomposed according to ∂ = ∂ D ∪ ∂ N where displacement and tractions are prescribed. In the elastic case the problem weak form associated with the strong form (43) results in looking for the displacement field u verifying the initial and Dirichlet boundary conditions, and fulfilling for any test function u * in an appropriate functional space.    An enriched stiffness can be derived at the element or patch levels, where the use of the the procedure based on the condensation allows properly addressing the forces p * j+1 , and therefor combining dynamics and plasticity.
For evaluating the performances of the proposed enrichment procedure, we consider again the problem defined in Fig. 5 and Table 7.
On the right face (the blue area in Fig. 1) a vertical traction is enforced, T(t) = (0, 0, 4.1 · 10 6 sin(2πωt))N/m 2 , whereas on the opposite face displacement is prevented. Without loss of generality homogeneous initial conditions u(x, t = 0) = 0 and v(x, t = 0) = 0 are assumed. Table 8 reports the considered coarse mesh whereas Table 9 defined the refined one from which the reference solution is calculated. Figures 18 and 19 depict respectively the solutions computed using the different techniques at four different times and the vertical displacement time evolution at the right border. Figures 20 and 21 present similar results but in a case where plasticity takes place. Again the best solutions are the ones related to a nine-patches discretization for the reasons widely exposed before.

Conclusions
This paper proposed two different procedures to capture 3D behaviors by enriching elements (or super-elements) connected to plate discretization. This goal is performed by integrating, in a non-intrusive manner, 3D elements (or patches) in plate or shell based commercial codes. In order to compute this enriched 2D element, two different techniques were presented: one based on the PGD in-plane-out-of-plane separated representation and the other on the static condensation. The method was firstly developed in linear elastic settings and then successfully extended to structures exhibiting inelastic behaviors or dynamics.