Reduced Order Modeling of Composite Laminates Through Solid-Shell Coupling

ABSTRACT:

Composite laminates display a complex mechanical behavior due to their microstructure, with a through-thickness variation of the displacement and stress fields that depends on the fiber orientation in each layer. Aiming to develop reduced-order numerical models mimicking the real response of composite structures, we investigated the capability and accuracy of finite element analyses coupling layered shell and solid kinematics. This study represents the first step of a work with the goal of accurately matching stress evolution in regions close to possible impact locations, where delamination is expected to take place, with reduced computational costs. Close to such locations, a 3-D modeling is adopted, whereas in the remainder of the structure, a less computationally demanding shell modeling is chosen. To test the coupled approach, results of numerical simulations are presented for a quasi-statically loaded cross-ply orthotropic plate, either simply supported or fully clamped along its boundary.

KEYWORDS:
Composite structures; Reduced-order modeling; Solid-shell coupling

INTRODUCTION

The numerical analysis of complex structures often requires very detailed space discretizations, especially in regions where high-field gradients are expected or where an enhanced solution accuracy is required for reliability issues. To speedup the analysis, reduced-order models can be formulated from a purely mathematical standpoint; alternatively, a wise coupling of refined and coarse-grained discretizations can be adopted. According to the notation proposed in Sellitto et al. (2011), a coupling of local (namely, accurate) and global (namely, time efficient) models can be envisaged.

The aforementioned problem can be relevant also in the analysis of composite laminates, which typically display a small thickness in comparison to their in-plane dimensions. Accordingly, a plate- or shell-like kinematics proves efficient to study problems mainly ruled by bending, if the structural response is of concern. If instead one focuses on decohesion mechanisms like delamination (Allix and Ladevèze 1992; Corigliano 1993; Abrate 1998; Schoeppner and Abrate 2000), a so-called meso-mechanical approach will be necessary. By allowing for the small thickness of the resin-enriched regions, simplified computational models resting on a lumping of the interlaminar decohesion onto 0-thickness surfaces have been developed in Corigliano and Mariani (2002), Corigliano et al. (2003), Mariani and Corigliano (2005), and Corigliano et al. (2006).

To model the decohesion processes, the through-thickness variation of the displacement field in displacement-based finite element (FE) simulations is required, through 3-D space discretizations. A fully solid discretization of the structural component would result in a heavy computational burden for real-life case studies. The other way around, a shell-like 2-D modeling proves sufficient for regions not exposed to delamination events. Here a hybrid approach is adopted, able to ensure accuracy with limited computational costs by coupling a solid kinematics to a shell kinematics away from the delaminating zones.

Dealing with ductile fracture processes in homogeneous metals, this approach was already adopted in Corigliano et al. (1999) to model the through-thickness propagation of a crack in pressurized pipelines, in Stringfellow and Paetsch (2009) to model the collision-induced failure of structural parts of a cab car, in Kim (2003) for the fatigue assessment of welded joints, and in Gong et al. (2016) to study the effect of the residual stress induced by welding on the buckling of storage tanks; in all the cases, stress triaxiality in the process regions plays a prominent role, and a shell-like model, basically missing the out-of-plane constraint on the stress field, does not necessarily prove accurate. In Reinoso et al. (2012), the envisioned approach was adopted to model the damage evolution in the skin-stringer joint of a complex composite sample for aeronautical applications (Li et al. 2013). In most of these studies, the commercial FE code Abaqus was adopted, and its solid-to-shell coupling feature was exploited to match the displacement fields at the boundaries of the facing solid and shell domains as well as to prevent unphysical kinking. This coupling is locally enforced by Abaqus in a strong form (Dassault Systems 2010); an alternative, variational enforcement was instead proposed in Blanco et al. (2008).

Focusing on coupled modeling for thin plates undergoing bending deformations induced by a distributed load, in this paper we provide an assessment of the hybrid approach in terms of plate deflection, elastic energy stored in the laminate, through-thickness variation of the stress fields, and speedup with respect reference solid models. Results are compared not only to those of fully solid 3-D simulations, but also to the outcomes (if available) of the Pagano approach (Pagano 1969, 1970) and of the first-order shear deformation theory (FSDT), (Reddy 2002; Wang et al. 2000).

In what follows, we highlight some computational features of the coupled procedure in the section "Solid-Shell Coupling and Computational Issues". In the subsequent section, results are discussed for a [90/0]_6s composite plate, either simply supported or fully clamped along its boundary. Finally, some concluding remarks are gathered.

SOLID-SHELL COUPLING AND COMPUTATIONAL ISSUES

Let us consider a flat rectangular plate, as depicted in Fig. 1. The plate is either simply-supported or fully-clamped along its boundary and is quasi-statically loaded over its top surface by a uniform transversely-distributed load q₀ , which is positive if pointing outward (the solution can be obviously generalized to the non-uniform load case). An orthonormal reference frame is introduced, with x and y axes located on the mid-plane of the plate and z axis pointing in the same direction of the load.

Figure 1
Plate bending problem, coupled solid-shell modeling: geometry and notation.

The central portion of the plate, whose side lengths are s_a and s_b , is modeled with a full 3-D kinematics (Fig. 1); the remainder of the plate is instead modeled with a 2-D shell-like kinematics. Accordingly, the former region is discretized with 3-D brick elements, whereas the latter region is discretized with 2-D shell elements. Numerical simulations here collected have been run using the commercial FE code Abaqus (Dassault Systems 2010).

In the remained of this section, we discuss 2 topics that affect the accuracy of the solutions provided by the hybrid 2- and 3-D model: the elements' kinematics as well as the coupling between the solid and shell structural regions.

As far as the shell kinematics is concerned, restricting the analysis to the small displacement regime, the 4-node (reduced-integration) S4R element has been adopted. The through-thickness variation of the strain and stress fields featured by such element is compliant with the FSDT; this kinematics implies that segments normal to the mid-plane of the plate retain their straight shape when deformed. Even if shear correction factors are adjusted in the case of anisotropic materials, and an inhomogeneous material response can be allowed for along the thickness direction, the transverse shear deformations are assumed constant throughout the whole thickness by this element kinematics. This assumption represents an approximation to the real composite behavior, which is instead characterized by the so-called zigzagging of the stress and strain fields (Bogdanovich and Pastore 1996; Reddy 2002; Carrera 2003). The adopted reduced-integration is developed based on an assumed strain approach via Hu-Washizu variational principle plus stabilization to avoid 0-energy deformation patterns and shear locking in case of thin plates (Bathe and Dvorkin 1984; Simo et al. 1989; Dassault Systems 2010).

In the solid region, linear brick elements have been adopted. Due to the considered bending-dominated deformation of the plate, the C3D8I elements featuring incompatible modes have been selected to avoid, or reduce as much as possible, parasitic shear and volumetric locking effects (Simo and Rifai 1990).

The coupling between solid and shell regions of the plate has been obtained through a surface-based interaction. This technique allows matching the displacements of nodes on the border of the solid region to the displacements and rotations of nodes on the border of the shell region, so as not only displacement jumps but also rotation jumps (i.e. kinks) are locally prevented (Fig. 2). Because of the S4R shell kinematics mentioned above, which cannot model the zigzagging of the fields along the through-thickness direction, this coupling introduces a local perturbation in the numerical solution. To assess the effects of such perturbation on the modeled structural response, in the next section the size ratios α = s_a /a and β = s_b /b (Fig. 1) are varied in the range of 0.25 ≤ α, β ≤ 0.75.

Figure 2
Example of the deformed confi guration of a bent composite plate modeled with the hybrid solid-shell approach (the plate defl ection has been magnifi ed to highlight the quality of coupling along the border between the solid and shell regions). Solid and shell elements are respectively depicted using light and dark grey colors.

RESULTS AND DISCUSSION: [90/0]_6S LAMINATE

A symmetric cross-ply [90/0]_6s laminate is considered; the side lengths of the composite plate are a = b = 127 mm, while the thickness is h = 3.36 mm. A uniform transversely-distributed load q₀ = 2 MPa is adopted. The elastic properties of each transversely isotropic layer, in a local x₁, x₂, x₃ reference frame, with axis xaligned with the fiber longitudinal direction, are: E₁₁ = 144,000 MPa, E₂₂ = E₃₃ = 9,690 MPa, G₁₂ = 55,385 MPa, G₁₃ = G₂₃ = 5,760 MPa, and ν₁₃ = ν₂₃ = 0.3 (Salerno 2009). In what precedes, E is a Young's modulus, G represents a shear modulus, and ν is a Poisson's ratio.

A preliminary investigation of the effects of meshing on the results has been carried out; what turned out is that the stacking sequence, and the relevant zigzagging of the in-plane strain and stress fields require a meshing with more than 1 element across the thickness of each layer to approach the FSDT solution (Salerno 2009). Anyhow, results are here presented in the case of 1 element discretizing each layer in the through-thickness direction to assess the accuracy of the numerical solutions when minimal meshing and, therefore, minimal computing time are enforced. To assess the relevant accuracy, the outcomes of the simulations will be compared with those obtained with a so-called overkill model featuring 4 elements to discretize across the thickness each layer of the composite panel; such model has been checked to provide mesh-independent results, hence no benefits would be obtained by further refining the space discretization in the through-thickness direction and over the mid-plane of the plate.

As far as the simply supported plate case is concerned, results are reported in Fig. 3 in terms of the through-thickness variation of the stress fields for α = β = 0.5; while stress components σ_xx and σ_yy are provided at the plate center, the shear component σ_xy is given as measured at the plate corner due to the symmetry in the solution (for cross-ply laminates) leading to a vanishing σ_xy close to the center. While the fully solid discretization provides results almost perfectly matching the FSDT ones, at least at the plate center, the solid-shell coupling introduces a perturbation, leading to reduced stress amplitudes close to the top and bottom surfaces of the plate. Since σ_xy is given at the plate corner, it must be noted that with the coupled solid-shell kinematics only 2 points are shown for this field in Fig. 3, as provided by the FE code for shell elements.

Figure 3
Simply supported laminate. Comparison among Pagano and FSDT solutions as well as results of the simulations in terms of: (a) σ_xx and (b) σ_yy at the plate center; (c) σ_xy at the plate corner.

Results in terms of central plate deflection and overall stored strain energy are shown in Fig. 4 to understand the effect of the modeling parameters α and β. The strain energy stored in the whole plate is here adopted to implicitly assess the in-plane size of the region affected by the distortion in the stress and strain fields due to the coupling of the solid and shell kinematics, keeping in mind that the zigzagging of the fields is not enforced in the computational model. In the investigated intervals, it is shown that the impact of α and β is marginal: the maximum discrepancy with respect to the reference overkill simulation (whose results are represented by the dashed lines in the graphs) amounts to 2.6% as for the central deflection and 7.6% as for the overall strain energy.

Figure 4
Simply-supported laminate. Effect of coupling coefficients α (with β = 0.5) (a and c) and β (with α = 0.5) (b and d) on: (a and b) central plate deflection; (c and d) stored strain energy. Dashed lines in the graphs represent the results of the reference overkill simulation.

If the laminate is instead assumed to be fully clamped along its boundary, as before a comparison, in terms of the through-thickness variation of the stress field, it is reported in Fig. 5 between the fully-solid and the coupled solid-shell (for α = β = 0.5) models featuring the minimal discretization in the out-of-plane direction, with a single element for each lamina. In this example, the coupled model provides a slight overestimation of the stress level at the center of the plate, amounting to around 1.8% in terms of the extreme values at the top and bottom surfaces of the plate as well as an underestimation at the plate corner, amounting to about 14.6% at the same top and bottom surfaces. As in Fig. 3, in these plots the stress components σ_xx and σ_yy are provided for all the nodes across the plate thickness in the region modeled with solid elements; the stress component σ_xy at the plate corner is instead provided for the coupled solid-shell model only at the top and bottom plate surfaces, as furnished by the FE code.

Figure 5
Clamped laminate. Comparison between results of the simulations, either allowing for solid-shell coupling or adopting a uniform 3-D space discretization in terms of: (a) σ_xx and (b) σ_yy at the plate center; (c) σ_xy at the plate corner.

Figure 6 shows the results in terms of central deflection and overall stored strain energy, compared to the reference overkill model (whose results are again represented by the dashed lines in the graphs). Even with a fully-clamped boundary, the α and β ratios do not affect much the solution; at most, a discrepancy of 5.4% with respect to the fine mesh solution is reported.

Figure 6
Clamped laminate. Effect of coupling coefficients α (with β = 0.5) (a and c) and β (with α = 0.5) (b and d) on: (a and b) central plate deflection; (c and d) stored strain energy. Dashed lines in the graphs represent the results of the reference overkill simulation.

Figures 4 and 6 show that parameters α and β might have a different effect on the accuracy of the solution, especially in terms of the elastic energy stored in the whole plate. This can be basically linked to the transversely isotropic behavior of each lamina, to the panel lay-up and so to the mechanical interaction among laminae when the plate is bent; since α and β move the boundary between the solid and shell regions (the latter is not able in the present model to capture the actual zigzagging of the stress and strain fields) along the 2 orthogonal in-plane directions, they can have a different impact on the solution, as shown here.

To understand the computational gain achieved through the hybrid kinematics, which looks necessary to motivate the adoption of a complicated modeling of the composite plates, outcomes are gathered in Fig. 7 in terms of the speedup as a function of the inverse of the relative in-plane size of the solid region, i.e. , for both the considered boundary conditions. The speedup is measured as the ratio between the CPU time of the fully-solid analysis and the CPU time of the hybrid one. The same in-plane discretization has been adopted in the 2 aforementioned solutions; in the hybrid analyses, the same through-thickness discretization of the fully-solid simulation has been adopted at the central plate portion. As a term of assessment, the dashed line in the plot shows how the speedup would scale if linearly depending on (αβ)⁻¹; as the computing time depends also on the number of layers in the composite, having used 1 element across each lamina in the considered minimal meshing strategy, such linear dependence is here reported only as a qualitative, although rather good, estimation of the computational gain attained with the hybrid kinematics. These results have been obtained by running Abaqus on a laptop with Windows 7-64 bit as operating system and an Intel Core I7-2620 M @ 2.70 GHz as CPU.

Figure 7
Effect of the relative in-plane size (αβ)^-1 of the solid region on the speedup.

CONCLUSIONS

In this paper, a numerical investigation to assess the accuracy and efficiency of coupled solid-shell modeling for composite laminates has been provided.

Quasi-static or low-velocity impact loadings may lead to delamination processes, located along the interlaminar regions; hence, standard shell elements cannot be used in finite element modeling, as they do not account for jumps in the displacement field across the thickness. To capture the post-impact response of composite structures for aeronautical applications, a full 3-D kinematics for the plate region around the impact location looks necessary; in the remainder of the structure, where delamination is not expected to occur, a shell kinematics can instead prove sufficient.

In case of thin laminated plates, we have accordingly approached the problem by coupling solid and shell elements. For a [90/0]_6s composite plate under a uniform distributed load, it has been shown that the through-thickness variation of the stress field is matched by the hybrid modeling with a rather good level of accuracy: a discrepancy of at most 10 - 15%, with respect to reference numerical solutions, has been reported. It has been also shown that the speedup (computed as the ratio between the CPU time required to run the 3-D solid analysis and the CPU time required to run the hybrid model) basically scales with the ratio between the volume of the whole composite and the volume of the 3-D modeled region. In case of real-life composite structures, such approach is therefore expected to provide an excellent reduction of the computational costs.

Maximum computational gain can therefore be attained if the 3-D region is assumed large enough to enclose all the possible delaminating regions. Since the extent of delamination cannot be known a priori, being dependent on the loading conditions, a mesh updating procedure would be also designed so as the kinematics is locally switched from a 2-D shell one to a 3-D solid one every time a critical stress threshold is approached throughout the laminate thickness.

REFERENCES

• Abrate S (1998) Impact on Composite Structures. Cambridge: Cambridge University Press.
• Allix O, Ladevèze P (1992) Interlaminar interface modelling for the prediction of delamination. Compos Struct 22(4):235-242. doi: 10.1016/0263-8223(92)90060-P
  » https://doi.org/10.1016/0263-8223(92)90060-P
• Bathe KJ, Dvorkin EN (1984) A continuum mechanics-based four-node shell element for general non-linear analysis. Eng Computation 1(1):77-88. doi: 10.1108/eb023562
  » https://doi.org/10.1108/eb023562
• Blanco PJ, Feijoo RA, Urquiza SA (2008) A variational approach for coupling kinematically incompatible structural models. Comput Meth Appl Mech Eng 197(17-18):1577-1602. doi: 10.1016/j.cma.2007.12.001
  » https://doi.org/10.1016/j.cma.2007.12.001
• Bogdanovich A, Pastore C (1996) Mechanics of textile and laminated composites with applications to structural analysis. London: Chapman & Hall.
• Carrera E (2003) Historical review of Zig-Zag theories for multilayered plates and shells. Appl Mech Rev 56(3):287-308. doi: 10.1115/1.1557614
  » https://doi.org/10.1115/1.1557614
• Corigliano A (1993) Formulation, identification and use of interface elements in the numerical analysis of composite delamination. Int J Solid Struct 30(20):2779-2811.
• Corigliano A, Maier G, Mariani S (1999) Analysis of ductile fracture in damaged pipelines by a geometric parameter method. Eng Struct 21(10):924-936. doi: 10.1016/S0141-0296(98)000406
  » https://doi.org/10.1016/S0141-0296(98)000406
• Corigliano A, Mariani S (2002) Identification of a constitutive model for the simulation of time-dependent interlaminar debonding processes in composites. Comput Meth Appl Mech Eng 191(1718):1861-1894. doi: 10.1016/S0045-7825(01)00357-7
  » https://doi.org/10.1016/S0045-7825(01)00357-7
• Corigliano A, Mariani S, Pandolfi A (2003) Numerical modeling of rate-dependent debonding processes in composites. Compos Struct 61(1-2):39-50. doi: 10.1016/S0263-8223(03)00030-8
  » https://doi.org/10.1016/S0263-8223(03)00030-8
• Corigliano A, Mariani S, Pandolfi A (2006) Numerical analysis of rate-dependent dynamic composite delamination. Compos Sci Tech 66(6):766-775. doi: 10.1016/j.compscitech.2004.12.031
  » https://doi.org/10.1016/j.compscitech.2004.12.031
• Dassault Systems (2010) Abaqus 6.10 Theory Manual. Vélizy-Villacoublay: Abaqus Software.
• Gong J, Yu L, Wang F, Xuan F (2016) Effect of welding residual stress on the buckling behavior of storage tanks subjected to harmonic settlement. J Pressure Vessel Technol 139(1):011401011401-9. doi: 10.1115/1.4033941
  » https://doi.org/10.1115/1.4033941
• Kim D (2003) Welding simulation of ship structures using coupled shell and solid volume finite elements (Master's thesis). Bethlehem: Lehigh University.
• Li D, Qing G, Liu Y (2013) A layerwise/solid-element method for the composite stiffened laminated cylindrical shell structures. Compos Struct 98:215-227. doi: 10.1016/j.compstruct.2012.11.013
  » https://doi.org/10.1016/j.compstruct.2012.11.013
• Mariani S, Corigliano A (2005) Impact induced composite delamination: state and parameter identification via joint and dual extended Kalman filters. Comput Meth Appl Mech Eng 194(5052):5242-5272. doi: 10.1016/j.cma.2005.01.007
  » https://doi.org/10.1016/j.cma.2005.01.007
• Pagano NJ (1969) Exact solutions for composites laminates in cylindrical bending. J Compos Mater 3(3):398-411. doi: 10.1177/002199836900300304
  » https://doi.org/10.1177/002199836900300304
• Pagano NJ (1970) Exact solutions for rectangular bidirectional composites and sandwich plates. J Compos Mater 4:20-34.
• Reddy JN (2002) Energy principles and variational methods in applied mechanics. New York: John Wiley and Sons.
• Reinoso J, Blázquez A, Estefani A, París F, Cañas J, Arévalo E, Cruz F (2012) Experimental and three-dimensional global-local finite element analysis of a composite component including degradation process at the interfaces. Compos B Eng 43(4):1929-1942. doi: 10.1016/j.compositesb.2012.02.010
  » https://doi.org/10.1016/j.compositesb.2012.02.010
• Salerno G (2009) Damage analysis of composite laminates subject to low velocity impacts (PhD thesis). Milan: Politecnico di Milano.
• Schoeppner GA, Abrate S (2000) Delamination threshold loads for low velocity impact on composite laminates. Compos Appl Sci Manuf 31(9):903-915. doi: 10.1016/S1359-835X(00)00061-0
  » https://doi.org/10.1016/S1359-835X(00)00061-0
• Sellitto A, Borrelli R, Caputo F, Riccio A, Scaramuzzino F (2011) Methodological approaches for kinematic coupling of non-matching finite element meshes. Procedia Engineering 10:421-426. doi: 10.1016/j.proeng.2011.04.071
  » https://doi.org/10.1016/j.proeng.2011.04.071
• Simo JC, Fox DD, Rifai MS (1989) On a stress resultant geometrically exact shell model. Part III: Computational aspects of the nonlinear theory. Comput Meth Appl Mech Eng 79(1):21-70. doi: 10.1016/0045-7825(90)90094-3
  » https://doi.org/10.1016/0045-7825(90)90094-3
• Simo JC, Rifai MS (1990) A class of mixed assumed strain methods and the method of incompatible modes. Int J Numer Meth Eng 29(8):1595-1638. doi: 10.1002/nme.1620290802
  » https://doi.org/10.1002/nme.1620290802
• Stringfellow R, Paetsch C (2009) Modeling material failure during cab car end frame impact. Proceedings of the 2009 ASME Joint Rail Conference; Pueblo, USA.
• Wang CM, Reddy JN, Lee KH (2000) Shear deformable beams and plates. New York: Elsevier Science.

Publication Dates

• Publication in this collection
  26 June 2017
• Date of issue
  Jul-Sep 2017

History

• Received
  19 July 2016
• Accepted
  21 Jan 2017