SciELO - Scientific Electronic Library Online

vol.13 número9Ballistic Limit of High-Strength Steel and Al7075-T6 Multi-Layered Plates Under 7.62-mm Armour Piercing Projectile ImpactStrain Rate Dependent Behavior and Modeling for Compression Response of Hybrid Fiber Reinforced Concrete índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados




Links relacionados


Latin American Journal of Solids and Structures

versão impressa ISSN 1679-7817versão On-line ISSN 1679-7825

Lat. Am. j. solids struct. vol.13 no.9 Rio de Janeiro set. 2016 


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

Flávio Luiz de Silva Bussamraa 

Eliseu Lucena Netob 

Marcos Antonio Campos Rodriguesc 

aDepartment of Aeronautical Engineering, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil. E-mail:

bDepartment of Civil Engineering, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil. E-mail:

cDepartment of Aeronautical Engineering, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil. E-mail:


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 Papkovitch-Neuber 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 p-refinement is exploited.

Keywords: Finite element analysis; Trefftz; composite structures; geometric discontinuities


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.


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.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 ur 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








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.


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 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.

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 [0o/90o/0o] 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.

Figure 1 Plate with circular hole: 44 × 3, 56 × 3 and 62 × 3 element meshes (meshes 1, 2 and 3, respectively). 

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

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

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

Table 1 Maximum σxx for the laminate [0o/90o/0o]. 

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 [45o]3. Relevant results of convergence for this analysis are summarized in Table 2.

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

Figure 6 Convergence of the maximum σ11 under p-refinement at 45( ply for mesh 2. 

Figure 7 Convergence of the maximum σ11 under p-refinement at 45( ply for mesh 3. 

Table 2 Maximum σ1 1 on the fiber direction for the laminate [45(]3

Mesh Nastran LHS(9,2) Diff.
1 3.839 0.29%
2 3.85 3.852 -0.05%
3 3.865 -0.39%

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 qx = 1 (force/area), as shown in Fig. 8.

Figure 8 Cracked plate and 48 element mesh. 

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.

Figure 9 Convergence of strain energy under p-refinement. 

Table 3 Stress intensity fator. 

a Nastran LHS(10,2) LHS(10,3) LHS(10,4)
0.5H 4.34 4.30 4.45 4.49

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 [0o/90o/45o/-45o/90o/0o] laminated clamped square plate, with side length L = 1 and thickness h = 0.06L, subject to a uniformly distributed transverse load qx (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.

Figure 10 Nonsymmetric laminated plate. 

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

σ xx σ yy σ xy
Layer Nastran LHS(7,3) Diff. Nastran LHS(7,3) Diff. Nastran LHS(7,3) Diff.
-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%
63.308 62.080 1.94% 3.218 3.289 -2.20% -0.248 -0.232 6.51%

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

σ xx σ yy σ xy
Layer Nastran LHS(8,4) Diff. Nastran LHS(8,4) Diff. Nastran LHS(8,4) Diff.
-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%
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 qx = 1(force/area). Three different stacking sequences for this structure were analyzed: 2-layer [0o/45o], 2-layer [45o/-45o], both with thickness [0.02/0.02], and 3-layer [45o/0o/-45o], 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.

Figure 11 Laminate and meshes. 

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.

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

Figure 13 Convergence of transverse displacement for the laminate [0o/45o]. 

Figure 14 Convergence of transverse displacement for the laminate [45o/-45o]. 

Figure 15 Convergence of transverse displacement for the laminate [45o/ 0 o/-45o]. 

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

Figure 17 σyy across the thickness at the plate center for the laminate [0o/45o]. 

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

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

Figure 20 σyy across the thickness at the plate center for the laminate [45o/-45o]. 

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

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

Figure 23 σyy across the thickness at the plate center for the laminate [45o/0o/-45o]. 

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

Table 6 Transverse displacement at x = 1, y = 0, z = 0.02. 

Laminate Mesh Nastran LHS (10,2) Diff
[0º/45º] 1 6.990 6.651 4.85%
2 7.010 -0.29%
[45º/-45º] 1 3.213 2.936 8.62%
2 3.235 -0.68%
[45º/0º/-45º] 1 0.390 0.367 5.90%
2 0.401 -2.82%


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.


This work has been partially supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).



Crack length


Body force

D, D*

Differential operator, and its adjoint


Young’s moduli


Shear moduli


Material stiffness matrix


Stress intensity factor

ns, nz

Number of stress and displacement modes


Unit outward normal vector matrix


Displacement weights


Uniform load (force/area) in x direction


Position vector


Stress basis for σ

t Γ

Specified surface traction


Displacement vector in V


Displacement vector on Γi and Γσ

u Γ

Specified surface displacements


Displacement basis for u


Strain energy


Domain of an element


Stress weights


Displacement basis for ũ


Boundary of an element

Γi, Γu, Γσ

Interelement, kinematic and static boundaries


Vector with strain components


Degree of the polynomial displacement ũ


Poisson’s ratios


Local coordinate system

σ, σr

Vectors with stress and residual stress components


Scalar harmonic displacement potentials


Vector harmonic displacement potentials

Gradient operator


Bathe, K.J. (1996). Finite Element Procedures. Prentice Hall, Englewood Cliffs. [ Links ]

Bussamra, F.L.S., Lucena Neto, E., Ponciano, W.M. (2014). Simulation of stress concentration problems by hexahedral hybrid-Trefftz finite element models. Computer Modeling in Engineering & Sciences, 99(3): 255-272. [ Links ]

Bussamra, F.L.S., Lucena Neto, E., Raimundo Jr., D.S. (2012). Hybrid quasi-Trefftz 3-D finite elements for laminated composite plates, Computers and Structures. 92-93:185-192. [ Links ]

Bussamra, F.L.S., Pimenta, P.M., Freitas, J.A.T. (2001). Hybrid-Trefftz stress elements for three-dimensional elastoplasticity. Computer Assisted Mechanics and Engineering Sciences, 8: 235-246. [ Links ]

Cook, R.D., Malkus, D.S., Plesha, M.E., Witt R.J. (2002). Concepts and Applications of Finite Element Analysis. John Wiley, New York. [ Links ]

Freitas, J.A.T., Bussamra, F.L.S. (2000). Three-dimensional hybrid-Trefftz stress elements. International Journal for Numerical Methods in Engineering, 47(5): 927-950. [ Links ]

Freitas, J.A.T., Ji, Z.Y. (1996). Hybrid-Trefftz finite element formulation for simulation of singular stress fields. International Journal for Numerical Methods in Engineering, 39(2), 281-308. [ Links ]

Fung, Y.C., Tong, P. (2001). Classical and Computational Solid Mechanics, World Scientific, Singapore. [ Links ]

Jin, W.G., Sheng, N., Sze, K.Y., Li, J. (2005). Trefftz indirect methods for plane piezoelectricity. International Journal for Numerical Methods in Engineering, 63(1): 139-158. [ Links ]

Qin, Q.H. (2003). Variational formulations for TFEM of piezoelectricity. International Journal of Solids and Structures, 40(23): 6335-6346. [ Links ]

Souza, C.O., Proença, S.P.B. (2009). Formulação híbrida-Trefftz com enriquecimento seletivo: aplicação a problemas bidimensionais da elasticidade. Caderno de Engenharia de Estruturas, São Carlos, 11: 125-157. [ Links ]

Funding Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).

Received: December 10, 2015; Revised: April 04, 2016; Accepted: April 04, 2016

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License