A minimally-intrusive fully 3D separated plate formulation in computational 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 the one of 2D simulations. In the present paper we propose an efficient integration of fully 3D descriptions into existing plate software.


Introduction
We consider the linear elastostatic problem defined in the plate domain = xy × z , with xy = [0, H x ] × [0, H y ] and z = [0, H z ] in which the thickness (out-of-plane) dimension is much lower than the other ones, i.e. H z H x , H y . The linear elastic behavior relating the Cauchy's stress σ and the strain ε using Voigt notation reads where C is the elasticity matrix. 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 the behavior writes With g(x) the body forces, the displacement field 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. In what follows and without loss of generality we assume T = 0 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 ε(u * ) · (C · ε(u)) dx = u * · g dx, (5) fulfills for any test function u * , with the trial and test fields defined in appropriate functional spaces.

Plate theory
The Reissner-Mindlin theory is based on the following fundamental hypotheses [1]: (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; (ii) the plate thickness remains unchanged; (iii) 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 whose standard finite element discretization leads to K xy U = B xy (9) where for notational simplicity the hat symbol (•) is omitted. In the previous expression (9), K xy is the stiffness matrix and U and B xy 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 seem compulsory, and the in-plane-out-of-plane separated representations become particularly suitable.

PGD-based in-plane-out-of-plane decomposition
The in-plane-out-of-plane separated representation was applied in our former works to efficiently solve 3D elastic problems in plate geometries [2][3][4]. 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:

Enriched formulations
As reported in the previous section plate kinematics can be written as a single-term separated decomposition with f x (z) = − z, f y (z) = − z and f z (z) = 1.
For the sake of generality we are considering generic functions f x (z), f y (z) and f z (z) assumed known, but than can be different to the ones related to the standard Reissner-Mindlin plate theory, and its associated 3D kinematics given by Eq. (12). Consequently θ x , θ y and w do not represent rotations and deflection anymore.

The displacements gradient becomes
that allows defining the strain separated form, that taking into account its symmetry reads In the case of a general material the Hooke tensor can also be written as For an homogenous material we have simply where C z is given by Eq. (3) and C xy is a tensor whose all the entries are 1, Taking this into consideration the method that we are going to explain can be used both for homogenous and not homogenous materials. For the sake of simplicity we are going to present it in the case where in expression (16) only one term appears in the sum, but it can be easily extended to involve more terms. The virtual work principle, expressed using a matrix notation, involves the internal work In the previous expression matricesĈ Now, the virtual work integral reads where Now, if we assume an approximation based on a piecewise linear interpolation on a triangular finite element, related to an in-plane mesh of xy = ∪ E e=1 e xy , with the shape Using that approximation we can express vectors i (x, y) in each element e from the generalized nodal displacements where B i,e (x, y) contains the shape functions and theirs derivatives, according to the expressions involved in the components of i (x, y), i = 1, 2.
Now, if we consider the virtual work of the body forces g(x), it involves where without loss of generality we assume , and the single-mode decomposition of the body forces given by ). The fact of considering a single mode in the decomposition of the body force is not restrictive as discussed later.
The virtual work (27) can be expressed as where matrixĴ readŝ with I the identity matrix. Now, the integral results where V e (x, y) and G e (x, y) are approximated respectively from and with R e containing the nodal values of G(x, y) and N(x, y) = [N 1 (x, y) N 2 (x, y) N 3 (x, y)], and from which, the principle of virtual works reads Remark 1 In general the displacement decomposition within the PGD rationale involves more than a single mode, however, within the updating process, when calculating the n mode, the n − 1 already computed move to the right hand member, acting as generalized body force.
Remark 2 Thus, the in-plane functions determining the kinematics can be obtained from a standard plate theory software by using the elementary rigidity and forces given respectively by K e xy and B e xy considered in expression (26) and (38).

Remark 3
If traction T = 0 the same procedure can be applied to treat the corresponding terms.

Calculation of the out-of plane functions
The expression of solution obtained in the previous section is given by where and Now, we proceed to updated the out-of-plane functions involved in W(z) from the just calculated in-plane functions V(x, y) by considering again the principle of virtual work where now in Eq. (40) the in-plane functions are assumed known and we look for the ones involved in the out-of-plane contribution W(z). Thus, the previous integral form can be integrated on xy , and then Eq. (43) reduced to a one dimensional problem in z involving as unknown functions f x (z), f y (z) and f z (z). The same rationale that was previously addressed when performing the in-plane calculations is considered again but now with the test functions given by and Now, from the updated out-of-plane functions in W(z), the in-plane functions in V(x, y) are again updated and the procedure repeats until reaching the convergence (fixed point). The procedure for computing the out-of-plane components in this minimally-intrusive way is detailed in Appendix A.

Numerical results
The problem taken into consideration is depicted in Fig. 1 Table 1 On the right boundary face of the domain (the blue zone in Fig. 1) a vertical traction is enforced, T = (0, 0, 8 · 10 9 )N/m 2 and on the opposite face homogeneous Dirichlet boundary conditions are imposed. No volumetric body forces are considered. As in the considered domain the thickness (out-of-plane) dimension is not much lower than the other ones (in-plane dimensions), the linear variation of the displacement field along the thickness described by (2) is not more true as we can notice in Fig. 2 that compares the plate solution from the fully 3D solution assumed as reference. However using the just proposed minimally-intrusive fully 3D separated plate formulation we can notice how the solution is improved. Figure 3 shows the error of the solution respect to the 3D FEM solution, computed as as a function of the number of modes. The error of the plate theory solution results ξ (u plate ) = 0.0633.
We consider now the same problem as in the previous example but this time we suppose that there is an hole in the domain. As in the previous example, in Fig. 4 it is depicted the solution computed using the different techniques. Moreover in Figs. 5, 6 and 7 different perspectives of the out-of-plane stress tensor components are depicted. Let's note how the proposed method is able to take into consideration the σ zz components, which is ignored in plate theory, and to obtain the parabolic evolution around the thickness for the σ xz and σ yz typical of a 3D solution (Fig. 7).
In Figs. 8, 9 and 10 the same quantities are computed using a 3D finite element method, proving the good accuracy of the proposed method.
Again, in Fig. 11 shows the error of the solution respect to the 3D FEM solution as a function of the number of modes. The error of the plate theory solution being ξ (u plate ) = 0.0638.

Extension of the method to elasto-plastic dynamics
In this section we extend the method to dynamics problem in which plastic behavior can occur. With g(x, t) the body forces, 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 T(t) are prescribed.
The behavior relating the Cauchy's stress σ and the elastic strain ε e reads [5] where C is the Hooke tensor, ε is total strain and ε p is the plastic strain. The problem weak form associated with the strong form (47) 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.
We consider at time t j+1 the standard explicit time integration [6] (widely considered in commercial codes) given by or, by putting all the explicit terms at the right hand side, (51) Fig. 5 Out-of-plane stress tensor components around the hole using the minimally-intrusive fully 3D separated plate formulation Fig. 6 Out-of-plane stress tensor components in the z = 65 mm plane using the minimally-intrusive fully 3D separated plate formulation Fig. 7 Out-of-plane stress tensor components around the hole for x = 146 mm and y = 97 mm using the minimally-intrusive fully 3D separated plate formulation Fig. 8 Out-of-plane stress tensor components around the hole using 3D FEM Fig. 9 Out-of-plane stress tensor components in the z = 65 mm plane using 3D FEM Recalling (28) we can write Supposing the out-of-plane functions known, the left hand side term in (51) can be expressed as where matrixĴ j+1 readŝ with I the identity matrix. Now, the integral results with Integrating in the mesh xy = ∪ E e=1 e xy , where V j+1,e (x, y) is approximated from with N(x, y) = [N 1 (x, y) N 2 (x, y) N 3 (x, y)], and Fig. 12 The elasto-plastic dynamical problem taken into consideration The different terms at the right hand side of (51) can be treated in a similar way, as already explained for the static case, so that at each time step j the virtual work principle reads Remark 4 As for the static case, the in-plane functions determining the kinematics can be obtained from a standard plate theory software by using the elementary mass and forces given respectively by M j+1,e xy and B j,e xy .
Remark 5 Again the out-of-plane functions can be obtained in a similar manner as already explained in the static case.
For evaluating the performances of the method we consider the problem defined in Fig. 12. The geometrical and mechanical properties of the plate domain are defined in Table 2. On the right boundary face of the domain (the blue zone in Fig. 12) an horizontal traction is enforced, T = (2.7 · 10 8 , 0, 0)N /m 2 and on the opposite face homogeneous Dirichlet boundary conditions are imposed. No volumetric body forces are considered. For the sake of simplicity, we use the Von Mises criterion [7], assuming a Krupkowski isotropic hardening [8] given by the formula: whereε 0 = 0.008 is the initial equivalent plastic strain,ε p is the equivalent plastic strain, K = 0.4619 GPa is a strength coefficient and p = 0.1 is the strain hardening exponent where u 0 , v 0 and w 0 are the displacements of the points on the middle plane along x, y and z respectively 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. Again the solution computed using the proposed method is able to get a 3D behavior (as the one of the 3D FEM solution) with the striction along the thickness in the zone with a smaller section, which is typical of a 3D plastic solution.

Conclusions
Here we proposed a minimally intrusive formulation of mechanical problems (linear, elasto-plastic, static and dynamic) defined in separable domains, enabling 3D solutions expressed as a finite sum decomposition involving the product of functions defined in the plane and in the thickness. The main contribution of the present work is that the calculation of functions defined in the plane, the most expensive computationally, can be ensured by a standard plate solver, whereas the solution of those defined in the thickness, defining The different numerical examples prove the procedure efficiency that allows computing 3D solutions while keeping the computational cost the one characteristic of standard 2D plate and shell calculations.

Thus, integral (66) reads
The virtual work (27) of the body forces can be expressed as where matrixÔ readŝ with I the identity matrix. Now, the integral results Integrating in the mesh z = ∪ Q q=1 q z , where W q (z) and H q (z) are approximated respectively from W q (z) = S(z)L q , and H q (z) = S(z)M q , with M q containing the nodal values of H(z) and S(z) = [N 1 (z) N 2 (z)], and Thus, it results Q q=1 q from which, the principle of virtual works reads Thus, the out-of-plane functions determining the kinematics can be obtained from a standard 1D software by using the elementary rigidity and forces given respectively by K