Three-dimensional stress analysis of variable angle tow composite laminate using hybrid brick elements

Thin–Walled Structures 148 (2020) 106587

Contents lists available at ScienceDirect

Thin-Walled Structures journal homepage: http://www.elsevier.com/locate/tws

Full length article

Three-dimensional stress analysis of variable angle tow composite laminate using hybrid brick elements Lala Bahadur Andraju, Gangadharan Raju * Department of Mechanical and Aerospace Engineering, Indian Institute of Technology Hyderabad, India

A R T I C L E I N F O

A B S T R A C T

Keywords: Composite laminate Stress analysis Continuous tow shearing technique Variable angle tow laminates Hybrid brick elements

Variable angle tow (VAT) composites manufactured using continuous tow shearing (CTS) technique exhibits smooth thickness variation across the plane of the laminate as a function of fiber angle variation. These VAT laminates have a flat surface on one side and a curved surface on the other side resulting in asymmetric crosssection about the mid-surface of the laminate. Mostly shell-based finite elements have been used for modeling the VAT panels, and they neglect the thickness variation introduced by manufacturing in the stress analysis of laminates. In this work, the use of hybrid brick elements for three-dimensional stress analysis of VAT laminates manufactured using CTS technique is studied. A user-defined element subroutine is written for the 8-node and 27-node hybrid brick elements in the ABAQUS finite element software. The numerical performance of the hybrid brick elements in computing the stress solution of composite laminates is evaluated by studying benchmark problems and the results are compared with the displacement brick element. The influence of mesh distortion, shear locking and boundary layer effect in the stress solution evaluated by hybrid elements are studied. VAT laminate with thickness variation is modeled using the hybrid brick elements and the effect of asymmetric crosssection on the through-thickness stress field is investigated. The influence of the fiber angle distribution as well as thickness variation on the phenomenon of stress redistribution across the plane of the laminate for different VAT configuration is studied. Further, the numerical results of VAT laminate with thickness variation obtained using hybrid brick elements are compared with the experimental results available in the literature. The results from this numerical study will help the engineer in designing efficient VAT laminates fabricated using CTS technique.

1. Introduction New advancements in the automation of composite manufacturing have enabled the aerospace, automobile and wind energy sectors to gradually move away from conventional composite fabrication methods for improved quality and high productivity. Automated fiber placement (AFP) and Continuous tow shearing (CTS) techniques have made possible the steering of fibers in the plane of the lamina which provides the engineer with more tailorability options in designing composite structures. As a result of variation of the fiber angle, there is a contin­ uous change in the stiffness properties, and such laminates are called variable angle tow (VAT) composites. Hyer et al. [1,2] applied the concept of tow steering for minimizing the stress concentration effects due to the presence of a hole in a composite laminate. Gurdal et al. [3] proposed a linear fiber angle definition for variable stiffness composites and demonstrated the improvement in the buckling performance of VAT laminates over straight fiber composites. The improvement is attributed

to the load redistribution from the center to the supported edges of the VAT plate by steering the fibers. Groh and Weaver [4,5] derived a mixed stress/displacement plate model based on Hellinger-Reissner formula­ tion to evaluate three-dimensional (3D) stress solution of VAT lami­ nates. They employed a differential quadrature method to solve the plate model, and the computed stress results were compared with the finite element (FE) analysis. Wu et al. [6,7] employed the Rayleigh-Ritz method to solve the buckling and post-buckling problem of VAT plates under compression and demonstrated the improvement in post-buckling performance over optimal straight fiber composites. Oliveri et al. [8–10] evaluated the buckling and post-buckling of VAT composite laminates using the Rayleigh-Ritz method. Also, they formulated an eXtended Ritz method (X-Ritz) to evaluate the buckling and post-buckling analysis of stiffened panels with embedded cracks. Raju et al. [11] studied the post-buckling performance of VAT plates under combined axial compression and shear loads and showed the improvement over straight fiber composites. Gunay and Timarci [12] carried out the stress analysis

* Corresponding author. Department of Mechanical and Aerospace Engineering, IIT Hyderabad, Kandi, Sangareddy, Telangana, 502285, India. E-mail address: [email protected] (G. Raju). https://doi.org/10.1016/j.tws.2019.106587 Received 28 June 2019; Received in revised form 19 December 2019; Accepted 22 December 2019 Available online 20 January 2020 0263-8231/© 2019 Elsevier Ltd. All rights reserved.

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

impact (CAI). Few works have been reported on the accurate 3D stress analysis of VAT laminates manufactured using CTS. Accurate 3D stress evaluation is necessary for the various damage modes at the micro and macro scale of the composite structure and also in the damage evolution studies under various loading conditions. Most of the work reported on the stress analysis of laminated composites are based on displacement brick elements. As there is high stress gradient in the thickness direction in laminated composites due to the lower transverse shear stiffness properties, more number of displacement brick elements are needed to obtain accurate stress results. The displacement brick elements with linear interpolation are susceptible to shear locking [17] and Poisson’s ratio locking [18]. To overcome the limitations of displacement brick elements, the hybrid brick elements are formulated based on the mixed formulation. Pian et al. [19–23] established the framework of hybrid elements based on two-field variational formulation involving displacement and stress interpolation. Many 2D and 3D hybrid elements for composite structures were developed, and their numerical performance was assessed by studying the benchmark problems. Mau et al. [24] derived a quadrilateral multilayer plate element based on the mixed formulation and included the effects of transverse shear deformation. Spilker et al. [25] formulated a couple of quadrilateral multilayered plate bend­ ing/stretching elements for the analysis of multilayered composite plates. These plate elements differ from each other based on the number of terms used in the stress and displacement interpolation functions. The hybrid plate elements are not able to compute the through-thickness stresses accurately in composite laminates. For accurate stress as well as failure analysis of composite structures, 3D brick elements are nor­ mally used in the finite element models. Jog [26] developed a 27-node hybrid brick element with 75 terms in the stress interpolation. The stress terms in the interpolation are chosen so that the element stiffness matrix has no spurious energy modes. Later, Jog [27] improved the perfor­ mance of the 27-node hybrid brick element by considering 90 stress terms for studying nonlinear transient problems. The main advantage of the hybrid element is attributed to the eval­ uation of accurate 3D stress at any point within the element rather than Gauss points (material points or integration points) in displacement brick element. The hybrid brick elements are less susceptible to the shear locking phenomenon problem present in thin laminates. Also, the hybrid element results are not influenced by the mesh distortion present in composite structures with complex shapes and discontinuities. Howev­ er, the challenge in usage of hybrid elements is attributed to the proper selection of terms in the stress interpolation to avoid numerical in­ stabilities [27]. In this paper, the 8-node and 27-node hybrid brick elements are used to model the VAT composite laminate manufactured using CTS tech­ nique. Initially, the benchmark composite laminate problems are stud­ ied using hybrid brick elements, and the computed results are compared with the analytical solution. Furthermore, the problems of shear locking, mesh distortion are investigated using hybrid brick elements. The effect of the boundary layer on the stress results near the clamped edges of unsymmetric laminates are studied using the hybrid brick elements. Sub­ sequently, the hybrid brick elements are used to study the influence of thickness variation on the 3D stress solution of the VAT laminate. Finally, the numerical performance of hybrid brick elements in analyzing VAT laminates manufactured using the CTS technique is compared with the experimental results available in the literature. The paper is organized as follows. Section 2 briefly discuss the modeling of VAT laminate with thickness variation. In Section 3, the finite element formulation of hybrid brick elements is presented. Section 4 presents the stress analysis results of VAT laminates using hybrid el­ ements and comparison with displacement brick elements. Also, the numerical performance of hybrid brick elements is compared with the experimental results of VAT laminates. Finally, conclusions are drawn in Section 5.

Fig. 1. CTS manufactured VAT lamina 0〈0j 70〉 with smooth thick­ ness variation.

on the thin-walled cantilever box beams of variable stiffness laminate subjected to transverse and torsional loading at the free end using the displacement-based FE method. They have investigated the effect of different fiber paths and evaluated the axial and transverse shear stress distributions on the beam at various locations. Most of the works on VAT laminates are based on higher-order plate/shell theories which does not take thickness variation introduced by manufacturing in the analysis. In this work, the use of 3D brick elements based on hybrid formulation is investigated for stress analysis of VAT laminates manufactured using the CTS technique. In the AFP technique, the steering of fiber is done by the in-plane bending action which causes local fiber buckling and introduces wrinkle defects in the structural component [13,14]. Besides, the AFP technique cannot steer the fibers around low radius of curvature and introduces unwanted overlaps/tow gaps in the VAT composite structure. These overlaps/gaps may affect the aerodynamic performance and also introduce unwanted stress raisers in composite structures. To overcome the limitations of AFP, Kim et al. [14] developed continuous tow shearing (CTS) technique to manufacture VAT laminates by applying in-plane shear deformation to the fibers and steer them in the plane of the laminate. In the CTS approach, the fibers can be steered around tighter radius compared to the AFP, and the resulting VAT laminates have a smooth continuous thickness variation. Groh and Weaver [15] employed a shell model to study the buckling behavior of VAT laminates manufactured using CTS technique. They compared the buckling results of VAT laminates computed using the shell model with 3D FE results. Dang et al. [16] developed a 3D finite element model to capture the features of manufactured VAT laminates using CTS. They studied experimentally and numerically the failure behavior of VAT laminates under sequential tests of the drop-weight impact and compression after 2

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

1.1. Modeling of VAT laminates manufactured by CTS

displacement. This is due to the stress being obtained directly from the minimization of functional rather than by differentiation of the displacement. In addition, the transverse stress and displacement con­ tinuity along the interfaces laminated composites can be easily enforced in the hybrid element formulation and allows accurate evaluation of the transverse stresses. Another advantage being the stresses can be calcu­ lated accurately at any point within the hybrid element, whereas the stresses are accurate only at the Gauss points in the displacement based elements. In the hybrid formulation, the displacements and the stresses are interpolated independently. The displacement field of the hybrid brick element is expressed as, 0 1 u1 1B v1 C 0 1 0 C u N1 0 0 N2 0 0 …:: B B w1 C C uj (3) u ðorÞ ui ¼Nij b u¼ @ v A ¼ @ 0 N1 0 0 N2 0 …:: AB B u2 C ¼Nb C 0 0 N1 0 0 N2 …:: B w @ : A :

Continuous tow shearing (CTS) is a recently developed technique to fabricate VAT laminates with minimal process-induced defects. The key feature of CTS is the application of in-plane shear deformation to the tow which helps in avoiding tow overlaps and gaps. This allows the fabri­ cation of VAT laminates with smooth thickness variation as shown in Fig. 1. To maintain the volume fraction of fiber, the thickness of tow inherently increases as it is sheared during the steering process. The relation between unsheared tow thickness t0 and sheared tow thickness tθ is given by Ref. [14], tθ ¼

t0 ; cosθ

(1)

where θ is the shearing angle of the tow. VAT composites manufactured by CTS have smooth thickness variation which mimics a flat surface on one side and a cylindrical surface of varying radius of curvature on the other side (see Fig. 1). In this work, VAT laminate with linear fiber angle variation is studied. The Lagrange polynomial is used to interpolate the fiber angle at pre-selected control points in the lamina. The linear fiber angle variation for each ply k is expressed as [3], � 2 T k1 T k0 θðkÞ ðxÞ ¼ φ þ (2) jxj þ T k0 ; a

0

∂u ∂x

1

0

∂N1 B B ∂x

B C B C B C B B C B B C B ∂ v B C B C B 0 0 1 B B ∂y C B εxx B C B B C B C B B εyy C B C B B C B ∂w C B B C B C B B C B ∂z C B 0 B εzz C B C B C B B C B εu ¼ B C ¼ B C¼B Bγ C B C B B xy C B ∂v ∂u C B ∂N1 B C B þ C B B C B ∂x ∂y C B ∂y B γ yz C B C B @ A B C B B C B B ∂w ∂v C B γ xz B C B B ∂y þ ∂z C B 0 B C B B C B B C B B C B @ ∂u ∂w A @ ∂N1 þ ∂z ∂x ∂z

εuij ¼

0

∂N2 ∂x

0

0

∂N1 ∂y

0

0

∂N2 ∂y

0

0

∂N1 ∂z

0

0

∂N2 ∂z

∂N1 ∂x

0

0

∂N2 ∂N2 ∂y ∂x 0

∂N1 ∂N2 ∂x ∂z

0 1

1

0

∂N1 ∂N1 ∂z ∂y

where, u1 ; v1 ; w1 ; … are the nodal degrees of freedom, Ni are the displacement interpolation functions and N is the displacement shape function matrix. The elastic linear strains are evaluated using the rela­ tion given by

0

∂N2 ∂N2 ∂z ∂y 0

∂N2 ∂x

…:: C C C C C C …:: C0 1 C u1 C CB C CB v1 C CB C B C …:: C CB C CB w1 C CB C u CB C ¼ Bb CB u C CB 2 C …:: CB C CB C CB : C C@ A C C : …:: C C C C C C A …::

B C B C B C B C B C B C B C B C B C B C B C B C B C B C B C B C BorC B C B C B C B C B C B C B C B C B C B C B C B C B C B C @ A

(4)

� 1 ui;j þ uj;i 2

where, B represents the strain-displacement matrix and the superscript u represents the strains obtained from displacements. The 3D stress field of the hybrid brick element is expressed as, 10 1 0 1 0 1 0 0 0 0 0 x 0 0 0 0 0 y …:: β1 σ xx B σ yy C B 0 1 0 0 0 0 0 x 0 0 0 0 0 …:: CB β2 C CB C B C B B σ zz C B 0 0 1 0 0 0 0 0 x 0 0 0 0 …:: CB β3 C b C B CB C σ¼B (5) B τxy C ¼ B 0 0 0 1 0 0 0 0 0 x 0 0 0 …:: CB β4 C ¼P β CB C B C B @ τyz A @ 0 0 0 0 1 0 0 0 0 0 x 0 0 …:: A@ : A 0 0 0 0 0 1 0 0 0 0 0 x 0 …:: : τxz

where, θðkÞ ðxÞ is the local fiber angle at coordinate x, Tk0 is the fiber angle Tk1

at center of the lamina x ¼ 0, is the fiber angle at the ends x ¼ � a= 2, a is the length of the lamina and φ is the angle of rotation of fiber path with respect to x-axis. The variation of fiber angle of a VAT lamina is represented by the notation θ ¼ φ �〈T0 jT1 〉. In the present work, it is assumed that the fiber angle is varying along the x-direction and remains constant in the y-direction. 2. Hybrid finite element formulation

where, P is the stress interpolation matrix and b β is a vector of stress parameters used in the interpolation. Using the constitutive relation, the strains are evaluated as,

Hybrid elements are developed to overcome the limitations of displacement based formulation elements for the accurate evaluation of stresses in laminated composites. In the hybrid/mixed formulation, both displacement and stress interpolation functions are chosen indepen­ dently. The accuracy of the stress is of the same order as the

ε ¼ Sσ ¼ SPbβ ðorÞ εij ¼ Sijkl σkl

3

(6)

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 2. Brick elements. a) 8-node brick element b) 27-node brick element.

here, S is the material compliance matrix. The variation of displacement and stress interpolation are given as δu ¼ Nδb u δσ ¼ Pδb β

Z ZΩ

(7)

Γt

ZΩ bf ¼

Ω

Γt

Kb u ¼ bf

σ ¼ Pbβ ¼ PH 1 Gb u

Ω

R

By applying the Green-Gauss theorem δui σij nj dΓ ¼ δui σij dΩ þ Γ Ω � σij δui;j dΩ , and using the symmetry of stress tensor (σ ij ¼ σ ji ), Eqn. (9)

can be simplified as Z Z Z Z h σ ij δui;j dΩ þ δui bi dΩ þ δui ti dΓ þ δσ ij εuij Ω

Ω

Γ ti

i Sijkl σkl dΩ ¼ 0

(10)

Ω

As the virtual variations δu and δσ are arbitrary and by considering two cases ðδu; δσ Þ ¼ ð0; δσ Þ and ðδu; δσ Þ ¼ ðδu; 0Þ, Eqn. (10) can be decomposed as Z δσ t ½εu Sσ �dΩ ¼ 0 (11)

the minimum number of stress interpolation terms works for static problems but results in numerical instabilities for transient problems [27,30]. Similar numerical issues related to the order of expansion of displacement and stress functions have been reported in mixed plate theory based on generalized unified formulation for static loading [31]. As expected, the solution converges in the initial load steps (or with small increments), but in the later load steps, the solution starts to deviate and results in instabilities. One reason being the stress interpo­ lation should have all the terms starting from the lowest to the highest order, which is the requirement for numerical convergence of FE solu­ tion. Another problem is the presence of common stress interpolation terms due to which the stress components are coupled and results in spurious stresses [27]. To overcome these issues, and to get the optimal stress interpolation, Jog [27] has formulated guidelines in choosing terms in the stress interpolation for the robustness of the hybrid brick elements by removing the terms associated with spurious energy modes.

Ω

Z

Z ðδεu Þt σ dΩ ¼

Ω

Z δut bdΩ þ

Ω

δut tdΓ

(16)

The challenging task in hybrid element formulation is choosing the number of terms in the stress interpolation to get accurate results. The minimum number of terms in the interpolation function must be � n r, where n is the number of degrees of freedom per element and r is the number of rigid-body motion. For example, consider an 8-node brick element which has 24 degrees of freedom and six rigid body motions. The minimum number of terms in the stress interpolation should be 24 6 ¼ 18. There is no constraint on the maximum number of terms in the stress interpolation, but adding more terms increases the stiffness of the hybrid brick element and also increases the computational cost [18,27]. The increase in the magnitude of stiffness matrix (K) results in less compliance of the hybrid element ( b u ¼ K 1 bf ). The strategy of choosing

(9) R

(15)

where, K ¼ GT H 1 G is the stiffness matrix of the hybrid element. The stress field of the hybrid brick element is evaluated using the relation

Ω

R

Nt bdΩ Ω

By eliminating the b β parameters from Eqn. (13), and simplifying gives

Ω

Γ ti

(14) Z

Nt tdΓ þ Γ

where, ti is the traction acting on the surface Γ, ti is the prescribed traction and bi is the body force acting on the domain Ω. Eqn. (8) can be written as Z Z Z Z Z h i δui σ ij;j dΩ þ δui bi dΩ þ δui ti dΓ δui ti dΓ þ δσij εuij Sijkl σ kl dΩ ¼0 Ω

Pt SPdΩ

H¼

The mixed variational form of the hybrid finite element is based on Hellinger-Reissner principle [26,28,29] and expressed as Z Z Z h i � (8) δui σij;j þ bi dΩ þ δui ðti ti ÞdΓ þ δσij εuij Sijkl σ kl dΩ ¼ 0 Ω

Pt BdΩ

G¼

(12)

Γt

where ðσ Þ6�1 ;ðεÞ6�1 , ðSÞ6�6 are the stress, strain and compliance tensors expressed in Voigt notation. After substituting Eqs. (3)–(5) and (7) in Eq. (11) and (12), we get the following matrix equations � �� � � � b H G β ¼ 0 (13) t b G 0 f b u where

4

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

smooth curved surface on the other side. For creating VAT laminate mesh with varying thickness, the bottom surface nodes along the length direction are created first and then the top surface nodes are created by using Eqn. (1) at each node point which gives the thickness information. Then, the nodes are created along the thickness direction as per the mesh information, and the elements are created which have smooth thickness variation across the lamina. Later, the fiber angle is evaluated at the center of each element using Eqn. (2), and assumed to be constant for all the Gauss points used in evaluating the element stiffness matrix. A Matlab code is written to generate the ABAQUS input file (.inp file) for hybrid brick elements to define the fiber angle and thickness values for each element of the VAT lamina. In this section, the stress fields for the composite beam problems are normalized using the following expression

Fig. 3. A simply supported composite beam subjected to sinusoidal distrib­ uted load. Table 1 Mesh density information of the brick elements. Element type

σ xx ¼

Number of elements along x, y and z-direction (x � y � z) Mesh 1

Mesh 2

Mesh 3

Mesh 4

C3D8R

20 � 1 � 6

20 � 1 � 12

20 � 1 � 18

20 � 1 � 24

C3D8RH

20 � 1 � 6

20 � 1 � 12

20 � 1 � 18

20 � 1 � 24

8H18

20 � 1 � 6

20 � 1 � 12

20 � 1 � 18

20 � 1 � 24

27H90

10 � 1 � 3

10 � 1 � 6

10 � 1 � 9

10 � 1 � 12

�a � t2 1 1 �a � σ xx ; z ; τxz ¼ τxz ð0; zÞ; σ zz ¼ σ zz ; z ; 2 q0 q0 2 q0 a2

and the stress solution for plate problems are normalized as � � � � t2 a b t2 a b σxx ¼ 2 σxx ; ; z ; σ yy ¼ 2 σ yy ; ; z ; 2 2 2 2 q0 a q0 b � � � � 1 b 1 a b τxz ¼ τxz 0; ; z ; σzz ¼ σzz ; ; z ; q0 2 q0 2 2

(17)

(18)

where q0 is the total pressure applied on the surface of the laminate. In this work, both 8-node and 27-node hybrid elements are used to perform stress analysis of VAT composite laminates manufactured using CTS technique. The stress interpolation for the 8-node and 27-node hybrid brick elements (see Fig. 2) is given in the appendix.

3.1. Composite beam subjected to sinusoidal loading Firstly, the benchmark problem of cylindrical composite beam bending proposed by Pagano [33] is studied. The material properties of the lamina are given by E1 ¼ 172:36 GPa, E2 ¼ E3 ¼ 6:89 GPa, ν12 ¼ ν23 ¼ ν13 ¼ 0:25, G12 ¼ G13 ¼ 3:44 GPa, G23 ¼ 1:37 GPa. Laminate with a stacking sequence [0/90/0] and thickness (t ¼ 125 mm) to length (a ¼ 1000 mm) ratio of 1:8 is considered. The plane strain condition was enforced by considering one element in the width direction and con­ straining the lateral edges from expansion. A sinusoidal pressure load of magnitude q0 ¼ 100 kPa is applied equally on the top and bottom surfaces of the composite beam as shown in Fig. 3. The composite beam is discretized using the C3D8R, C3D8RH, 8H18 and 27H90 brick elements and the details of the mesh sizes are given in Table 1. C3D8RH is the hybrid brick element available in the ABAQUS element library with reduced integration option. C3D8RH is based on the mixed formulation (assumed displacements and constant pressure) by writing the equilibrium equations in terms of deviatoric stresses and pressure to model incompressible materials [18,32]. The mesh sizes are chosen in such a way that the number of nodes along the x and z-di­ rections are the same for all the element models. However, the 27H90 element is quadratic and contains an additional mid-node in the width direction compared to the C3D8R, C3D8RH and 8H18 element, but the results are not affected because of the plane strain condition enforced in

3. Results and discussion In this section, the stress analysis of straight fiber and VAT laminates is performed using the hybrid brick elements and their results are compared with the displacement brick elements. A user-defined element (UEL) subroutine for the 8-node and 27-node hybrid brick elements is developed in the ABAQUS [32] finite element software. In this paper, the 8-node and 27-node hybrid brick elements are abbreviated as 8H18 and 27H90, respectively. The first number indicates the number of nodes, the letter ‘H’ represents hybrid element formulation and the last number designates the number of polynomial terms used in the stress interpolation. For the 8-node and 27-node hybrid brick elements, full integration of 2 � 2 � 2 and 3 � 3 � 3 Gauss quadrature, respectively are used to evaluate all the integrals over the element domain. The nu­ merical performance of the hybrid brick elements is compared with the 8-node displacement brick element C3D8R available in ABAQUS finite element software. As already discussed in Section 2, the VAT laminates manufactured using the CTS technique are having a flat surface on one side and a

Fig. 4. Normalized through-thickness stress distribution evaluated using C3D8R element and comparison with Pagano’s analytical solution. a) σxx b) τxz c) σ zz (z= t is normalized through-thickness location). 5

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 5. Normalized through-thickness stress distribution evaluated using C3D8RH element and comparison with Pagano’s analytical solution. a) σxx b) τxz c) σzz .

Fig. 6. Normalized through-thickness stress distribution evaluated using 8H18 element and comparison with Pagano’s analytical solution. a) σxx b) τxz c) σzz .

Fig. 7. Normalized through-thickness stress distribution evaluated using 27H90 element and comparison with Pagano’s analytical solution. a) σ xx b) τxz c) σ zz .

the y-direction. Figs. 4–7 show the through-thickness stress distribution results of the C3D8R, C3D8RH, 8H18 and 27H90 brick elements for different mesh density. The computed stress results are converging toward the Pagano’s [33] analytical solution with the increase in number of elements in the thickness direction. In the case of C3D8R brick element, the stress solutions are approaching towards the analytical solution with the increase in mesh density as shown in Fig. 4. However, the axial stress σxx does not match with the analytical solution at the top and bottom surface of the beam. Also, the transverse shear stress τxz is not converging towards zero on the top and bottom surfaces of the composite beam (see Fig. 4(b). Fig. 5 shows that results of C3D8RH element are similar to the C3D8R element. Also, the results of C3D8H, ABAQUS hybrid brick element with full integration (not presented here) is poor when compared to the reduced integration. Fig. 6(a) shows that the axial stress σ xx distribution evaluated using 8H18 element match the analytical solution for a coarse mesh. Similarly,

the τxz , σ zz stress distributions obtained using 8H18 element are converging towards the analytical solution with increase of mesh den­ sity as shown in Fig. 6(b) and (c). In the case of 27H90 element, the axial stress σ xx results matches with the analytical solution for all the mesh sizes as shown in Fig. 7(a). In the case of coarse mesh, the τxz solution is discontinuous at the interface of 00 and 900 plies as shown in Fig. 7(b). With the increase of mesh density in the thickness direction, the discontinuities in the τxz solution vanish and match the analytical solution. Fig. 7(c) shows that the σ zz stress solution approaches the analytical solution for a coarse mesh. Thus from the numerical study, the hybrid elements are better in capturing the through-thickness stress distribution results with fewer elements when compared to the displacement brick elements. 3.2. Composite laminate subjected to sinusoidal pressure load A symmetric cross-ply [0/90/0] laminate subjected to sinusoidal � � � � pressure load b p t ¼ q0 sin πax sin πby as shown in Fig. 8 is studied. The 6

Thin-Walled Structures 148 (2020) 106587

L.B. Andraju and G. Raju

phenomenon. It is observed that all the stress distributions computed using the C3D8R, 8H18 and 27H90 brick elements matches with the analytical solution. However, the C3D8R element requires fine mesh when compared to the 8H18 and 27H90 elements. In the case of thin laminate (t=a¼1:100), the σ xx , τxz stresses are accurately captured by all the brick elements (see Fig. 10(a) and (b)). The σzz stress distribution evaluated using the 27H90 element is closer to the analytical solution. However, the 8H18 and C3D8R elements are not able to capture the σzz stress distribution in the thickness direction for a fine mesh as shown in Fig. 10(c). Accurate evaluation of σzz is necessary as it plays an impor­ tant role in the evaluation of mode-I and mixed-mode delamination failure initiation criteria in composite laminates. Thus, the hybrid ele­ ments are free from shear locking problems and can be used for accurate stress evaluation in thin as well as thick composite laminates. 3.3. Effect of mesh distortion Fig. 8. A simply supported composite plate subjected to sinusoidal pressure load on the top surface of the plate.

In this section, the effect of mesh distortion on the computation of stress solution using the hybrid brick elements is investigated. The mesh distortion is considered in the thickness direction as well as the in-plane direction of the laminate. Initially, the beam problem in Section 4.1 is considered with mesh distortion in the thickness direction as shown in Fig. 11. The dotted lines depict the intermediate nodes of the 27H90 element. Fig. 12 shows the σxx stress distribution on the top surface along the length of beam for the regular and distorted mesh of brick elements. In the case of 27H90 element, the stress solution computed is the same for both regular and distorted mesh. In the case of C3D8R and 8H18 ele­ ments, the mesh distortion has some influence on the stress solution. However, the distorted mesh results of 8H18 are better when compared with a fine regular mesh of C3D8R element. Next, the composite plate subjected to sinusoidal loading on the top

plate is simply supported on all the edges. The material properties of the lamina considered are same as the previous beam problem. Laminates with the thickness to length ratio of 1:100 and 1:10 are studied to un­ derstand the shear locking behavior of the hybrid brick elements. The plate is modeled with in-plane dimensions a ¼ b ¼ 1000 mm and thicknesses t ¼ 10 mm, 100 mm. The magnitude of the pressure load q0 ¼ 100 kPa is applied on the top surface of the laminate. Fig. 9 shows the normalized stress distribution results evaluated using the brick elements for a thick laminate (t=a ¼ 1:10) and compared with the Pagano’s analytical solution [34]. Full integration is performed for both the hybrid elements, whereas reduced integration is considered for the displacement brick elements to alleviate the shear locking

Fig. 9. Normalized through-thickness stress distribution in the laminate (t=a ¼ 1 : 10) subjected to sinusoidal loading: Comparison of the brick elements (C3D8R, 8H18, 27H90) result with Pagano’s analytical solution. a) σxx b) τxz c) σzz .

Fig. 10. Normalized through-thickness stress distribution in the laminate (t=a ¼ 1 : 100) subjected to sinusoidal loading: Comparison of the brick elements (C3D8R, 8H18, 27H90) result with Pagano’s analytical solution. a) σxx b) τxz c) σzz . 7

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 11. Composite beam with distorted mesh in the thickness direction.

surface (see Section 4.2) is studied for the distorted mesh in the plane of the laminate as shown in Fig. 13 with t=a ¼ 1 : 10. Normalized stresses � � � � σ xx x; 2b; t and σyy 2a; y; t distribution at the top surface of the laminate along the x and y-direction is shown in Fig. 14. In the case of C3D8R element, mesh sizes of 20 � 1 � 6 and 20 � 1 � 24 are considered and the stress solution is not influenced by the in-plane mesh distortion. Simi­ larly, the hybrid element results are also not influenced by mesh distortion. However, from Fig. 14 it is clear that C3D8R is under pre­ dicting the results with a fine mesh and 8H18 and 27H90 elements are good in capturing the stress results with a coarse mesh. Thus, the nu­ merical study shows that the hybrid element results are not influenced by the mesh distortion and is suitable for studying complex geometries with unstructured mesh distribution. 3.4. Effect of boundary layer near the supported edges In this section, the effect of boundary layer on the stress solution near the clamped edges is studied. A fine mesh is required to capture the effect of the boundary layer, and the behavior of the brick elements in capturing the stresses near the clamped ends is investigated. A multi­ layered composite beam clamped at both the ends with length ¼ 1000 mm and thickness ¼ 100 mm is studied. A uniformly distributed load of bt ¼ P b b ¼ 50 kPa is applied equally over the top and bottom surfaces of P

Fig. 12. Normalized axial stress σxx distribution at the top surface of the beam along the length for distorted and regular mesh.

the beam as shown in Fig. 15. The laminate stacking sequence and material properties details are given in Table 2, where Laminate 2 is a un-symmetric sandwich laminate. Material P [33] (same as used in Section 4.1) is used to model the composite lamina and material PVC is used to model the isotropic foam in the sandwich panel with properties of E1 ¼ E2 ¼ E3 ¼ 1:72 GPa, ν12 ¼ ν23 ¼ ν13 ¼ 0:3, G12 ¼ G13 ¼ G23 ¼ 0:663 GPa [5]. As there is no analytical solution available for this problem, the hybrid brick elements results are validated using the converged solution of C3D8R element given by Groh and Weaver [5]. They have used 96000 C3D8R elements by considering 160 elements along the thickness, 600 elements along the length and one element across the width. Plane strain condition is enforced in the width direction by using a single element and restraining the lateral edges from expansion. The mesh size of the hybrid elements is given in Table 2. The normalized stress fields σxx ,τxz and σzz are studied at four different locations 5%,10%,15% and 20% of the beam length from the clamped end. The through-thickness stress distribution results of Laminate 1 evaluated using the hybrid brick ele­ ments is shown in Fig. 16 and they match the converged displacement brick element solution presented in the Ref. [5]. Fig. 16(a) shows the through-thickness profile of σxx at different lo­ cations and the phenomenon of stress channeling is captured as reported in the reference [5]. The transverse shear stress τxz distribution (Fig. 16 (b)) corresponding to the 5% location from clamped edge shows that the shear stress is not maximum near to the center of the laminate. This phenomenon is attributed to the effect of the clamped boundary con­ ditions. However, at 20% location, τxz is maximum at the center of the laminate which reflects the absence of boundary layer effect on the

Fig. 13. Composite laminate with in-plane distorted mesh.

8

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 14. Normalized stress distribution at the top surface of the laminate (t=a ¼ 1 : 10) along the length and width direction for distorted and regular mesh. a) σxx b) σ yy .

sandwich panel. The total number of nodes in the case of C3D8R model is 1,93,522, whereas 51,858 and 77,787 for 8H18 and 27H90 elements respectively which is 27% and 40% of the C3D8R element model. Thus the hybrid elements are able to capture the boundary layer effect with less number of elements when compared to the displacement brick element. 3.5. VAT composite beam In this section, the behavior of VAT composite beam subjected to uniform pressure is studied using the hybrid brick elements. The VAT laminates considered with dimensions of length ¼ 250 mm and thick­ ness ¼ 12.5 mm (t=a ¼ 1:20) are pinned at both the edges. Table 3 summarizes the stacking sequence (taken from Ref. [5]) of the sym­ metric variable stiffness laminates VAT 1 and VAT 2 and the mesh size considered for different brick elements. The laminae are modeled using IM7 material with a thickness of 1.5625 mm. The material properties of IM7 is given as E1 ¼ 163 GPa, E2 ¼ E3 ¼ 12 GPa, ν12 ¼ ν23 ¼ ν13 ¼ 0:3, G12 ¼ 5 GPa, G13 ¼ 4 GPa, G23 ¼ 3:2 GPa. A uniformly distributed bt ¼ P b b ¼ 50 kPa is applied on the top and bottom surfaces. pressure of P

Fig. 15. A composite beam subjected to uniformly distributed load with clamped boundary conditions at both the ends.

stress results. Fig. 16(c) shows the behavior of σzz stress distribution at different locations and undergo changes in shape as we move away from the fixed end due to the boundary layer effect. Fig. 17 shows the through-thickness stress distribution results of sandwich laminate (Laminate 2) evaluated using hybrid brick elements. Hybrid brick element results are in good agreement with the converged results of C3D8R and are able to capture the boundary layer effects in the

Plane

strain

condition

is

enforced

in

the

width

direction.

Table 2 Composite laminate stacking sequence and mesh size information. Laminate

t=a

Thickness ratio

Material

Laminate 1

1:10

Laminate 2

1:10

[ð1=4Þ4 ]

[P4 ]

[ð1=8Þ2 =0:5]

[P2 =PVC=P2 ]

Stacking sequence

8H18

27H90

[0/90/0/90]

200 � 1 � 128

100 � 1 � 64

½0=90�

200 � 1 � 128

100 � 1 � 64

Fig. 16. Normalized through-thickness stress distribution in Laminate 1 with clamped boundary conditions at 5%, 10%, 15% and 20% location from the fixed end. a) σxx b) τxz c) σzz . (legend is same for all the plots). 9

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 17. Normalized through-thickness stress distribution in sandwich laminate (Laminate 2) with clamped boundary conditions at 5%, 10%, 15% and 20% location from the fixed end. a) σ xx b) τxz c) σ zz . (legend is same for all the plots). Table 3 VAT laminate stacking sequence and mesh size information. Laminate

Stacking sequence

VAT 1

½〈90j0〉=〈

VAT 2

½〈90j20〉=〈45j

90j0〉=〈45j 25〉=〈

45〉=〈 90j

45j45〉�s

20〉=〈

45j25〉�s

C3D8R

8H18

27H90

799 � 1 � 120

120 � 1 � 80

240 � 1 � 24

799 � 1 � 120

120 � 1 � 80

240 � 1 � 24

Fig. 18. Normalized through-thickness stress distribution in the VAT 1 laminate under pinned boundary condition. a) σ xx b)τxz c) σzz .

Fig. 19. Normalized through-thickness stress distribution in the VAT 2 laminate under pinned boundary condition. a) σ xx b)τxz c) σzz .

Through-thickness stress distribution of the VAT laminates are evalu­ ated using the hybrid brick elements and the stress results evaluated at the quarter span (x ¼ a=4) of the beam are shown in Figs. 18 and 19. The axial stress σxx computed for VAT 1 and VAT 2 laminate requires few hybrid elements to match the converged C3D8R solution (18(a) and 19(a)). The transverse shear stress τxz distribution computed using the hybrid element is slightly different from the C3D8R element in both the laminates (see Figs. 18(b) and 19(a)). In the case of VAT 1, the τxz dis­ tribution obtained using hybrid brick elements capture stress reversal in the outermost plies, whereas the C3D8R element fails to capture it. In

the case of VAT 2, the τxz stress reversal is captured by displacement as well as hybrid elements, but the magnitudes are different in the inner­ most plies. The difference is attributed to the use of reduced integration in displacement brick elements. The displacement brick element can attain the hybrid brick element results by using full integration. The transverse normal stress σzz distribution (Figs. 18(c) and 19(c)) computed using the 27H90 element match the converged C3D8R solu­ tion for both the VAT laminates. This numerical study shows that the hybrid brick elements are able to capture the through-thickness stress distribution of VAT laminates more accurately compared to 10

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 20. VAT laminate with thickness variation subjected to pressure on the bottom surface of the beam under clamped boundary conditions. a) VAT A: ½0 � 〈0j70〉�s b) VAT B: ½0 � 〈70j0〉�s .

Table 4 Tow-steered laminate (beam) stacking sequence and the mesh size information. Laminate

t=a

Material

Thickness ratio

Stacking sequence

C3D8R

8H18

27H90

VAT A

1:20

IM7

200 � 1 � 40

IM7

½0 � 〈0j70〉�s

400 � 1 � 80

1:20

[ð1=4Þ4 ]

799 � 1 � 120

VAT B

799 � 1 � 120

200 � 1 � 40

100 � 1 � 24

[ð1=4Þ4 ]

½0 � 〈70j0〉�s

Fig. 21. Normalized through-thickness stress distribution in the VAT laminate ½0 � 〈0j70〉�s without thickness variation under clamped boundary condition at different locations in the beam. a) σ xx b) τxz c) σ zz . (legend is same for all the plots).

Fig. 22. Normalized through-thickness stress distribution in the VAT laminate ½0 � 〈0j70〉�s with thickness variation under clamped boundary condition at different locations in the beam. a) σxx b) τxz c) σzz . (legend is same for all the plots).

displacement brick elements with less number of degrees of freedom.

clamped at both the ends, and the schematic of the thickness distribution of the VAT laminates is shown in Fig. 20. A uniform pressure of b b ¼ 100 kPa is applied on the flat bottom surface of the magnitude P

3.6. VAT laminated beam with variable thickness

laminate. The VAT laminate of length ¼ 250 mm and thickness ¼ 12.5 mm is considered for the numerical study and IM7 material properties are used for the lamina. Plane strain condition is enforced along the width direction. The stacking sequence and mesh size information of the displace­ ment and hybrid brick elements for VAT A and VAT B laminates are given in Table 4. For VAT A laminate (see Fig. 20(a)), the normalized stress distributions σxx , τxz and σ zz evaluated using the brick elements at

In this section, the effect of thickness variation on the stress solution of VAT laminate is investigated. Understanding the through-thickness stress distributions can provide insights into the design of VAT com­ posite laminate with thickness variation and their failure behavior under different loading conditions. In this section, the stress results of VAT laminate with thickness variation are compared with the same VAT laminate neglecting thickness variation. The VAT laminates studied are 11

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 23. Comparison of through-thickness stress contours of VAT laminated (½0 � 〈0j70〉�s ) beam evaluated using 27H90 element under clamped boundary condition with constant and variable thickness. a) σxx (constant thickness) b) σxx (variable thickness) c) τxz (constant thickness) d) τxz (variable thickness) e) σ zz (constant thickness) f) σ zz (variable thickness). (Comparison with C3D8R element is given in supplementary material).

Fig. 24. Normalized through-thickness stress distribution in the VAT laminate ½0 � 〈70j0〉�s without thickness variation under clamped boundary condition at different locations in the beam. a) σ xx b) τxz c) σ zz . (legend is same for all the plots).

Fig. 25. Normalized through-thickness stress distribution in the VAT laminate ½0 � 〈70j0〉�s with thickness variation under clamped boundary condition at different locations in the beam. a) σxx b) τxz c) σzz . (legend is same for all the plots).

different locations (5%, 20% and 50%) of the beam for constant thick­ ness is shown in Fig. 21. At 5% location from the clamped end, the axial stress σ xx changes from tension (bottom surface) to compression (top surface) across the thickness of the laminate as shown in Fig. 21(a). But as we move grad­ ually away from the clamped end towards the middle of the beam (20% from the clamped end), σ xx changes from compression at the bottom surface to tension at the top surface of the laminate (Bending behavior). The transverse shear stress τxz distribution is parabolic and their

magnitude change as we move from the fixed end as shown in Fig. 21(b). The magnitude of τxz distribution is high near the clamped ends and low near the mid-span of the beam. This behavior is attributed to the boundary layer effect on the stress distribution (σxx , τxz ). Fig. 21(c) shows that the transverse normal stress σ zz distribution remains the same at all the locations. Next, the VAT A laminate with thickness variation is considered. Through-thickness stress distribution computed using the brick elements at different locations along the span of the beam is shown in Fig. 22. 12

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 26. Comparison of through-thickness stress contours of VAT laminated (½0 � 〈70j0〉�s ) beam evaluated using 27H90 element under clamped boundary condition with constant and variable thickness. a) σxx (constant thickness) b) σxx (variable thickness) c) τxz (constant thickness) d) τxz (variable thickness) e) σ zz (constant thickness) f) σ zz (variable thickness). (Comparison with C3D8R element is given in supplementary material).

Fig. 27. VAT laminate with thickness variation subjected to pressure on the bottom surface of the beam under simply supported boundary conditions. a) VAT A: ½0 � 〈0j70〉�s b) VAT B: ½0 � 〈70j0〉�s .

Fig. 28. Comparison of through-thickness stress contours of VAT laminated (½0 � 〈0j70〉�s ) beam evaluated using 27H90 element under simply supported boundary condition with constant and variable thickness. a) σxx (constant thickness) b) σ xx (variable thickness) c) τxz (constant thickness) d) τxz (variable thickness) e) σ zz (constant thickness) f) σzz (variable thickness). (Comparison with C3D8R element is given in supplementary material).

The thickness variation and the curvature effect in the VAT laminate results in a totally different through-thickness stress distribution when compared to constant thickness laminate. At 5% from the clamped end, the axial stress σ xx on the top surface is compressive and changes to

tensile at the bottom surface (Fig. 21(a)). At 20% from the clamped edge the axial stress is compressive throughout the thickness of the laminate. Around the mid-span of the beam, the σ xx results exhibits the regular behavior of changing from tensile to compressive (top surface to bottom 13

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 29. Comparison of through-thickness stress contours of VAT laminated (½0 � 〈70j0〉�s ) beam evaluated using 27H90 element under simply supported boundary condition with constant and variable thickness. a) σxx (constant thickness) b) σ xx (variable thickness) c) τxz (constant thickness) d) τxz (variable thickness) e) σ zz (constant thickness) f) σzz (variable thickness). (Comparison with C3D8R element is given in supplementary material).

Fig. 30. VAT composite panel design with a stacking sequence of ½θ1 = (Reproduced from Ref. [16]).

θ1 =

θ2 =θ2 �s (a) Reference tow paths (b) fibre angle distributions across the half width

surface). The transverse shear stress τxz is high at the top surface of the laminate at 5% location from the clamped edge as shown in Fig. 22(b). At 20% from the clamped end, the τxz is non-zero near the top surface of the beam due to the curvature introduced by fiber angle change. The transverse normal stress σzz distribution shows high magnitude at the top surface of the laminate near to the clamped end (see Fig. 22(c). As we move away from the clamped end, the σ zz stress value is decreasing and attains zero value on the top surface at the mid-span of the beam. Stress contour plots of the VAT A beam with constant and variable thickness are shown in Fig. 23. The σxx stress distribution is symmetric about the z-axis and shows that the bending effect is observed near the mid-span of the beam and the boundary layer effect is observed near the clamped ends. The bending effect is observed in the beam where the fiber angle is close to 00 . The τxz stress distribution is maximum near the clamped ends and they are anti-symmetric about the z-axis. The σ zz distribution attains maximum variation near the clamped ends in the case of thickness variation and they are symmetric about the z-axis. The

stress distribution for the VAT A laminate without and with thickness variation are entirely different and their damage initiation and pro­ gression may be dissimilar. Next, VAT B laminate (½0 � 〈70j0〉�s ) is considered for the numerical study. In this case, the laminate is thick at the center and thin at the ends as shown in Fig. 20(b). The distribution of the through-thickness stress fields at different locations (5%, 20% and 50%) are evaluated using the brick elements and is shown in Figs. 24 and 25 for constant and variable thickness laminates, respectively. In the case of VAT B laminate with constant thickness, the σxx stress distribution shown in Fig. 24(a) shows that the magnitude of stress is high at the fixed ends where the fiber angle is 00 and low at the mid-span of the beam where the fiber angle is 700 . The τxz stress distribution is parabolic and its magnitude decreases as we move away from the clamped end as shown in Fig. 24(b). The σ zz stress distribution has considerable variation near to the clamped ends and as we move away from the edges, the variation is minimal as shown in Fig. 24(c). 14

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Fig. 31. FE analysis of VAT Laminate under compression. a) FE Model with Boundary conditions b) Validation of global force-displacement response obtained from FE model with experiments.

span of the beam as shown in Fig. 25(c). Fig. 26 shows the through-thickness stress contour plots for VAT B laminate with constant and variable thickness laminate evaluated using the 27H90 brick element. The σ xx stress distribution clearly depicts the reverse bending behavior due to boundary layer effect near the clamped ends for constant and variable thickness laminates. For constant thick­ ness laminate, the σxx stress distribution is showing bending behavior (tension on the top surface and compression on the bottom surface) at the mid-span of the beam. For variable thickness laminate, the σ xx stress distribution at the mid-span of the beam is of low magnitude. The τxz stress distribution is high near the clamped ends for constant and vari­ able thickness laminate. But, in the case of thickness variation, the transverse shear stress τxz attains localized high value near the top sur­ face as shown in Fig. 26(d). The transverse normal stress σ zz is symmetric about z-axis for variable thickness laminate. Thickness variation in the VAT A and VAT B introduces high compressive transverse normal stress σ zz in the thick region of the laminate which is not present in the case of constant thickness laminates. For the given loading direction, the σ zz is compressive and they mitigate mode I delamination. However, if the loading direction is reversed, the interlaminar transverse normal stresses σ zz leads to the mode I delamination near the thicker regions of the VAT laminate.

Fig. 32. Comparison of stress-strain results of VAT laminate under compression between experiments and FE simulation.

In the case of VAT B laminate with thickness variation, the σ xx stress distribution is high near the clamped edges, but away from the edges, the magnitude of stress is low as shown in Fig. 25(a). At 20% location from the clamped end, the σxx stress distribution is always positive (tension) throughout the thickness of the laminate which is uncharac­ teristic for a beam subjected to transverse load. At the mid-span of beam (x ¼ a=2), the σ xx distribution is changing from compression to tension from bottom surface to mid plane of the laminate and the distribution is changing again from tension to compression from mid plane to the top surface of the laminate. The transverse shear stress τxz distribution is not parabolic and non-zero shear is observed on the top surface of the laminate near the clamped ends as shown in Fig. 25(b). The transverse normal stress σ zz distribution is influenced by the thickness variation and exhibit high magnitude of compressive stress in the center plies at mid-

3.6.1. Effect of simply supported boundary conditions In this section, VAT A (½0 � 〈0j70〉�s ) and VAT B (½0 � 〈70j0〉�s ) lam­ inates with constant and variable thickness are evaluated under simply supported boundary conditions. A uniform pressure load is applied on the flat bottom surface of the laminate in the downward direction as shown in Fig. 27. The 27H90 brick element with a mesh size of 200 � 1 � 40 and 400 � 1 � 60 for VAT A and VAT B, respectively are used for evaluating the through-thickness stress distribution. Fig. 28 shows the through-thickness stress contour plots of VAT A laminate with constant and variable thickness configuration. In the case of simply supported boundary conditions, the results show the minimal influence of the boundary layer effect on the stress solution of VAT laminate. The axial stress σ xx distribution shows the bending behavior 15

L.B. Andraju and G. Raju

Fig. 33. Through-thickness stress contours of VAT laminate ½θ1 = 27H90 element. a) σ zz b) τxz c) τyz .

Thin-Walled Structures 148 (2020) 106587

θ1 =

θ2 =θ2 �s with constant thickness and thickness variation under compression evaluated using

regions of constant and variable thickness VAT laminate (see Fig. 28(a) and (b)). The transverse shear stress τxz is zero on the top and bottom surfaces of the beam in the case of constant thickness laminate, and the τxz distribution is high near the supported edges as shown in Fig. 28(c). For variable thickness laminate, τxz is non-zero on the top surface of the beam and the τxz distribution is high at locations away from the sup­ ported ends as shown in Fig. 28(d). The transverse normal stress σ zz distribution remains the same across the beam for constant thickness laminate. For variable thickness laminate, the σzz stress distribution is of tensile nature near the thick regions of the beam and might lead to

delamination failure. Fig. 29 shows the through-thickness stress distribution in VAT B laminate evaluated using the 27H90 element. The influence of thickness variation on the axial stress σ xx distribution of VAT B laminate is shown in Fig. 29(a) and (b). A non-zero transverse shear stress τxz of high magnitude is observed on the top surface of the beam as shown in Fig. 29 (d) in the case of variable thickness laminate. The transverse normal stress σ zz distribution is constant throughout the span of the beam except at the mid-span location for constant thickness laminate as shown in Fig. 29(e). In the case of thickness variation, interlaminar tensile stress 16

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

of high magnitude is observed in the center plies at the mid-span of the beam (see Fig. 29(f)). For the laminates (VAT A and VAT B) with thickness variation, the magnitude of the transverse stresses τxz and σ zz are high near the thick region of the laminate. Through-thickness stress results show the effect of thickness variation in the modeling of VAT laminates and also provides insight into the mechanics behind damage initiation and propagation.

with constant and variable thickness. The results also shows the differ­ ence in through-thickness stress distribution of VAT laminates with and without thickness variation along the width of the laminate. The stress distributions σ zz , τxz , and τyz are considered because they are responsible for the delamination initiation and propagation. It is reported in the reference [16], the delamination was observed near the edges, which supports the stress distribution evaluated using the FE models. The transverse normal stress σzz is maximum near the edges for the case of variable thickness VAT laminate and might lead to Mode-I delamination failure. But, for constant thickness VAT laminate, σzz attains maximum value away from the edges and highlights the importance of considering thickness variation in the FE model. The transverse shear stress τxz dis­ tribution is discontinuous with stress reversal observed near the edges for both the constant and variable thickness VAT laminate. The trans­ verse shear stress τyz variation is negligible, but exhibits variation on the edges x ¼ �a=2 which is shown in the isometric view in Fig. 33(c). Thus, the FE results of the CTS fabricated VAT laminates evaluated using hybrid brick elements are validated using experimental results. There­ fore, the hybrid brick elements can be used for accurate stress analysis of VAT laminates as well as progressive damage analysis under different loading conditions.

3.7. VAT laminate under compression In this section, the numerical performance of hybrid brick elements is validated with the experimental results available in the literature for the VAT laminates. In the Ref. [16,35], the experimental results of CTS manufactured VAT laminates under compression are presented. The VAT laminate stacking sequence considered for the study is ½θ1 = θ1 = θ2 =θ2 �s . Fig. 30(a) and (b) shows the fiber trajectories and the non-linear fiber angle variations of the VAT lamina, respectively. A plate with length of 150 mm along y-direction and width of 100 mm in the x-direction is considered. The VAT laminate with thickness variation is modeled with the reference lamina thickness of t0 ¼ 0.13 mm, and the variation in the thickness is evaluated using Eqn. (1). Here, the thickness at the two edges (x ¼ �a/2) of the laminate is � 4.1 mm and �1.2 mm at the middle of the laminate. The material properties of the lamina is given as E1 ¼ 163 GPa, E2 ¼ E3 ¼ 6:8 GPa, ν12 ¼ ν13 ¼ 0:28, ν23 ¼ 0:4 GPa, G12 ¼ G13 ¼ 3:4 GPa, G23 ¼ 2:4 GPa [16]. A uni­ form compressive displacement in the y-direction is applied at y ¼ b= 2 by fixing the other end at y ¼ b=2. Simply supported boundary conditions are applied at x ¼ �b=2 and are not constrained in the x-direction. The out-of-plane displacement along the edges y ¼ � a= 2 and z ¼ 0 is constrained. Fig. 31(a) shows the FE model and the boundary conditions of the VAT laminate under compression. The laminate is evaluated using C3D8R (Mesh-200 � 300 � 8), 8H18 (Mesh-200 � 300 � 8), and 27H90 (Mesh-100 � 150 � 8) elements. The global force-displacement response of the VAT laminate under compression obtained from the experiments and FE analysis are shown in Fig. 31(b) As the damage analysis is beyond the scope of this paper, the response in the linear regime is considered to validate the present FE model results. The global force-displacement response matches well for all the element models (see Fig. 31(b)). To compare the effect of thick­ ness variation, the global force-displacement of VAT laminate without thickness variation is also shown in Fig. 31(b). The stiffness of VAT laminate with thickness variation is �2.8 times higher than the stiffness of VAT laminate with constant thickness of the same configuration (fiber angle and material properties) and shows the importance of considering thickness variation in the FE model of VAT laminate. In Ref. [35], the compression strain field of the VAT laminate measured using the 3D digital image correlation (DIC) technique is presented. The DIC strains were captured corresponding to the area of � 64 mm � 115 mm. The compressive strain is obtained by taking average over the area of the DIC acquisition, and the stress is evaluated by dividing the load with the cross-sectional area of the laminate. The stress and strain from the FE model are evaluated using the same procedure and the comparison is shown in Fig. 32. It is reported in the reference [35], that the VAT laminate exhibits a non-linear behavior at around 0.14% of the applied strain without any damage in the panel. Fig. 32 shows the good correlation between the experimental results and the FE results until 0.3% of compressive strain and deviates beyond it. This is attributed to the linear analysis of FE models and the nonlinear geo­ metric strain terms being not considered in the FE analysis. Fig. 33 shows the through-thickness stress distribution evaluated using the 27H90 hybrid brick element at y ¼ b=2 for the VAT laminate

4. Conclusion In this paper, the effect of thickness variation on the 3D stress dis­ tribution of VAT laminates manufactured using the CTS technique has been studied using the hybrid brick elements. A UEL subroutine for the 8-node and 27-node hybrid brick elements has been implemented in ABAQUS finite element software. Benchmark problems in composites were studied and the 3D stress solution computed using the hybrid brick elements match the analytical solution requiring fewer elements than the displacement brick element. The hybrid element results were shown to be free from the effects of shear locking and mesh distortion. Also, the hybrid elements were able to capture the boundary layer effect in the stress solution requiring fewer elements. Through-thickness stress dis­ tribution results of VAT laminates manufactured using CTS are completely different from VAT laminates with constant thickness. The FE solution obtained using the hybrid brick elements were validated with the experimental results of VAT laminate under compression available in literature. The combined effect of the thickness variation as well as the fiber angle variation results in a different stress redistribution in the plane of the VAT laminate and also, the failure initiation behavior may be different from the constant thickness VAT laminate is reported. The stress results show the significance of considering manufacturing effects in the design of VAT laminates using CTS technique. In future, the hybrid brick element will be extended to perform progressive damage analysis of VAT composite laminates. Acknowledgements The authors gratefully acknowledges the Department of Science and Technology (DST) – India for partial funding of this research work through the Fast Track Scheme for Young Scientists (Project No: YSS/ 2015/000677). The authors would like to thank Dr. Syed Nizamuddin Khaderi, Assistant professor in the Department of Mechanical & Aero­ space Engineering, IIT Hyderabad for his help in the implementation of UEL subroutine in ABAQUS. The authors wish to thank Dr. Rainer Groh, Senior Research Associate, University of Bristol for providing the analytical solutions of benchmark problems.

17

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.tws.2019.106587. Appendix Stress interpolation Stress interpolation for the 8-node and 27 node hybrid brick elements are as follows: 8-Node hybrid brick element: Stress interpolation functions with 18 terms are given as [ [23],

σξξ ¼ β1 þ β2 η þ β3 ζ þ β4 ηζ; σηη ¼ β5 þ β6 ξ þ β7 ζ þ β8 ξζ; σζζ ¼ β9 þ β10 ξ þ β11 η þ β12 ξη; τξη ¼ β13 þ β14 ζ; τηζ ¼ β15 þ β16 ξ; τξζ ¼ β17 þ β18 η: 27-Node hybrid brick element: Stress interpolation functions with 18 terms are given as [27],

σξξ ¼ β1 þ β2 ξ þ β3 η þ β4 ζ þ β5 ξη þ β6 ηζ þ β7 ξζ þ β8 ξηζ þ β9 ξη2 þβ10 ξζ2 þ β11 ξηζ2 þ β12 ξη2 ζ þ β13 ξη2 ζ2 þ β14 η2 þ β15 ζ2 þβ16 η2 ζ þ β17 ηζ2 þ β18 η2 ζ2 ; σηη ¼ β19 þ β20 ξ þ β21 η þ β22 ζ þ β23 ξη þ β24 ηζ þ β25 ξζ þ β26 ξηζ þβ27 ξ2 η þ β28 ηζ2 þ β29 ξ2 ηζ þ β30 ξηζ2 þ β31 ξ2 ηζ2 þ β32 ξ2 þβ33 ζ2 þ β34 ξ2 ζ þ β35 ξζ2 þ β36 ξ2 ζ2 ; σζζ ¼ β37 þ β38 ξ þ β39 η þ β40 ζ þ β41 ξη þ β42 ηζ þ β43 ξζ þ β44 ξηζ þβ45 ξ2 ζ þ β46 η2 ζ þ β47 ξ2 ηζ þ β48 ξη2 ζ þ β49 ξ2 η2 ζ þ β50 ξ2 þβ51 η2 þ β52 ξ2 η þ β53 ξη2 þ β54 ξ2 η2 ; τξη ¼ β55 þ β56 ξ þ β57 η þ β58 ζ þ β59 ξη þ β60 ηζ þ β61 ξζ þ β62 ξηζ þβ63 ξζ2 þ β64 ηζ2 þ β85 ζ2 þ β88 ζ2 ξη; τηζ ¼ β65 þ β66 ξ þ β67 η þ β68 ζ þ β69 ξη þ β70 ηζ þ β71 ξζ þ β72 ξηζ þβ73 ξ2 η þ β74 ξ2 ζ þ β86 ξ2 þ β89 ξ2 ηζ; τξζ ¼ β75 þ β76 ξ þ β77 η þ β78 ζ þ β79 ξη þ β80 ηζ þ β81 ξζ þ β82 ξηζ þβ83 ξη2 þ β84 η2 ζ þ β87 η2 þ β90 ξη2 ζ: 0

The transformation matrix between the stresses with respect to the Cartesian and natural coordinate system is given by Ref. [26]. 1 0 1

σ ξξ σxx B σyy C B σ ηη C B C B C B σzz C B C B C ¼ T B σ ζζ C; B τxy C B τ ξη C B C B C @ τyz A @ τηζ A τξζ τxz where. 0

1

B a21 B B b2 B 1 B 2 B c1 T¼B Ba b B 1 1 B B b1 c1 B @a c 1 1

a22 b22 c22

a23 b23 c23

2a1 a2 2b1 b2

a2 b2

a3 b3

ða2 b1 þ a1 b2 Þ

b2 c2

b3 c3

ðb2 c1 þ b1 c2 Þ

a2 c2

a3 c3

ða1 c2 þ a2 c1 Þ

2c1 c2

C C C C C C 2c2 c3 2c1 c3 C; ða2 b3 þ a3 b2 Þ ða1 b3 þ a3 b1 Þ C C C ðb2 c3 þ b3 c2 Þ ðb1 c3 þ b3 c1 Þ C C ða3 c2 þ a2 c3 Þ ða3 c1 þ a1 c3 Þ A 2a2 a3 2b2 b3

2a1 a3 2b1 b3

where a1 ; b1 ; … are evaluated from the Jacobian matrix given by 0 ∂x ∂y ∂z 1 B ∂ξ ∂ξ B B ∂x ∂y B J¼B B ∂η ∂η B @ ∂x ∂y ∂ζ ∂ζ

∂ξ C 0 1 C a1 b1 c1 ∂z C C @ C ¼ a2 b2 c2 A: ∂η C a b c C

∂z A ∂ζ

3

3

3

18

L.B. Andraju and G. Raju

Thin-Walled Structures 148 (2020) 106587

References [17]

[1] M.W. Hyer, R. Charette, Use of curvilinear fiber format in composite structure design, AIAA J. 29 (6) (1991) 1011–1015. [2] M.W. Hyer, H. Lee, The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes, Compos. Struct. 18 (3) (1991) 239–261. [3] Z. Gurdal, R. Olmedo, In-plane response of laminates with spatially varying fiber orientations-variable stiffness concept, AIAA J. 31 (4) (1993) 751–758. [4] R. Groh, P. Weaver, A computationally efficient 2d model for inherently equilibrated 3d stress predictions in heterogeneous laminated plates. part ii: model validation, Compos. Struct. 156 (2016) 186–217. [5] R. Groh, P. Weaver, Deleterious localized stress fields: the effects of boundaries and stiffness tailoring in anisotropic laminated plates, Proc. R. Soc. A 472 (2194) (2016) 1–22. [6] Z. Wu, P.M. Weaver, G. Raju, B.C. Kim, Buckling analysis and optimisation of variable angle tow composite plates, Thin-Walled Struct. 60 (2012) 163–172. [7] Z. Wu, P.M. Weaver, G. Raju, Postbuckling optimisation of variable angle tow composite plates, Compos. Struct. 103 (2013) 34–42. [8] V. Oliveri, A. Milazzo, P.M. Weaver, Thermo-mechanical post-buckling analysis of variable angle tow composite plate assemblies, Compos. Struct. 183 (2018) 620–635. [9] A. Milazzo, V. Oliveri, Buckling and postbuckling of stiffened composite panels with cracks and delaminations by ritz approach, AIAA J. 55 (3) (2016) 965–980. [10] V. Gulizzi, V. Oliveri, A. Milazzo, Buckling and post-buckling analysis of cracked stiffened panels via an x-ritz method, Aero. Sci. Technol. 86 (2019) 268–282. [11] G. Raju, Z. Wu, P.M. Weaver, Buckling and postbuckling of variable angle tow composite plates under in-plane shear loading, Int. J. Solids Struct. 58 (2015) 270–287. [12] M.G. Günay, T. Timarcı, Stresses in thin-walled composite laminated box-beams with curvilinear fibers: antisymmetric and symmetric fiber paths, Thin-Walled Struct. 138 (2019) 170–182. [13] A. Beakou, M. Cano, J.B. Le Cam, V. Verney, Modelling slit tape buckling during automated prepreg manufacturing: a local approach, Compos. Struct. 93 (10) (2011) 2628–2635. [14] B.C. Kim, K. Potter, P.M. Weaver, Continuous tow shearing for manufacturing variable angle tow composites, Compos. Appl. Sci. Manuf. 43 (8) (2012) 1347–1356. [15] R. Groh, P. Weaver, Buckling analysis of variable angle tow, variable thickness panels with transverse shear effects, Compos. Struct. 107 (2014) 482–493. [16] T.D. Dang, S.R. Hallet, B.C. Kim, Y.L. Cahain, R. Butler, W. Liu, Modelling of as manufactured geometry for prediction of impact and compression after impact

[18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]

19

behaviour of variable angle tow laminates, J. Compos. Mater. 49 (12) (2015) 1423–1438. E.Q. Sun, Shear locking and hourglassing in msc nastran, abaqus, and ansys, in: Msc Software Users Meeting, 2006. M.A. Bhatti, Advanced Topics in Finite Element Analysis of Structures: with Mathematica and MATLAB Computations, John Wiley, New York, 2006. T.H. Pian, Derivation of element stiffness matrices by assumed stress distributions, AIAA J. 2 (7) (1964) 1333–1336. T. Pin, T.H. Pian, A variational principle and the convergence of a finite-element method based on assumed stress distribution, Int. J. Solids Struct. 5 (5) (1969) 463–472. T.H. Pian, D.P. Chen, Alternative ways for formulation of hybrid stress elements, Int. J. Numer. Methods Eng. 18 (11) (1982) 1679–1684. T.H. Pian, K. Sumihara, Rational approach for assumed stress finite elements, Int. J. Numer. Methods Eng. 20 (9) (1984) 1685–1695. T.H. Pian, P. Tong, Relations between incompatible displacement model and hybrid stress model, Int. J. Numer. Methods Eng. 22 (1) (1986) 173–181. S. Mau, P. Tong, T. Pian, Finite element solutions for laminated thick plates, J. Compos. Mater. 6 (2) (1972) 304–311. R. Spilker, S. Chou, O. Orringer, Alternate hybrid-stress elements for analysis of multilayer composite plates, J. Compos. Mater. 11 (1) (1977) 51–70. C.S. Jog, A 27-node hybrid brick and a 21-node hybrid wedge element for structural analysis, Finite Elem. Anal. Des. 41 (11–12) (2005) 1209–1232. C.S. Jog, Improved hybrid elements for structural analysis, J. Mech. Mater. Struct. 5 (3) (2010) 507–528. E. Reissner, On a variational theorem in elasticity, J. Math. Phys. 29 (1–4) (1950) 90–95. E. Reissner, On a certain mixed variational theorem and a proposed application, Int. J. Numer. Methods Eng. 20 (7) (1984) 1366–1368. C.S. Jog, P. Motamarri, An energy-momentum conserving algorithm for nonlinear transient analysis within the framework of hybrid elements, J. Mech. Mater. Struct. 4 (1) (2009) 157–186. L. Demasi, Mixed plate theories based on the generalized unified formulation. part v: Results, Compos. Struct. 88 (1) (2009) 1–16. ABAQUS/Standard User’s Manual, Version 2017, 2016. Simulia. N.J. Pagano, Exact solutions for composite laminates in cylindrical bending, J. Compos. Mater. 3 (3) (1969) 398–411. N.J. Pagano, Exact solutions for rectangular bidirectional composites and sandwich plates, J. Compos. Mater. 4 (1) (1970) 20–34. Y.M. Le Cahain, D.S. Ivanov, Tow-scale analysis of structures lacking an elementary representative element: application to variable angle tow composites, in: Proceedings of V ECCOMAS Thematic Conference on the Mechanical Response of Composites, University of Bristol, UK, 2015.