1 INTRODUCTION
The important role played by transverse shear and normal stresses in identifying damage mechanisms in laminated composite plates, allied to requirements of modern industries in structural optimization, has stimulated the application of models that contain full 3-D kinematics and constitutive relations. Also, the analysis of plates with cracks, holes, fillets, notches and other geometric discontinuities demands a high quality in the prediction of stress fields (^{Bathe, 1996}; ^{Cook et al., 2002}). Stress-based finite elements, derived from mixed, hybrid and hybrid-mixed formulations, are suitable 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.
2 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 Γ_{σ}.
3 ELEMENT FORMULATION
3.1 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 u_{r} collects the rigid-body displacement components and U collects functions related to the displacement basis, on gets from (7) and (10)
3.2 Stress Function Approximation
Each column of U is stated as a Papkovitch-Neuber solution of (9) for an (hypothetic) isotropic material:
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 similar 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.
3.3 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
3.4 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
The problem can be written as the symmetric linear system
The optimal relation between the numbers n_{s} and n_{z} 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 n_{s} < n_{z} - 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.
4 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 4-node 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 E_{i} are Young's moduli, ν _{ij} are Poisson's ratios and G_{ij} are shear moduli referred to the material coordinate system. For simplicity, E = 1.
The quasi-Trefftz elements are identified by LHS(d_{s} , d_{z} ), where the symbols d_{s} and d_{z} denote the degrees of the polynomial approximations adopted for stress and displacement, respectively.
4.1 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 q_{x} = 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.
Layer 0º | Layer 90º | |||||
---|---|---|---|---|---|---|
Mesh | Nastran | LHS(10,2) | Diff. | Nastran | LHS(9,2) | Diff. |
1 | 14.07 | -0.86% | 0.5074 | 8.90% | ||
2 | 13.95 | 14.09 | -1.00% | 0.557 | 0.5080 | 8.80% |
3 | 14.35 | -2.87% | 0.5265 | 5.48% |
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.
4.2 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 left edge is fixed and the right one is submitted to the action of a uniformly distributed load q_{x} = 1 (force/area), as shown in Fig. 8.
The stress intensity factor in a state of plane stress can be calculated by
where the strain energy release rate
has the derivative ∂U/∂a numerically replaced by the ratio between the strain energy change ΔU 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.3 Nonsymmetric Laminated Plate
(Bussamra et al., 2011) presented good results for the analysis of cross-ply laminate using Quasi-Trefftz element. In order to expand the application of the element, a nonsymmetric laminated plate is analyzed.
The structure showed in Figure 10 is a 6-layer [0^{o}/90^{o}/45^{o}/-45^{o}/90^{o}/0^{o}] laminated clamped square plate, with side length L = 1 and thickness h = 0.06L, subject to a uniformly distributed transverse load q_{x} (force/area). A 5×5×6 element mesh is used in LHS tests. Tables 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.
σ _{xx} | σ _{yy} | σ _{xy} | |||||||
---|---|---|---|---|---|---|---|---|---|
Layer | Nastran | LHS(7,3) | Diff. | Nastran | LHS(7,3) | Diff. | Nastran | LHS(7,3) | Diff. |
0º | -63.308 | -62.330 | 1.55% | -3.218 | -3.302 | -2.61% | -0.248 | -0.232 | 6.47% |
90º | -1.893 | -1.893 | -0.02% | -39.248 | -38.730 | 1.32% | -0.248 | -0.241 | 2.66% |
-45º | -3.766 | -3.633 | 3.53% | -3.783 | -3.660 | 3.25% | 2.951 | 2.850 | 3.43% |
45º | 3.766 | 3.502 | 7.01% | 3.783 | 3.530 | 6.69% | 2.951 | 2.731 | 7.46% |
90º | 1.893 | 1.879 | 0.72% | 39.248 | 38.490 | 1.93% | -0.248 | -0.241 | 2.54% |
0º | 63.308 | 62.080 | 1.94% | 3.218 | 3.289 | -2.20% | -0.248 | -0.232 | 6.51% |
σ _{xx} | σ _{yy} | σ _{xy} | |||||||
---|---|---|---|---|---|---|---|---|---|
Layer | Nastran | LHS(8,4) | Diff. | Nastran | LHS(8,4) | Diff. | Nastran | LHS(8,4) | Diff. |
0º | -63.308 | -62.640 | 1.06% | -3.218 | -3.309 | -2.83% | -0.248 | -0.232 | 6.35% |
90º | -1.893 | -1.866 | 1.41% | -39.248 | -38.940 | 0.78% | -0.248 | -0.243 | 2.10% |
-45º | -3.766 | -3.647 | 3.16% | -3.783 | -3.670 | 2.99% | 2.951 | 2.862 | 3.02% |
45º | 3.766 | 3.514 | 6.69% | 3.783 | 3.540 | 6.43% | 2.951 | 2.743 | 7.05% |
90º | 1.893 | 1.852 | 2.15% | 39.248 | 38.700 | 1.40% | -0.248 | -0.243 | 1.94% |
0º | 63.308 | 62.390 | 1.45% | 3.218 | 3.298 | -2.48% | -0.248 | -0.232 | 6.35% |
4.4 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 q_{x} = 1(force/area). Three different stacking sequences for this structure were analyzed: 2-layer [0^{o}/45^{o}], 2-layer [45^{o}/-45^{o}], both with thickness [0.02/0.02], and 3-layer [45^{o}/0^{o}/-45^{o}], with thickness [0.01/0.02/0.01].
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 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.
5 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.