Simulation of Stress Concentration Problems in Laminated Plates by Quasi-Trefftz Finite Element Models

Hybrid quasi-Trefftz finite elements have been applied with success to the analysis of laminated plates. Two independent fields are approximated by linearly independent, hierarchical polynomials: the stress basis in the domain, adapted from PapkovitchNeuber solution of Navier equations, and the displacement basis, defined on element surface. The stress field that satisfies the Trefftz constraint a priori for isotropic material is adapted for orthotropic materials, which leads to the term “quasi”. In this work, the hexahedral hybrid quasi-Trefftz stress element is applied to the modeling of nonsymmetric laminates and laminated composite plates with geometric discontinuities. The hierarchical prefinement is exploited.

Latin American Journal of Solids and Structures 13 (2016) 1677-1694 ble for these applications, since it is well known to present improvements of the stress prediction when compared to conventional displacement-based FEM.In particular, the hybrid-Trefftz formulation approximates stresses within the element and displacements on its boundary.The Trefftz constraint consists of assuming that the stress approximation basis is the local solution of the Navier equation (Freitas and Bussamra, 2000).Freitas and Ji (1996) presented highly accurate finite element solution procedures for simulation of singular stress fields with two-dimensional hybrid-Trefftz element.Souza and Proença (2009) presented a hybrid-Trefftz formulation for 2-D elasticity with enrichment of the initial approximations.
Trefftz method has been applied to piezoelectricity by Qin (2003), where four modified variational functionals have been presented.Jin et al. (2005) have applied Trefftz collocation and Trefftz Galerkin methods to the analysis of plane piezoelectricity.Freitas and Bussamra (2000) have formulated hybrid-Trefftz element for 3-D elasticity.The element presented low sensitivity to geometric irregularities, incompressibility, and produced good estimates for stresses and displacements.In particular, very low sensitivity to mesh distortion was reported for the brick element.The resulting linear systems are highly sparse.Bussamra et al. (2001) expanded the use of the element for elastoplastic analysis.Motivated by these features, Bussamra et al. (2014) showed that the hexahedral 3-D hybrid-Trefftz stress elements was capable to accurately model thin or thick plates made of isotropic material, with geometric discontinuities caused by cracks, holes and notches.
The above element attributes, in addition to the symmetry, high sparsity and well conditioning of the solving system, have motivated the use of the hexahedral quasi Trefftz elements presented by Bussamra et al. (2011) to handle thin laminated plates.Solution of bending and stress concentration problems is carried out and comparison of the obtained stress concentration factors and stress intensity factors is made with available results.The hierarchical p-refinement strategy is exploited in the numerical tests.It is shown that good results can be directly computed with relatively coarse meshes.

EQUATIONS FOR 3-D ELASTICITY
Let V be the domain of a typical finite element and Γ its boundary.Let Γu and Γ σ be the portions of Γ on which displacements and tractions are specified, respectively, where The field equations below (Eqs.(2-4)) summarize the equilibrium equation, strain-displacement and constitutive relations of the linear elastic response in the domain V: where vectors σ and ε collect the independent components of the stress and strain tensors, respectively, and u is the displacement vector.As a geometrically linear model is assumed, the differential equilibrium and compatibility operators D and D * are linear and adjoint.In Eq. ( 4), the stiffness matrix k collects the elastic constants.For simplicity, body forces b and residual stresses σr are not taken into account.On the boundary, where matrix N contains the components of the unit outward normal vector to Γ vectors u Γ and t Γ are the specified surface displacements on Γu and tractions on Γ σ .

Approximations
The Hybrid-Trefftz element independently approximates stress σ in V and displacement ũ on Γ σ and also on the interelement boundary Γi, respectively: where S and Z are matrices that collect the approximation functions, and vectors X and q collect the corresponding weights.
The equilibrium equation can be expressed in terms of the displacement field by known as Navier equation.Trefftz constraint requires a choice of S such that u within the element satisfies (9).The stress is then determined from the displacement by means of If the displacement is written as where ur collects the rigid-body displacement components and U collects functions related to the displacement basis, on gets from ( 7) and ( 10)

Stress Function Approximation
Each column of U is stated as a Papkovitch-Neuber solution of (9) for an (hypothetic) isotropic material: Latin American Journal of Solids and Structures 13 (2016) 1677-1694 where ψ and ϕ are vector and scalar harmonic displacement potentials, respectively, r the position vector, the gradient operator,  the Poisson's ratio of the isotropic (hypothetic) material.A simi- lar solution seems not to be available for anisotropic material.It has been shown that the strain energy clearly converges under h-refinement for each value chosen for  in Eq. ( 13), so numerical results presented from now onwards are for  = 0 (Bussamra et al., 2011).Because of that this element has been named quasi-Trefftz.
Legendre and Chebychev modified polynomials are attributed to ψ and ϕ, as explained by Freitas and Bussamra (2000).They form a complete stress approximation basis for degrees lower than seven.Whenever a linearly dependent mode is detected, it is suppressed by numerical procedures.

Boundary Displacement Function Approximation
The displacement approximation functions are formed by hierarchical monomial bases, in a local coordinate system (ξ1, ξ2) in each face of the element.The polynomial of degree ζ is defined according to the Pascal's triangle scheme.As there are three displacement components, the number of independent displacement modes on each face is

Element Matrices
The element is based on the variational expressions Substitution of the Eqs.(7-8) into Eqs.(15-16) yields the system of equations where Latin American Journal of Solids and Structures 13 (2016) 1677-1694 NS u The problem can be written as the symmetric linear system The optimal relation between the numbers ns and nz of independent functions in S and Z should be found by numerical experimentation.The rank of the matrix will be deficient and there will be spurious modes if ns < nz -6 (Fung and Tong, 2001).
Advantage has been taken of the fact that the system ( 23) is symmetric and very sparse, typically greater than 99%, and identical elements, wherever are placed, have the same submatrix F, thus, computational efforts are reduced.

NUMERICAL APPLICATIONS
Some problems have been chosen to test the quasi-Trefftz element for stress concentration in plane stress and membrane-bending coupling in laminated plates.Results are compared with Nastran 4node laminated element models.
All the examples are three-dimensional models of orthotropic plates.The mechanical properties adopted are the same ones used by Bussamra et al. (2011): where Ei are Young's moduli, νij are Poisson's ratios and Gij are shear moduli referred to the material coordinate system.For simplicity, E = 1.The quasi-Trefftz elements are identified by LHS(ds, dz), where the symbols ds and dz denote the degrees of the polynomial approximations adopted for stress and displacement, respectively.

Plate with a Circular Hole
The first set of tests is concerned with stress concentration factor calculation in a plate with a circular hole composed of a 3-layer laminate [0 o /90 o /0 o ] with thickness [0.01/0.03/0.01].The plate has total thickness h = 0.05, height H = 10h, length L = 30h and circular hole with diameter d = 0.5H.The left edge is fixed and the right one is submitted to the action of a uniformly distributed load qx = 1 (force/area)The global coordinate system is located in the plate bottom surface and the x axis follows the direction of the length of the plate, y axis the height and z axis the thickness.
Figure 1 shows the structure and three meshes used in the analyses.Figure 2 to 4 present the convergence under p-refinement for the maximum σxx tension found at the first 0 ply, where σxx is plotted against the number of degrees of freedom NDF.Relevant results of this numerical analysis are summarized in Table 1.Figures 5 to 7 present the convergence of the maximum σ11 tension on the fiber direction for mesh 1, 2 and 3 when the laminate is replaced by [45 o ]3.Relevant results of convergence for this analysis are summarized in Table 2.   Results with the Nastran laminated element are also shown for the plate composed of both laminates.Even for simple meshes, the LHS elements have presented good convergence patterns.

Cracked Plate
In this test, the stress intensity factor K is evaluated for a cracked plate.The orthotropic plate has thickness h = 1, height H = 10h, length L = 15h and one straight crack of length a = 0.5H.The where the strain energy release rate has the derivative ∂ /∂a numerically replaced by the ratio between the strain energy change Δ and the small value Δa assigned to the crack extension.In the present finite element formulation, the strain energy is evaluated from in which F refers to the assembled equation that comes from each element contribution ( 19) and X results from the equation solution ( 23).
The 48-element mesh used in the analysis is illustrated in Figure 8.The strain energy evaluated under p-refinement is presented in Figure 9 and the stress intensity factor is summarized in Table 3.     4 and 5 show for the middle of each layer stresses at x = y = L/2 using LHS(7,3) and LHS(8,4) elements, respectively.

Membrane-Bending Coupling
The membrane-bending coupling behavior is tested in a square plate, with side length L = 1.The left edge is fixed and the right one is submitted to the action of a uniformly distributed load qx = 1(force/area).Three different stacking sequences for this structure were analyzed: In this test, the deformed shapes of the plates are analyzed.To evaluate h-refinement, the quasi-Trefftz elements are used in two different meshes, with 9 and 25 elements per ply (meshes 1 and 2, respectively), as depicted in Figure 11.The origin of the global coordinate system xyz is on the plate bottom surface.Convergence of the transverse displacement at the free corner x = 1, y = 0, z = 0.02 using LHS(10,2) element is summarized in Table 6 and compared with results from a fine mesh Nastran model.Deformed shapes using Mesh 2 are displayed in Figure 12

CONCLUSIONS
The results reported in this paper clearly illustrate the efficiency of the three-dimensional model of the hybrid quasi-Trefftz stress finite element formulation in the solution of laminated plates.Many degrees of freedom are involved, but substantial memory and computing time are saved in the assemblage and problem solution by taking numerical advantage of the high sparsity.Accurate estimates for stress concentration factors and stress intensity factors in plates with geometric discontinuities are found.The tests reported also show that the element can produce good estimates for both stresses and displacements in plane stress and plate membrane-bending coupling problems.Although good convergence rates are observed even for coarse meshes around the discontinuities, the p-refinement procedure is still restricted to a complete polynomial stress basis of sixth degree.
Research on the extension of this basis is being developed to fully exploit the hierarchical nature of the formulation.

Figure 2 :
Figure 2: Convergence of the maximum σxx under p-refinement at 0 ply for mesh 1.

Figure 3 :Figure 4 :
Figure 3: Convergence of the maximum σxx under p-refinement at 0 ply for mesh 2.

Figure 5 :Figure 6 :
Figure 5: Convergence of the maximum σ11 under p-refinement at 45 ply for mesh 1.

Figure 7 :
Figure 7: Convergence of the maximum σ11 under p-refinement at 45 ply for mesh 3.
fixed and the right one is submitted to the action of a uniformly distributed load qx = 1 (force/area), as shown in Fig.8.The stress intensity factor in a state of plane stress can be calculated by

Figure 9 :
Figure 9: Convergence of strain energy under p-refinement.
and the convergence curves are shown in Figures 13 to 15.The stresses across the thickness at the plate center are depicted in Figures 16 to 24.

Figure 12 :
Figure 12: Deformed shapes of the plates with membrane-bending coupling.

Figure 16 :Figure 17 :
Figure 16: xx across the thickness at the plate center for the laminate [0 o /45 o ].

Figure 18 :
Figure 18: xy across the thickness at the plate center for the laminate [0 o /45 o ].

Figure 19 :Figure 20 :
Figure 19: xx across the thickness at the plate center for the laminate [45 o /-45 o ].

Figure 21 :
Figure 21: xy across the thickness at the plate center for the laminate [45 o /-45 o ].

Figure 22 :
Figure 22: xx across the thickness at the plate center for the laminate [45 o /0 o /-45 o ].0 2

Figure 24 :
Figure 24: xy across the thickness at the plate center for the laminate [45 o /0 o /-45 o ].

Table 2 :
Maximum σ11 on the fiber direction for the laminate [45]3.

Table 4 :
Stresses at the middle of each layer at plate center for element LHS(7,3).

Table 5 :
Stresses at the middle of each layer at the plate center for element LHS(8,4).