On the coupling of local 3D solutions and global 2D shell theory in structural mechanics
Giacomo Quaranta1,
Mustapha Ziane1,
Fatima Daim1,
Emmanuelle Abisset-Chavanne2,
Jean-Louis Duval1 &
Francisco Chinesta2
Advanced Modeling and Simulation in Engineering Sciences volume 6, Article number: 1 (2019) Cite this article
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.
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, \(\Omega = \Omega _{xy} \times \Omega _z\), with \(\Omega _{xy} = [0,H_x] \times [0,H_y]\) and \(\Omega _z = [0,H_z]\) in which the thickness (out-of-plane) dimension is much lower than the other ones, i.e. \(H_z \ll H_x, H_y\).
Example of a plate domain
The linear elastic behavior relating the Cauchy's stress \(\varvec{\sigma }\) and the strain \(\varvec{\varepsilon }\) tensors reads
$$\begin{aligned} \varvec{\sigma } = \mathbb {C} \, \varvec{\varepsilon }, \end{aligned}$$
where \({\mathbb {C}}\) is the Hooke's fourth order tensor. The relation between strain \(\varvec{\varepsilon } \) and displacement \({\mathbf {u}}\) (with components \({\mathbf {u}}=(u,v,w)\)) writes
$$\begin{aligned} \varvec{\varepsilon } = \nabla _s {\mathbf {u}} = {\mathbf {G}} {\mathbf {u}}, \end{aligned}$$
where \({\mathbf {G}} = \nabla _s \bullet = \frac{1}{2} (\nabla \bullet + \nabla ^T \bullet )\) is the symmetric gradient operator. Considering an homogeneous and isotropic material and using the Voigt notation, the Hooke's tensor can be written as
$$\begin{aligned} {\mathbf {C}} = \frac{E}{(1+\nu )(1-2\nu )} \begin{bmatrix} 1-\nu&\nu&\nu&0&0&0 \\ \nu&1-\nu&\nu&0&0&0 \\ \nu&\nu&1-\nu&0&0&0 \\ 0&0&0&\frac{(1-2\nu )}{2}&0&0 \\ 0&0&0&0&\frac{(1-2\nu )}{2}&0 \\ 0&0&0&0&0&\frac{(1-2\nu )}{2} \end{bmatrix} . \end{aligned}$$
In absence of volumetric body forces, the displacement field evolution \({{\mathbf {u}}} ({\mathbf {x}})\) for \({\mathbf {x}} \in \Omega \) is described by the linear momentum balance equation
$$\begin{aligned} \nabla \cdot \varvec{\sigma } = {\mathbf {0}}. \end{aligned}$$
The domain boundary \(\partial \Omega \) is partitioned into Dirichlet, \(\Gamma _D\), and Neumann, \(\Gamma _N\), boundaries, where displacement \({\mathbf {u}}_g\) and tractions \({\mathbf {T}}\) are enforced respectively.
The problem weak form associated to the strong form (4) lies in looking for the displacement field \({\mathbf {u}}\) verifying the Dirichlet boundary conditions such that the weak form
$$\begin{aligned} \int _{\Omega } \varvec{\varepsilon }({\mathbf {u}}^*) \cdot \left( {\mathbf {C}}\, \varvec{\varepsilon }({\mathbf {u}}) \right) \, d{\mathbf {x}} = \int _{\Gamma _N} {\mathbf {u}}^* \cdot {\mathbf {T}} \, d{\mathbf {x}} \end{aligned}$$
fulfills for any test function \({\mathbf {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. \(\sigma _{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:
$$\begin{aligned} {\left\{ \begin{array}{@{}l@{}} u(x,y,z) = -z \theta _x(x,y)\\ v(x,y,z) = -z \theta _y(x,y)\\ w(x,y,z) = w(x,y) \end{array}\right. } \end{aligned}$$
where w is the vertical displacement (deflection) of the points on the middle plane and the rotations \(\theta _x\) and \(\theta _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 \(\mathbf {\hat{u}} \)
$$\begin{aligned} \mathbf {\hat{u}} = [\theta _x, \theta _y, w]^T \end{aligned}$$
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]
$$\begin{aligned} \int _{\Omega _{xy}} \varvec{\hat{\varepsilon }}(\mathbf {\hat{{\mathbf {u}}}^*}) \cdot \left( \hat{{\mathbf {C}}} \varvec{\hat{\varepsilon }}(\mathbf {\hat{u}}) \right) \, d{\mathbf {x}} = \int _{\partial _N \Omega _{xy}} \mathbf {\hat{{\mathbf {u}}}^*} \cdot \hat{{\mathbf {T}}} \, d{\mathbf {x}}, \end{aligned}$$
whose standard finite element discretization leads to
$$\begin{aligned} {\mathbf {K}} {\mathbf {U}} = {\mathbf {F}} \end{aligned}$$
where for notational simplicity the hat symbol (\(\hat{\bullet }\)) is omitted. In the previous expression (9), \({\mathbf {K}}\) is the stiffness matrix and \({\mathbf {U}}\) and \({\mathbf {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 \(\Omega = \Omega _{xy} \times \Omega _z\) with \((x,y)\in \Omega _{xy}\), \(\Omega _{xy} \subset {\mathbb {R}}^2\) and \(z \in \Omega _z\), \(\Omega _z \subset {\mathbb {R}}\). The separated representation of the displacement field \({\mathbf {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.:
$$\begin{aligned} {\mathbf {u}} (x,y,z)=\left( \begin{array}{c} u(x,y,z)\\ v(x,y,z)\\ w(x,y,z) \end{array}\right) \approx \sum _{i=1}^{N}{\left( \begin{array}{c} u^{i}_{xy}(x,y) \cdot u^{i}_z(z)\\ v^{i}_{xy}(x,y) \cdot v^{i}_z(z)\\ w^{i}_{xy}(x,y) \cdot w^{i}_z(z) \end{array}\right) }. \end{aligned}$$
Expression (10) can be written in a more compact form by using the Hadamard (component-to-component) product:
$$\begin{aligned} {\mathbf {u}}(x,y,z) \approx \sum _{i=1}^{N} {\mathbf {U}}^i_{xy} (x,y) \circ {\mathbf {U}}^i_z(z). \end{aligned}$$
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.
Enriched formulations
PGD-based enriched elements
We assume the 2D mesh defined on the middle plane of a plate geometry \(\Omega \), 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 \(\Omega _e\), by using a fully 3D description in order to extract its homogenized \(9\times 9\) element stiffness matrix corresponding to \(\Omega _e\), and therefore fully compatible with the plate kinematics enforced on the element boundary \(\partial \Omega _e\).
For that purpose, \(\Omega _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 \(\Omega \) depicted in Fig. 2 are described using the standard plate theory, the only 3D resolved is just the \(\Omega _e\) whose 3D mesh is depicted in Fig. 3.
Plate domain \(\Omega \) and its associated mesh
3D mesh of element \(\Omega _e\)
Nodal degrees of freedom defining the kinematics on \(\partial \Omega _e\)
The separated representation of the 3D displacement field in \(\Omega _e\) is denoted by \({{\mathbf {u}}^e} = ( u^e, v^e,w^e)\), whose first mode consists of the standard plate kinematics
$$\begin{aligned} {\left\{ \begin{array}{@{}l@{}} u^{e,1}_{xy}(x,y) = \theta _1^x N_1(x,y) + \theta _2^x N_2(x,y) + \theta _3^x N_3(x,y)\\ v^{e,1}_{xy}(x,y) = \theta _1^y N_1(x,y) + \theta _2^y N_2(x,y) + \theta _3^y N_3(x,y)\\ w^{e,1}_{xy}(x,y) = \omega _1 N_1(x,y) + \omega _2 N_2(x,y) + \omega _3 N_3(x,y)\\ u^{e,1}_z(z) = -z\\ v^{e,1}_z(z) = -z\\ w^{e,1}_z(z) = 1 \end{array}\right. } \end{aligned}$$
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 \(\Omega _e\), and
$$\begin{aligned} \mathbf {U^e} = (\omega _1, \theta _1^x, \theta _1^y, \omega _2, \theta _2^x, \theta _2^y, \omega _3, \theta _3^x, \theta _3^y) \end{aligned}$$
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 \(\partial \Omega _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 \({\mathcal {B}}(x,y)= N_1(x,y) N_2(x,y) N_3(x,y)\) (even if other alternatives exist for that purpose [10]), i.e. \(\forall i>1\)
$$\begin{aligned} {\left\{ \begin{array}{@{}l@{}} u^{e,i}_{xy}(x,y) = {\mathcal {B}}(x,y) P^i_x(x,y)\\ v^{e,i}_{xy}(x,y) = {\mathcal {B}}(x,y) P^i_y(x,y)\\ w^{e,i}_{xy}(x,y) = {\mathcal {B}}(x,y) P^i_z(x,y)\\ u^{e,i}_z(z) = T^i_x(z)\\ v^{e,i}_z(z) = T^i_y(z)\\ w^{e,i}_z(z) = T^i_z(z) \end{array}\right. } \end{aligned}$$
Thus, the solution on the triangle boundary \(\partial \Omega _e\) is given by the first mode, that is the standard plate kinematics, whereas inside \(\Omega _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\times 9\) element stiffness \(\mathbf {{K}}^e\) (related to element \(\Omega _e\)) to be assembled into the global stiffness matrix involved in the algebraic system (9), we consider the element average
$$\begin{aligned} \left\langle \bullet \right\rangle = \int _{\Omega _e} \bullet \ d{\mathbf {x}}, \end{aligned}$$
from which the element elastic energy \({\mathcal {U}}^e\) becomes
$$\begin{aligned} {\mathcal {U}}^e = \left\langle {\varvec{ \epsilon }}^{e^T} \varvec{ \sigma }^e \right\rangle = \left\langle \varvec{ \epsilon }^{e^T} {\mathbf {C}} \varvec{ \epsilon }^e \right\rangle . \end{aligned}$$
Now, with the strain defined from (2) we assume the existence of a localization tensor \({\mathbf {L}}\) such that
$$\begin{aligned} {\mathbf {u}}^e = {\mathbf {L}} {\mathbf {U}}^e, \end{aligned}$$
that using Eq. (2) results
$$\begin{aligned} \varvec{ \epsilon }^e = \nabla _s {\mathbf {u}}^e = {\mathbf {G}} {\mathbf {L}} \mathbf U^e, \end{aligned}$$
where \({\mathbf {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 \(\Omega _e\) by prescribing the canonical boundary displacements, i.e. \({\mathbf {U}}^e_1 = (1,0,0,0,0,0,0,0,0)\), \({\mathbf {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 \({\mathbf {M}} = {\mathbf {G}} {\mathbf {L}} \) we obtain
$$\begin{aligned} {\mathcal {U}}^e = \left\langle {\mathbf {U}}^{e^T} {\mathbf {M}}^T {\mathbf {C}} {\mathbf {M}} {\mathbf {U}}^e \right\rangle = {\mathbf {U}}^{e^T} \left\langle \mathbf M^T {\mathbf {C}} {\mathbf {M}} \right\rangle {\mathbf {U}}^e, \end{aligned}$$
which allows defining the element stiffness matrix \({\mathbf {K}}^e\) as
$$\begin{aligned} \mathbf {{K}}^e = \left\langle {\mathbf {M}}^T {\mathbf {C}} {\mathbf {M}} \right\rangle . \end{aligned}$$
In order to check the procedure, the stiffness matrix \({\mathbf {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 \({\mathbf {u}}^e\), the resulting stiffness matrix \({\mathbf {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 \({\mathbf {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 \(\Omega _e\) depicted in Fig. 3. The nodal degrees of freedom used for discretizing the displacement field \({\mathbf {u}}^e\), \({\mathbb {U}}^e\), can be decomposed in the ones related to internal nodes \(\mathbb U^e_i\) and the ones related to nodes located on \(\partial \Omega _e\), \({\mathbb {U}}^e_b\).
Thus, the standard 3D finite element formulation
$$\begin{aligned} {\mathbb {K}}^e {\mathbb {U}}^e = {\mathbb {F}}^e, \end{aligned}$$
$$\begin{aligned} {\mathbb {K}}^e&= \begin{pmatrix} {\mathbb {K}}^e_{ii} &{}{\mathbb {K}}^e_{ib}\\ {\mathbb {K}}^e_{bi} &{} {\mathbb {K}}^e_{bb} \end{pmatrix} , \end{aligned}$$
can be expressed as
$$\begin{aligned} \begin{pmatrix} {\mathbb {K}}^e_{ii} &{}{\mathbb {K}}^e_{ib}\\ {\mathbb {K}}^e_{bi} &{} {\mathbb {K}}^e_{bb} \end{pmatrix} \begin{pmatrix} {\mathbb {U}}^e_i\\ {\mathbb {U}}^e_b \end{pmatrix} = \begin{pmatrix} {\mathbb {F}}^e_i\\ {\mathbb {F}}^e_b \end{pmatrix}. \end{aligned}$$
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,
$$\begin{aligned} {\mathbb {U}}^e_i = {\mathbb {K}}^{e^{-1}}_{ii} \left( {\mathbb {F}}^e_i - {\mathbb {K}}^e_{ib} {\mathbb {U}}^e_b \right) , \end{aligned}$$
that inserted in the second leads to
$$\begin{aligned} {\mathcal {K}}^e_{bb} {\mathbb {U}}^e_b = {\mathcal {F}}_b^e, \end{aligned}$$
$$\begin{aligned} {\mathcal {K}}^e_{bb}&= {\mathbb {K}}^e_{bb} - {\mathbb {K}}^e_{bi} {\mathbb {K}}^{e^{-1}}_{ii} {\mathbb {K}}^e_{ib} \end{aligned}$$
$$\begin{aligned} {\mathcal {F}}^e_b&= {\mathbb {F}}^e_b - {\mathbb {K}}^e_{bi} \mathbb K^{e^{-1}}_{ii} {\mathbb {F}}^e_i . \end{aligned}$$
Now we enforce on the element border a kinematic compatible with the one of the plate theory, that is, the border nodal displacements \({\mathbb {U}}^e_b\) are expressed form \({\mathbf {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 \({\mathbb {U}}^e_b = {\mathbf {B}} \mathbf U^e\).
Thus, Eq. (25) can be rewritten as
$$\begin{aligned} {\mathcal {K}}^e_{bb} {\mathbf {B}} {\mathbf {U}}^e = {\mathcal {F}}_b^e, \end{aligned}$$
that premultiplying by \({\mathbf {B}}^T\) (looking for the Galerkin based discrete system) results
$$\begin{aligned} {\mathbf {B}}^T {\mathcal {K}}^e_{bb} {\mathbf {B}} {\mathbf {U}}^e = {\mathbf {B}}^T {\mathcal {F}}_b^e, \end{aligned}$$
that allows extracting the expression of the effective plate element stiffness matrix \(\tilde{{\mathbf {K}}}^e\) and the effective plate nodal forces \(\tilde{{\mathbf {F}}}^e\)
$$\begin{aligned} \tilde{{\mathbf {K}}}^e&= {\mathbf {B}}^T {\mathcal {K}}^e_{bb} {\mathbf {B}} \end{aligned}$$
$$\begin{aligned} \tilde{{\mathbf {F}}}^e&= {\mathbf {B}}^T {\mathcal {F}}^e . \end{aligned}$$
Our numerical simulation allowed proving that the effective plate element stiffness matrix \(\tilde{{\mathbf {K}}}^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 \({\mathbf {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 \(\Omega = [0,H_x] \times [0,H_y] \times [0,H_z]\) is defined in Table 1 whereas Table 2 specifies the considered mesh.
On the domain right face (blue area in Fig. 5) a vertical traction is applied, \({\mathbf {T}} = (0,0,10000) \mathrm {N/m^2}\) whereas the displacement is prevented on the opposite face.
Elastic problem defined in a plate domain
Table 1 Problem parameters
Table 2 Parameters values used in the simulation
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).
Displacement of the plate middle plane using plate theory (a), all plate theory elements except the homogenized (enriched) element (b), half plate theory elements (the triangles of the mesh in Fig. 2 with the right angle on the bottom) and half homogenized (enriched) elements (the triangles of the mesh in Fig. 2 with the right angle on the top) (c), all homogenized (enriched) elements (d), fully 3D FEM on a coarse mesh (e), fully 3D FEM on fine mesh (f)
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.
Table 3 Parameters values used in the simulation for the reference solution
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 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.
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).
Holed 3D triangular element and its associated mesh
Displacement of the plate middle plane using plate theory (a), all plate theory elements except the homogenized (enriched) holed element (b), half plates theory elements (the triangles of the mesh in Fig. 2 with the right angle on the bottom) and half homogenized elements (the triangles of the mesh in Fig. 2 with the right angle on the top), one of them holed (c), all homogenized elements, one of them holed (d), fully 3D theory on the fine mesh (e)
Out-of-plane components of the stress tensor around the hole
Out-of-plane components of the stress tensor on the plane \(z = 0.51\) mm
Out-of-plane components of the stress tensor through the domain thickness in the hole neighborhood
Plate mesh
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 \({\mathbf {K}}^e\) is computed using the procedure presented in the previous section, the problem is solved and the generalized displacement field \({\mathbf {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 \(\Omega _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 \(\sigma _{zz}\) that is neglected when using standard plate theory. Moreover, a parabolic evolution of the components \(\sigma _{xz}\) and \(\sigma _{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 patch delimited by the red line in Fig. 12, \(\Omega _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 \((3 {\mathcal {N}}_b)\times (3 {\mathcal {N}}_b)\) patch stiffness matrix (being \({\mathcal {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\times 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
$$\begin{aligned} {\mathbf {U}} = (\theta _1^x, \theta _1^y, \omega _1, \theta _2^x, \theta _2^y, \omega _2, \dots , \theta _{{\mathcal {N}}_b}^x, \theta _{{\mathcal {N}}_b}^y, \omega _{{\mathcal {N}}_b}), \end{aligned}$$
with \({\mathcal {N}}_b=28\). Figures 12, 13 and 14 depict respectively the patch location, its 3D representations and its connection with the surrounding plate elements.
3D patch mesh
Patch discretization
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) \mathrm {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 solution has been computed using a mesh as fine as the one considered to describe the enriched solution at the patch level.
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 one related to enriched elements where the kinematic constrains are applied to all the elements boundaries, its accuracy is superior.
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 \varvec{\sigma }\) and the elastic strain \(d \varvec{ \epsilon }^e\) increment reads [13]
$$\begin{aligned} d \varvec{\sigma } = {\mathbf {C}} d \varvec{ \epsilon }^e = {\mathbf {C}} (d \varvec{ \epsilon }- d \varvec{ \epsilon }^p), \end{aligned}$$
where \({\mathbf {C}}\) is the Hooke tensor, \(d \varvec{ \epsilon }\) is total strain and \(d \varvec{ \epsilon }^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 \(\Delta {\mathbf {u}}\) verifying the Dirichlet boundary conditions, verifying
$$\begin{aligned} \int _{\Omega } \varvec{\varepsilon }({\mathbf {u}}^*) \cdot \left( {\mathbf {C}} (\Delta \varvec{ \epsilon }(\Delta {\mathbf {u}}) - \Delta \varvec{ \epsilon }^p(\Delta {\mathbf {u}}))\right) \, d{\mathbf {x}} = \int _{\Gamma _N} {\mathbf {u}}^* \cdot \Delta {\mathbf {T}} \, d{\mathbf {x}} \end{aligned}$$
for any test function \({\mathbf {u}}^*\) in an appropriate functional space.
Table 4 Model parameters
Displacement using plate theory (a), all plate theory elements except the patch super-element (b), all patches as depicted in Fig. 16 (c), fully 3D refined FEM (d)
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
$$\begin{aligned} \int _{\Omega } \varvec{\varepsilon }({\mathbf {u}}^*) \cdot \left( {\mathbf {C}} \Delta _n ^1 \varvec{ \epsilon }\right) \, d{\mathbf {x}} = \int _{\Gamma _N} {\mathbf {u}}^* \cdot \Delta _n {\mathbf {T}} \, d{\mathbf {x}} \end{aligned}$$
from \(\Delta _n ^1 \varvec{ \epsilon }\), the stress increment \(\Delta _n ^1 \varvec{\sigma }\) results
$$\begin{aligned} \Delta _n ^1 \varvec{\sigma } = {\mathbf {C}} \Delta _n ^1 \varvec{ \epsilon }\end{aligned}$$
that allows updating the stress
$$\begin{aligned} \varvec{\sigma }_n ^1 = \varvec{\sigma }_{n-1} + \Delta _n ^1 \varvec{\sigma }, \end{aligned}$$
and from it computing the plastic strain increment \(\Delta _n ^1 \varvec{ \epsilon }^p\) using well known procedures. Then, the stress is updated according to
$$\begin{aligned} \varvec{\sigma }_n ^1 = \varvec{\sigma }_{n}^1 - {\mathbf {C}} \Delta _n ^1 \varvec{ \epsilon }^p. \end{aligned}$$
Table 5 Coarse 3D FEM
Table 6 Refined 3D FEM considered for defining the reference solution
Then, \(\Delta _n ^2 \varvec{ \epsilon }\) is calculated from
$$\begin{aligned} \int _{\Omega } \varvec{ \epsilon }({\mathbf {u}}^*) \cdot \left( {\mathbf {C}} \Delta _n ^2 \varvec{ \epsilon }\right) \, d{\mathbf {x}} = \int _{\Omega } \varvec{ \epsilon }({\mathbf {u}}^*) \cdot \left( {\mathbf {C}} \Delta _n ^1 \varvec{ \epsilon }^p \right) \, d{\mathbf {x}}, \end{aligned}$$
and the stress update from
$$\begin{aligned} \varvec{\sigma }_n ^2 = \varvec{\sigma }_{n}^1 + \Delta _n ^2 \varvec{\sigma }, \end{aligned}$$
At iteration j the problem to be solved reads
$$\begin{aligned} \int _{\Omega } \varvec{ \epsilon }({\mathbf {u}}^*) \cdot \left( {\mathbf {C}} \Delta _n ^j \varvec{ \epsilon }\right) \, d{\mathbf {x}} = \int _{\Omega } \varvec{ \epsilon }({\mathbf {u}}^*) \cdot \left( {\mathbf {C}} \Delta _n ^{j-1} \varvec{ \epsilon }^p \right) \, d{\mathbf {x}}, \end{aligned}$$
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) \mathrm {N/m^2}\), and on the opposite face displacement is prevented. The considered uniaxial stress yield is \(\sigma _0 = 250 \cdot 10^6 \mathrm {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).
Super-element discretization
Plate theory displacement (a), all patch elements according to Fig. 16 (b) and refined 3D FEM (reference solution) (c)
Extension to structural dynamics
In this section we address structural dynamic, again ignoring body forces without loss of generality. Thus, the displacement field evolution \({{\mathbf {u}}} ({\mathbf {x}},t)\) in the domain \(\Omega \) and time interval \(t \in I = (0,T]\) is described by the linear momentum balance equation
$$\begin{aligned} \rho \ddot{{\mathbf {u}}} ({\mathbf {x}},t) = div \, \varvec{\sigma }, \end{aligned}$$
where \(\rho \) is the density (\(\mathrm {kg/m^3}\)).
The boundary \(\partial \Omega \) is decomposed according to \(\partial \Omega = \partial _D \Omega \cup \partial _N \Omega \) 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 \({\mathbf {u}}\) verifying the initial and Dirichlet boundary conditions, and fulfilling
$$\begin{aligned} \rho \int _{\Omega } {\mathbf {u}}^* \cdot \ddot{{\mathbf {u}}} \, d{\mathbf {x}} + \int _{\Omega } \varvec{\varepsilon }({\mathbf {u}}^*) \cdot \left( {\mathbf {C}} \varvec{\varepsilon }({\mathbf {u}}) \right) \, d{\mathbf {x}} = \int _{\Gamma _N} {\mathbf {u}}^* \cdot {\mathbf {T}}(t) \, d{\mathbf {x}} \end{aligned}$$
The standard FEM space discretization reads
$$\begin{aligned} {\mathbf {M}} {\mathbf {a}}(t) + {\mathbf {K}} {\mathbf {u}}(t) = {\mathbf {F}}(t), \end{aligned}$$
where \({\mathbf {a}}(t)\) represents the acceleration. The time stepping consists in calculating the acceleration and displacement at each time step \(t_{j+1} = (j+1) \Delta t\) from the ones existing at the previous time step \(t_j = j \Delta t\).
A widely considered choice consists of the Newmark method [23] in which the velocity and the displacement fields at time \(t_{j+1}\) read
$$\begin{aligned} {\mathbf {v}}_{j+1} = {\mathbf {v}}_{j} + \Delta t \left( (1-\gamma )\mathbf a _{j} + \gamma {\mathbf {a}} _{j+1} \right) , \end{aligned}$$
$$\begin{aligned} {\mathbf {u}}_{j+1} = {\mathbf {u}}_{j} + \Delta t {\mathbf {v}}_{j} + \dfrac{\Delta t ^2}{2} \left( (1- 2 \beta ){\mathbf {a}} _{j} + 2 \beta {\mathbf {a}} _{j+1} \right) , \end{aligned}$$
so that at each time step the acceleration field can be obtained by solving the linear system
$$\begin{aligned} \mathbf {K^*} {\mathbf {a}}_{j+1} = {\mathbf {p}}^*_{j+1} \end{aligned}$$
$$\begin{aligned} \mathbf {K^*} = {\mathbf {M}} + \beta \Delta t ^2 {\mathbf {K}} \end{aligned}$$
$$\begin{aligned} {\mathbf {p}}^*_{j+1} = {\mathbf {F}}_{j+1} - {\mathbf {K}} {\mathbf {u}}_{j} - \Delta t {\mathbf {K}} {\mathbf {v}}_{j} - \dfrac{\Delta t ^2}{2} (1 - 2\beta ) {\mathbf {K}} {\mathbf {a}}_{j}. \end{aligned}$$
Table 8 Coarse FEM mesh
Table 9 Refined FEM mesh considered for computing the reference solution
According to the values of \(\beta \) and \(\gamma \) the method can be explicit (conditionally stable) or implicit (unconditionally stable). In our numerical examples we considered \(\beta = 1/4\) and \(\gamma = 1/2\) defining an implicit discretization. As soon as the acceleration is available, the velocity and displacement field can be computed using Eqs. (46) and (47) respectively.
Elastic displacement using plate theory, nine-patches and 3D FEM on the refined mesh at four different time-steps
Vertical elastic displacement evolution in time at the right border
Elasto-plastic displacement using plate theory, nine-patches and 3D FEM on the refined mesh at four different time-steps
Vertical elasto-plastic displacement evolution in time at the right border
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 \({\mathbf {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, \({\mathbf {T}}(t) = (0,0,4.1 \cdot 10^6 \sin (2 \pi \omega t)) \mathrm {N/m^2}\), whereas on the opposite face displacement is prevented. Without loss of generality homogeneous initial conditions \({\mathbf {u}}({\mathbf {x}},t=0) = 0\) and \({\mathbf {v}}({\mathbf {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.
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.
Oñate E. Structural analysis with the finite element method. Linear statics. Volume 2: Beams, plates and shells. Berlin: Springer; 2010.
Ahmad S, IB M, ZO C. Analysis of thick and thin shell structures by curved finite elements. Int J Numer Methods Eng. 1970;2(3):419–51.
Bognet B, Bordeu F, Chinesta F, Leygue A, Poitou A. Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity. Comput Methods Appl Mech Eng. 2012;201–204(Supplement C):1–12.
Chinesta F, Leygue A, Bordeu F, Aguado JV, Cueto E, Gonzalez D, Alfaro I, Ammar A, Huerta A. PGD-based computational vademecum for efficient design, optimization and control. Arch Comput Methods Eng. 2013;20:31–59.
Kirchhoff G. Uber das Gleichqewicht und die Bewegung einer elastichen Scheibe. J Reine und Angewandte Mathematik. 1850;40:51–88.
Reissner E. The effect of transverse shear deformation on the bending of elastic plates. J Appl Mech. 1945;12:69–76.
Mindlin RD. Influence of rotatory inertia and shear in flexural motions of isotropic elastic plates. J Appl Mech. 1951;18(1):31–8.
Bognet B, Leygue A, Chinesta F. Separated representations of 3D elastic solutions in shell geometries. Adv Model Simul Eng Sci. 2014;1:4.
Quaranta G, Bognet B, Ibañez R, Tramecon A, Haug E, Chinesta F. A new hybrid explicit/implicit in-plane-out-of-plane separated representation for the solution of dynamic problems defined in plate-like domains. Comput Struct. 2018;210:135–44.
Chinesta F, Keunings R, Leygue A. The proper generalized decomposition for advanced numerical simulations. A primer. Berlin: Springer; 2014.
Borzacchiello D, Aguado JV, Chinesta F. Arch Computat Methods Eng (2017). https://doi.org/10.1007/s11831-017-9241-4
Wilson LE. The static condensation algorithm. Int J Numer Methods Eng. 1974;8(1):198–203.
Owen DRJ, Hinton E. Finite elements in plasticity: theory and practice. Swansea: Pineridge Press; 1980.
Rv Mises. Mechanik der festen Krper im plastisch- deformablen Zustand. Nachrichten von der Gesellschaft der Wissenschaften zu Gottingen, Mathematisch-Physikalische Klasse. 1913;1913:582–92.
C Zienkiewicz O, Valliappan S, King I. Elasto-plastic solutions of engineering problems 'initial stress', finite element approach. J Numer Methods Eng. 1969;1:75–100.
Gallagher RH. Stress analysis of heated complex shapes. ARS J. 1962;32(5):700–7.
Argyris JH. Elasto-plastic matrix displacement analysis of three-dimensional continua. J R Aeronaut Soc. 1965;69(657):633–6.
Swedlow JL, Williams ML, Yang WH. Elasto-plastic stresses and strains in cracked plates. California Institute of Technology: Firestone Flight Sciences Laboratory; 1965.
Pope GG. A discrete element method for the analysis of plane elasto-plastic stress problems. Aeronaut Q. 1966;17(1):83–104.
Reyes SF, Deere DU. Elastic-plastic analysis of underground openings by the finite element method. Lisbon: Int Soc Rock Mech Rock Eng; 1966.
Marcal PV, King IP. Elastic-plastic analysis of two-dimensional stress systems by the finite element method. Int J Mech Sci. 1967;9(3):143–55.
Popov EP, Khojasteh-Bakht M, Yaghmai S. Bending of circular plates of hardening material. Int J Solids Struct. 1967;3(6):975–88.
Newmark NM. A method of computation for structural dynamics. J Eng Mech Div. 1959;85(3):67–94.
All the authors participated in the definition of techniques and algorithms. G. Quaranta, developed moreover the simulation software, assisted by M. Ziane. All authors read and approved the final manuscript.
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 675919.
ESI Group, Parc Icade, Immeuble le Seville, 3 bis, Saarinen, CP 50229, 94528, RUNGIS CEDEX, France
Giacomo Quaranta
, Mustapha Ziane
, Fatima Daim
& Jean-Louis Duval
PIMM, ENSAM ParisTech ESI GROUP Chair on Advanced Modeling and Simulation of Manufacturing Processes, 151 Boulevard de l'Hopital, 75013, Paris, France
Emmanuelle Abisset-Chavanne
& Francisco Chinesta
Search for Giacomo Quaranta in:
Search for Mustapha Ziane in:
Search for Fatima Daim in:
Search for Emmanuelle Abisset-Chavanne in:
Search for Jean-Louis Duval in:
Search for Francisco Chinesta in:
Correspondence to Francisco Chinesta.
Quaranta, G., Ziane, M., Daim, F. et al. On the coupling of local 3D solutions and global 2D shell theory in structural mechanics. Adv. Model. and Simul. in Eng. Sci. 6, 1 (2019) doi:10.1186/s40323-019-0125-z
Received: 08 November 2018
Plate and shells theories
In-plane-out-of-plane separated representations
Advanced Plate and Shell Models