On the modeling and design of composite multilayered structures using solid-shell finite element model

Finite Elements in Analysis and Design 70–71 (2013) 1–14

Contents lists available at SciVerse ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

On the modeling and design of composite multilayered structures using solid-shell finite element model H. Naceur a,n, S. Shiri b, D. Coutellier a, J.L. Batoz c a

Lab. LAMIH, UMR8201 CNRS, Université Lille Nord de France, 59313 Valenciennes, France Lab. I2M, UMR 5295 CNRS, Arts & Métiers ParisTech, Université Bordeaux I, 33607 Pessac, France c Lab. Roberval, UMR 6253 CNRS, Université de Technologie de Compiègne, 60205 Compig`ne, France b

art ic l e i nf o

a b s t r a c t

Article history: Received 19 July 2012 Received in revised form 31 December 2012 Accepted 20 February 2013 Available online 1 April 2013

In this investigation a coupling between a 3D solid-shell element for the analysis of multilayered composite shell structures and a specific response surface method is proposed. The first part of the paper is dedicated to the finite element formulation of a developed composite 8-node solid-shell element called SCH8γ7, based only on translational degrees of freedom. The basis of the present finite element formulation is the standard 8-node brick element with tri-linear shape functions. A particular attention is given to alleviate shear, trapezoidal and thickness locking, without resorting to the classical plane-stress assumption. Assumed natural strain method and enhanced assumed strain method are used to improve the relative poor element behavior of a standard hexahedral displacement element. The anisotropic material behavior of layered shells is modeled using a fully three dimensional elastic orthotropic material law in each layer, including the thickness stress component. The second part of the paper will focus on an adaptive response surface method for the structural optimization problem. The response surfaces are built using moving least squares approximations and design of experiments by means of a specific method called Diffuse Approximation. Several numerical applications to composite multilayered shell structures are studied to show the applicability and effectiveness of the proposed procedure. Good results of analysis and optimization using the developed SCH8γ7 solid-shell element have been obtained in comparison with reference analytical solutions and with those obtained using the SC8R solid-shell finite element available in Á ABAQUS code. & 2013 Elsevier B.V. All rights reserved.

Keywords: Composite structures Solid-shell Optimization Design of experiment

1. Introduction In the modeling of shell structures various problems may appear especially when shell finite elements are used in combination with solid elements. Therefore special connection elements are necessary to link shell elements with solids having different degrees of freedom. It turns obvious to develop general-purpose brick elements, which are able to deal with any type of structures (solid, shell, and/or their combination). 3D-continuum elements are a variety of Finite Element (FE) models halfway between solid elements and thin shells. They have the same freedom configuration of solid elements but account for shell-like behavior in the thickness direction. They are very attractive for modeling shell-like regions of a 3D structure without the need of special elements to connect solid elements to shell nodes. 3D-continuum elements also called solid-shell elements

n

Corresponding author. Tel.: þ33 3 2751 1412; fax: þ 33 0 3 2751 1316. E-mail address: [email protected] (H. Naceur).

0168-874X/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.finel.2013.02.004

have many benefits compared to the degenerated shell elements, mainly because of the simplicity of their kinematics, their ability in modeling complex structures generally composed of massive and thin-walled regions and also special rotations treatment in geometric nonlinear analysis may be avoided. The 3D-continuum FE concept has attracted many researchers since the late 1990s. Domissy et al. [15] and Sze et al. [37] are probably the first authors to investigate a solid-shell approach for the analysis of plates and shells. In 1998, Hauptmann and Schweizerhof [20] propose an extension to the original solid-shell theory called double-node models where the position of an arbitrary point in the element is assumed to depend either linearly or quadratically, on the thickness coordinate. In 1998, Cho et al. [12] proposed a solid-shell element model based on the assumed strain formulation for buckling and postbuckling analysis of shell structures. During the last decade, solid-shell FE models for thin shell structures have attracted considerable attention. Numerous models have been proposed in the literature and can be found in Sze and Yao [38], Vu-Quoc and Tan [40], Areias et al. [3], Harnau and

2

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

Schweizerhof [19], Alves de Sousa et al. [2], Hannachi et al. [18], Quy and Matzenmiller [32], Nguyen et al. [29], Shiri et al. [35] and more recently the work of Schwarze and Reese [34], Moreira et al. [25]. Solid-shell element properties make them appropriate also for the modeling of laminated structures as can be seen in some representative references [3,21,25,33,40], to name just a few. In all of these works it has been shown that it is possible to use a solid-shell finite element models which possess no rotational degrees of freedom and obtain good results for the resolution of various structural problems. However, development of these elements is not straightforward, transverse shear, trapezoidal and thickness locking phenomena must be alleviated and therefore special treatments have to be included to suppress the numerical locking effects become crucial. Transverse shear locking is characterized by an overestimation of stiffness associated to transverse shear strain energy [2]. The most common methods that have been used to solve the shear locking effect are the Selective Reduced Integration (SRI) scheme [2] and also the Assumed Natural Strain (ANS) which has been applied firstly on shells by Dvorkin and Bathe [17] and as can be seen in the works of [3,11,15,35] for a fully integrated solid-shell element. Trapezoidal locking is only found in structures where the directors of the element edges are not perpendicular to the midplane. One method to resolve this pathology is by using the ANS interpolation of the transverse strain in thickness direction as performed in [8,9,37,38], where it was proposed to avoid artificial thickness straining. Thickness locking, caused by Poisson's ratio coupling to the in-plane and transverse normal stress and normal strain responses [14,15] can be overcome by plane-stress assumption, or more generally by the Enhanced Assumed thickness Strain (EAS), where a 7th parameter is added [10,21] and suppressed by condensation technique. This locking pathology usually appears when a solidshell element with only translational degrees of freedom is to be used in problems involving bending. The 7th parameter as proposed by Büchter and Ramm [10] is an extension of classical shear deformation theory. The interest of this model is especially useful when full three-dimensional constitutive law is used which allows solving problems involving large strains. In their work, Büchter and Ramm [10] describe the 7th parameter model along with a FE formulation and they introduce it on the element level by means of the hybrid-mixed formulation. In the present solid-shell model, the 7th parameter is resulting independently from the FE formulation, i.e. it can be understood as semi-discretization of the solid-shell element through the thickness. Based on this technique, the obtained 7th parameter model is considered simply as a two-dimensional continuous concept with 7 degrees of freedom per node of the reference surface. The 7th parameter is then eliminated using the static condensation procedure on the element level. Another alternative, has been proposed recently by Professor Brunet and his collaborators [4], concerning a solid-shell element with 9 nodes: 8 are located at the element vertices and the 9th is located at the centroid. The authors used classically a reduced integration with one in-plane quadrature point with an assumed shear-strain field to avoid locking phenomena. The centroid node is used as an extra parameter to enhance the displacement in the thickness direction and thus obtain a linear normal strain, allowing the use of full three-dimensional constitutive strain–stress behavior. During last few years, Response Surface Method (RSM) [24] has gained more and more importance in the optimization of general shell structures [30]. RSM has the advantage of replacing a complex response model by an approximate one based on results calculated at various points in the design space. The optimization

is then performed at a lower cost over such response surfaces. Two important issues when applying RSM to a particular problem concern the Design of Experiments (DOE) and construction of accurate approximation functions [13,22] so that rapid convergence may be achieved. In the present work, we exploit the RSM based on Diffuse Approximation (DA) [7,26] and particularly the notion of pseudoderivative to design a specific optimization technique, custom built for this regression model. This new method is an extension of pattern search in two aspects: (1) accommodation of arbitrary regular and irregular patterns; (2) design points eligible for inclusion in any pattern instance belong to a predefined set in the design where experiments are authorized. The outline of the paper is as follows. In Section 2 we present the kinematics, the material law, and the variational equations of a 3D-shell formulation with convective coordinates with the associated FE description. A special technique, which is necessary for integration through the layers in the 3D-case, is also given. In Section 3 the optimization procedure dedicated for the design of material parameters is given, a special attention is set to the built of response surface model, we use to approximate the implicit objective function. In the fourth section we discuss the applicability of the proposed procedure which consists in coupling between RSM, DOE and solid-shell solution for the modeling of several standard benchmarks which have been given by NAFEMS [28]. The obtained results are compared with those obtain using Á ABAQUS commercial software via its solid-shell element called SC8R. Finally some conclusions are drawn.

2. FE formulation of the multilayered solid-shell element 2.1. Kinematics To overcome the difficulties of rotational degrees of freedom (dof) in conventional shell elements, the shell kinematics of deformation is described using the position vectors of a pair of material points at the top and at the bottom of the shell midsurface (Fig. 1). In this kinematic description, a straight transverse fiber before deformation remains straight after deformation. Such transverse fiber does not need necessarily to be normal to the shell midsurface before deformation, as well as after deformation. With respect to nodal designation (Fig. 2), the coordinate vector X and displacement vector u of the element, are [5,6,16]
 4 1−ζ − 1 þ ζ þ Xðξ,η,ζÞ ¼ X 0 ðξ,ηÞ þ ζX n ðξ,ηÞ ¼ ∑ N i ðξ,ηÞ Xi þ Xi ð1Þ 2 2 i¼1 4

uðξ,η,ζÞ ¼ u0 ðξ,ηÞ þ ζun ðξ,ηÞ ¼ ∑ N i ðξ,ηÞ i¼1




1−ζ − 1 þ ζ þ u þ ui 2 i 2

Fig. 1. An 8-node solid-shell element.

 ð2Þ

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

3

Using Eqs. (2) and (3), the displacement gradient vector in the curvilinear coordinate system can be expressed as du ¼ L ds

ð12Þ

with L ¼ F T Lζ

ð13Þ

and du ¼ Lζ dξ is the displacement gradient tensor in the parametric coordinate system. In order to calculate the strain tensor directly in the curvilinear coordinate system, we have to define a tensor which is given by [6,15,41] C ¼ F −1 Q

ð14Þ

where 2 1 aζ  t 1ζ 6 2 a t C ¼6 4 ζ 1ζ 0

Fig. 2. Position vectors and covariant basis.

a1ζ  t 2ζ a2ζ  t 2ζ 0

a1ζ  nζ

3

7 a2ζ  nζ 7 5 a3ζ  nζ

ð15Þ

where Ni are the two-dimensional 4-node Lagrangian interpolation functions, X −i , u−i and X iþ , uiþ are respectively, the coordinate and displacement vectors of the ith node on the bottom and top shell surfaces respectively (Fig. 1). Using Eq. (1), the gradient vector dX can be expressed in terms of dξ ¼ fdξ, dη, dζg

E c ¼ 12ðL þ LT Þ

dX ¼ F dξ

where components of E c are as follows:

ð3Þ

and F −1 ¼ ½a1ζ ; a2ζ ; a3ζ T

ð16Þ

The covariant strain tensor can be calculated using ð17Þ

with F the covariant basis, relative to the point q, it is given by

Eξξ ¼ aT1ζ  u,ξ ,

2Eξη ¼ aT1ζ  u,η þ aT2ζ  u,ξ

F ¼ ½a1ζ ; a2ζ ; a3ζ 

Eηη ¼ aT2ζ Eζζ ¼ aT3ζ

 u,η ,

2Eξζ ¼ aT1ζ  u,ζ þ aT3ζ  u,ξ

 u,ζ ,

2Eηζ ¼ aT2ζ  u,ζ þ aT3ζ  u,η

ð4Þ

where a1ζ and a2ζ (Fig. 2) are the covariant or natural basis vectors defined by [6,15]

Therefore, we obtain a simple relationship between the curvilinear and covariant strain tensors

a1ζ ¼ X p,ξ þ 12 ζV ,ξ a2ζ ¼ X p,η þ 12 ζV ,η a3ζ ¼ V 1 2

ð5Þ

The element volume dV of the hexahedron can be obtained using dV ¼ J dξ dη dζ

ð6Þ

with J ¼ detðFÞ ¼ jða1ζ  a2ζ Þa3ζ j

ð7Þ

An orthonormal basis can be constructed for each point q and defined by Q ¼ ½t 1ζ ; t 2ζ ; nζ 

ð8Þ

where t 1ζ and t 2ζ are unit vectors in the plane Aζ for ζ ¼ Cst. The normal unit vector is defined by nζ ¼

a1ζ  a2ζ ja1ζ  a2ζ j

ð9Þ

Several techniques may be used to define the tangent unit vectors t 1ζ , t 2ζ , following Batoz and Dhatt [6], one can define 2 3 1 2 1 6 β þ 1 þ βnζy −1 þ βnζx nζy nζx 7 6 7 6 7 1 1 2 ð10Þ Q ¼6 7 nζx nζy β þ nζx nζy 7 6− 4 1þβ 5 1þβ −nζx −nζy nζz with β ¼ nζ  k (see Fig. 2). Finally the relationship between the curvilinear basis and the global cartesian framework can be expressed by dX ¼ Q ds

ð11Þ

E ¼ C T Ec C

ð18Þ

Now, by using the engineering notation, relation of Eq. (18) may be rewritten in a simple form E ¼ C~ E c with 2

C 211 6 6 C 212 6 6 2C C 6 ~ C ¼ 6 112 12 6 C 13 6 6 2C C 4 11 13 2C 12 C 13

ð19Þ

C 221

C 11 C 21

0

0

C 222 2C 21 C 22

C 12 C 22 C 12 C 21 þ C 11 C 22

0 0

0 0

C 223

C 13 C 23

C 233

C 13 C 33

2C 21 C 23

C 23 C 11 þ C 21 C 13

0

C 11 C 33

2C 22 C 23

C 23 C 12 þ C 22 C 13

0

C 12 C 33

0

3

7 7 7 7 7 7 C 23 C 33 7 7 C 21 C 33 7 5 C 22 C 33 0 0

T

and E c ¼ fEξξ Eηη 2Eξη Eζζ′ 2Eξζ′ 2Eηζ′ g , E ¼ fExx Eyy 2Exy Ezz′ 2Exz′ 2Eyz′ gT . 2.2. Material law in convected basis The constitutive relation of laminated composites [39] can be described using an orthotropic material law. For that purpose, we express the components of the tangent elastic moduli tensor relative to the fiber reference axis fm1 ; m2 ; m3 g of lamina [23] 2 3 H 1111 H 1112 0 H 1113 0 0 6 1122 7 6H H 2222 0 H 2223 0 0 7 6 7 6 0 1212 0 0 0 7 0 H 6 7 H ¼ 6 1133 ð20Þ 7 2233 3333 6H 7 H 0 H 0 0 6 7 6 0 7 1313 0 5 0 0 0 H 4 0 0 0 0 0 H 2323 Each layer is constituted of an orthotropic elastic material directions fm1 ; m2 ; m3 g (see Fig. 3). The non-zero coefficients of

4

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

Fig. 4. Orthonormal basis for composite shell.

Fig. 3. Fiber-reinforced lamina orientation axes.

the behavior matrix are given by the following expressions: H 1111 ¼ E1 ð1−ν23 ν32 Þ=K, H

2222

H 1122 ¼ E1 ðν21 þ ν23 Þ=K H 1113 ¼ E1 ðν31 þ ν21 ν32 Þ=K

¼ E2 ð1−ν13 ν31 Þ=K,

H 2223 ¼ E2 ðν32 þ ν12 ν31 Þ=K, H 1212 ¼ G12 , H

2323

¼ G23 ,

H 3333 ¼ E3 ð1−ν12 ν21 Þ=K

for i,j ¼ 1,2,3 and i≠j

H 1313 ¼ G13 νij Ej ¼ νji Ei

with K ¼ 1−ν12 ν21 −ν13 ν31 −ν23 ν32 −ν12 ν23 ν31 −ν13 ν21 ν32 , E1, E2 and E3 are the elastic moduli in the principal material directions fm1 ; m2 ; m3 g, νij and Gij are respectively the Poisson's ratio and the shear modulus. Since matrix is associated with the principal material directions, we need to transform it from the lamina coordinate axes fm1 ; m2 ; m3 g to the global cartesian coordinate axes fX; Y; Zg. With θ being the fiber orientation angle relative to the global cartesian system (Fig. 3), the relationship between the lamina coordinate system and the global Cartesian system is given by m1 ¼ cos θX þ sin θY m2 ¼ −sin θX þ cos θY m3 ¼ Z

ð21Þ

The final constitutive tensor H~ can be expressed in the convective coordinates as H~ ¼ RT HR with 2

ðr 1 Þ2 6 1 2 6 ðr 1 Þ 6 2 6 1 1 6 2r 1 r 2 R¼6 6 ðr 1 Þ2 6 3 6 1 1 6 2r r 4 2 3 2r 11 r 13

ð22Þ

ðr 21 Þ2

r 11 r 21

ðr 31 Þ2

r 21 r 31

ðr 22 Þ2 2r 21 r 22

r 12 r 22 1 2 r 1 r 2 þ r 21 r 12

ðr 32 Þ2 2r 31 r 32

r 22 r 32 2 3 r 1 r 2 þ r 22 r 31

ðr 23 Þ2

r 13 r 23 1 2 r 2 r 3 þ r 13 r 22 r 11 r 23 þ r 13 r 21

ðr 33 Þ2

r 23 r 33 2 3 r 2 r 3 þ r 23 r 32 r 21 r 33 þ r 23 r 31

2r 22 r 23 2r 21 r 23

2r 32 r 33 2r 31 r 33

r 11 r 31

3

7 7 r 22 r 32 7 1 3 1 37 r 1 r 2 þr 2 r 1 7 7 7 r 13 r 33 7 7 1 3 1 37 r 2 r 3 þr 3 r 2 5 r 11 r 33 þr 13 r 31 ð23Þ

r ji

Q jζ mi .

and ¼ The stresses can be evaluated in the curvilinear coordinate system (Fig. 4) as S ¼ H~ Z where the stress tensor S ¼ fSxx Syy Sxy Szz′ Sxz′ Syz′ g

ð24Þ T

2.3. Variational formulation using solid-shell element The principal of virtual work can be stated as Wðu,δuÞ ¼ W int −W ext ¼ 0

ð25Þ

∀δu≠0 and u ¼ 0, u ¼ u on Su. With Su is the shell contour where displacements are imposed. The Principal of Virtual Work can be expressed in the curvilinear

Fig. 5. Composites modeling using solid-shell elements.

coordinate system as Z Z Z W ¼ W int −W ext ¼ δE T S dV− δuT f v dV− δuT f s dS V

V

ð26Þ

Sf

where f V and f s are volume and surface traction forces respectively. In order to deal with the several locking phenomena separately, we need to split the expression of virtual internal work by separating the membrane/bending, thickness and transverse shear effects [6,15]. W int ¼ W mb þ W tr þ W sh

ð27Þ

with

Z W mb ¼ δE Ts H~ 1 E s dV Z V W tr ¼ δE Tz H~ 2 E z dV ZV W sh ¼ δGTs H~ 3 Gs dV

ð28Þ

V

and E s ¼ fExx Eyy 2Exy gT , E z ¼ fExx Eyy Ez′z′ gT , Gs ¼ f2Exz′ 2Eyz′ gT . H~ 1 , H~ 2 and H~ 3 are sub-matrices extracted from the global material matrix H~ according to the different strain components. In the context of modeling multilayered composite structures using solid-shell elements, there exist two possibilities regarding the numerical implementation:  Case of one element per layer: This is the easiest way for modeling the entire thickness of the structure using several elements (1 element per layer) as indicated in Fig. 5. In this case the numerical procedure of integration is straightforward and does not require any efforts in the implementation compared to the case of isotropic material modeling. The user has to provide the following basic properties: 1. Declaration of n groups of different elements in the FE mesh, these groups correspond to the n different material layers constituting the laminate. 2. Provide the physical characteristics of each layer.  Case of several layers per element: This second technique consists in stacking the different material layers within the same element (see Fig. 5). Each FE is a stack of several layers,

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

therefore stress calculation and numerical integration of the stiffness matrix is carried out using the “single-layer approach” identical to the one commonly used for the integration of plasticity, by using either one or several Lobatto integration points through the thickness for each layer. In the present investigation, this second method has been implemented because it is more general and more convenient since it requires only one element in the thickness of the structure (gain of CPU time in the case of industrial applications). On the other side the drawback the chosen element kinematic will not be able to capture delamination effects. In this case, the numerical implementation of the FE model requires some modifications of the stresses calculation and the integration of the stiffness matrix. For each element the user has to provide the following properties: 1. Declaration of n groups of materials, these n groups correspond to the different layers constituting the multilayered structure. 2. Provide the physical characteristics for each material. i 3. Provide for each layer i, the thickness hi (with ∑ni¼ 1 h ¼ h total thickness of the structure). For instance, the numerical implementation for the in-plane membrane/flexion effect is done in the following way: ! Z Z Z Z W mb ¼

z ¼ þ h=2

A

z ¼ −h=2

δE Ts S s

z ¼ þ hi =2

nl

dz dA ¼ ∑

i¼1

Ai

z ¼ −hi =2

δE Ts S s

dz dA ð29Þ

where nl is the number of layers, hi represents the layer thickness and h is the total thickness of the structure. Eq. (29) implies that for each layer of material, all operators are reported to the mid-plane of every layer i then the numerical integration is performed using Lobatto integration scheme with N points through the thickness direction. In the present model, 2 Lobatto integration points are used in the thickness direction for each layer, since all applications involve only elastic material behavior, for the case of material nonlinearities, at least 5 integration points will be necessary to achieve a correct stress integration. A solid-shell element formulated using Eqs. (26)–(28) with standard integration based on 2  2  N Gauss schema will fail because of numerous locking phenomena. 2.4. Remedies for shear locking An effective method of resolving shear locking is the Assumed natural Strain method in which the natural transverse shear strains are sampled and then interpolated at some discrete element points. The transverse shear strains 2Eξζ and 2Eηζ are calculated according to the average surface plan (ζ ¼ 0), assuming that they vary linearly (Fig. 6), and are function of γ ξ and γ η at the mid-side points: 2EANS ξζ ¼

1−η A1 1 þη A2 γ þ γ 2 ξ 2 ξ

ð30Þ

2EANS ηζ ¼

1−ξ B1 1 þ ξ B2 γ þ γ 2 η 2 η

ð31Þ

5

Fig. 6. Assumed natural transverse shear strains.

Fig. 7. ANS interpolation for transverse normal strain.

element mid-surface (Fig. 7), namely 4 1 ANS E~ z ¼ ∑ ð1 þ ξi ξÞð1 þ ηi ηÞEζ ðξi ,ηi Þ i¼14

ð32Þ

Poisson's ratio coupling requires the thickness strain to be a linear function of ζ. Because our solid-shell element has only two layers, as consequence the thickness strain does not vary with ζ thus the element fail in reproducing the plane-stress condition. In order to obtain a linear distribution of the normal strain in thickness direction, we enhance the thickness strain field by adding an internal degree of freedom as E EAS ¼ E~ z z

ANS

þ αζeζ

ð33Þ

where α is the 7th independent internal parameter. Eq. (33) is known as the Enhanced Assumed Strains (EAS), in which the enrichment variable α will be eliminated by static condensation technique within the element level. By substituting Eqs. (30–31) and (33) into Eq. (28), the virtual internal work takes the final expression Z Z Z ζi þ 1 nl ~ EAS W int ¼ ∑ ðδETs H~ 1 E s þ δE EAS z T H 2 Ez ξ

i¼1

ζi

η

þ δGANS T H~ 3 GANS ÞJ dζ dη dξ s s

ð34Þ

where nl is the number of element layers and ζ i is the transverse coordinate position of the ith layer along the thickness. After applying the FE discretization with the use of Eqs. (30), (31) and (33) into the previous equation, one can obtain T

ANS T W int ¼ δuTn ðK mb þ K EAS tr þ K sh Þun þ δαkα α þ δαk αu  un þ δun kαu α

ð35Þ

2.5. Remedies for thickness and trapezoidal locking Similar to shear locking, trapezoidal locking occurs when lower order elements such as 8-node hexahedral elements are used to model curved shells so that their cross-sections assume the trapezoidal shape these excessive number of sampled thickness strains can be reduced by using a bilinear interpolation of the transverse normal strains sampled at the four corners of the

The expressions of kαu and kα for the composite multilayered element are given by Z Z Z ζi þ 1 nl 33 kαu ¼ ∑ ζ H~ 2 BTz eζ J dζ dξ dη ð36Þ ξ

i¼1

nl

kα ¼ ∑

i¼1

ζi

η

Z Z Z ξ

η

ζi þ 1 ζi

33 ζ 2 H~ 2 eζ J dζ dξ dη

ð37Þ

Finally we carry out a static condensation, in order to eliminate the 7th unknown parameter α in the previous expression which T leads to α ¼ ð−1=kα Þkαu  un . Then the final stiffness matrix takes

6

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

given by (Fig. 9)

the following expression: 1 T ANS kαu ⊗kαu K ¼ K mb þK EAS tr þ K sh − kα

ð38Þ

The resulting 8-node hexahedral solid-shell FE model with 24 degrees of freedom is called SCH8γ7.

3. Optimization using RSM based on Diffuse Approximation The optimization problem can be stated as [30] Minimize f ðxÞ,

x∈Rn

ð39Þ

subject to a set of m þ2n constraints: g i ðxÞ≤0, j ¼ 1,…,m Li ≤xi ≤U i , i ¼ 1,…,n

ð40Þ

where f is the objective function (OF), xi are the design variables, gj is the jth constraint. The region of interest is defined by Li and Ui which are respectively the lower and upper bounds on the design variables. The RSM approach consists in solving a problem where the OF is replaced by its approximation f~ . This new problem may be written as Minimize f~ ðxÞ, x∈Rn Subject to g~ i ðxÞ≤0, j ¼ 1,…,m Li ≤xi ≤U i , i ¼ 1,…,n

wðxi ,xÞ ¼ wðrÞ ¼

expðβ2 r 2 Þ−expðβ2 Þ 1−expðβ2 Þ

ð42Þ

where r ¼ ∥xi −x∥=di , with di being the radius of a spherical support. As we can observe in Fig. 9, the parameter β affects the shape of the weight function. Generally, small values of the parameter β lead to a smooth and diffuse approximation but with the cost of a slow convergence process, and for larger values of β the optimization process converges more rapidly but increases the risk of divergence. Through the present investigation, it has been found by the authors, that β≥2 (Fig. 9) is a good compromise since it insures a good convergence rate while maintaining a certain amount of diffusive character on the approximation. The approximation is local, which means that only the points closest to the current optimum are taken into account (Fig. 10). The approximation coefficients are continuous when panning and/or zooming of the region of interest is performed. Given the function values for a set of experimental points xi distributed according to a chosen Design of Experiment, the function f~ can be defined in terms of basis functions p and some adjusting coefficients a as f~ ðxÞ ¼ pT ðxÞ  aðxÞ

ð43Þ

ð41Þ

A common choice for the basis functions p are linear and quadratic monomials

The approximation Eq. (41) is based on a set of numerical experiments with the function f. Generally, the approximate functions encountered in RSM rely on second-order models, over a given region of interest, with the constant regression coefficients fitted by means of least squares. The idea in this work is to apply minimization algorithm progression by building new response surfaces centered each successive solution (Fig. 8). During the progression of the process, the region of interest moves and new numerical experiments are performed at each iteration, which is known as the Moving Least Square (MLS) approximation. In the present investigation, we explore the application of DA regression for building the response surface during the iterative process. In the DA Method used in this work [26,27], the approximating functions are polynomials fitted to the nodal values of each local domain by a weighted least squares approximation. Belytschko et al. [7] developed an Element Free Galerkin (EFG) method which is an alternative implementation using Moving Least Square (MLS) approximation [24]. In this paper we adopt the moving least squares (MLS) interpolation. In the DAM the idea is to replace the OF f(x) computed using the FE method, for a local moving square fitting. The resulting function f~ is more regular that the function of the FE method, since the discontinuous coefficients are placed, by continuous functions of weight, which gives a continuity C m ðm≥1Þ. The approximate function becomes smooth by using continuous weighting functions. Different weight functions have been proposed in the literature, they differ in both the shape of the domain of influence, and in functional form. The truncated Gaussian spherical weight function used in this work is

 T x2 x2 pðxÞ ¼ 1 x1 x2 x3 … xn x1 x2 x1 x3 … xi xi þ 1 … 1 … n 2 2

Fig. 8. Surface response optimization using the MLS approximation.

Fig. 10. Local moving response surface.

The size N of the vector p depends on the number of variables n and on the polynomial degree of the approximation, for instance if

Fig. 9. Gaussian weight function: influence of parameter β.

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

quadratic approximation is used, then N ¼ ðn þ 1Þðn þ 2Þ=2. The generalized coefficients a are determined by performing a weighted least squares fit for the local approximation, which is obtained by minimizing the error JðaÞ between the experimental and approximated values of the OF

7

where M is the number of performed experiments and xi are the experimental designs. The weights wi insure the continuity and the locality of the approximation and are defined wi > 0, decreasing within a fixed region around the point i called domain of influence of xi and vanish outside. The weight functions play a crucial role by influencing the way that the coefficients a depend on the location of the design point x. The minimization of JðaÞ leads to the coefficients a and consequently the function to the expression of the function f~ . Thus, to obtain nonlinear parameters vector a, we use

The structure is a two-layer, composite, orthotropic, square plate that is simply supported on its bottom edges. The layers are oriented at 7 451 with respect to the plate edges. Fig. 11 shows the plate dimensions and the loading pressure which is applied on the top surface of the plate. Each layer has the material properties given in Table 1. Table 2 summarizes the results of deflection by comparing the normal displacement at the center of the plate and the in-plane displacement at (x ¼ 0,y ¼ b=2) to the analytical solution and to Á those obtained using ABAQUS by means of the SC8R solid-shell FE model. We can observe from Table 2; that the present model gave more precise results than those obtained with the solid-shell Á model of ABAQUS with the same number of elements. For instance, with only 10 elements, we obtain a maximum deflection of 22.461 mm which represents 3.4% of error when compared to the analytical solution 23.250 mm, while the solution obtained Á using the SC8R solid-shell element of ABAQUS is 28.543 mm which presents and error of 22.7%.

∂J ¼ AðxÞa−BðxÞZ ¼ 0 ∂a

4.2. Laminated strip under 3-point bending

M

JðaÞ ¼ ∑ wðxi ,xÞðpT ðxi Þ  aðxÞ−f ðxi ÞÞ2

ð44Þ

i¼1

ð45Þ

with ANN ðxÞ ¼ P T WðxÞP

ð46Þ

BNM ðxÞ ¼ P T WðxÞ

ð47Þ

and

2

6 6 W MM ¼ 6 6 4

wðx1 Þ

0 wðx2 Þ ⋱ wðxM Þ

0 2 6 P MN ¼ 4

pT ðx1 Þ ⋮ pT ðxM Þ

3 7 5

and

One-quarter of the laminated strip is modeled. The same problem is analyzed using different meshes of the developed 3D solid-shell FE model. Various options are used to model the

3 7 7 7 7 5

8 9 > < f ðx1 Þ > = ⋮ Z MN ¼ > : f ðx Þ > ; M

ð48Þ

ð49Þ

Once the objective function is estimated, then an optimization algorithm can be used to minimize the resulting OF. In this work we have adopted an algorithm based on a Sequential Quadratic Programming technique [31] and a second algorithm based on the Simplex method. Both algorithms are robust enough and suitable to deal with constrained nonlinear optimization problems. But any Minimization algorithm may be used at this stage. We have to notice, that as the minimization procedure is an iterative process, therefore at least few iterations (generally less than 10) are often needed for the SQP algorithm in order to find the optimal solution of the metamodel obtained using the RSM based on DA.

Fig. 11. Geometry and loading of the anisotropic layered plate.

Table 1 Material properties of the anisotropic layered plate. Material property

Value

Units

E11 E22 ¼ E33 G12 ¼ G13 ¼ G23 ν13 ¼ ν23 ν12

27,6000 6900 2000 0 0.25

MPa MPa MPa

4. Numerical applications In order to model defined with the aim phenomena, as

evaluate the effectiveness of the solid-shell FE previously, several benchmarks are carried out of studying its behavior regarding the locking well as the global rate of convergence.

4.1. Analysis of an anisotropic layered plate The problem considered is the analysis of a flat plate made from two layers oriented at 7451, subjected to a uniform pressure loading. The example verifies simple laminated composite plate analysis. Our numerical results are compared with the analytical solution given in Spilker et al. [36] and a second verification is Á done by a comparison to the results obtained using ABAQUS SC8R solid-shell FE model. In this benchmark, the cross-section is not balanced, so the response includes membrane-bending coupling.

Table 2 Plate deflection convergence. No. of elements

uA (mm) Present model

uA (mm)

uC (mm) Present model

uC (mm)

ABAQUS

2 4 6 8 10 20 30 40 50 100

0.335 0.360 0.365 0.368 0.370 0.371 0.371 0.371 0.371 0.371

0.423 0.414 0.416 0.417 0.417 0.418 0.418 0.419 0.419 0.419

17.960 20.806 21.766 22.204 22.461 22.956 23.114 23.191 23.236 23.241

29.435 28.649 28.584 28.555 28.543 28.460 28.475 28.480 28.485 28.495

Analytical

0.376

0.376

23.250

23.250

Á

Á

ABAQUS

8

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

laminated strip through the thickness. The models using the developed SCH8γ7 solid-shell employ a single element in the thickness direction using a composite section based on 7 singlelayers through the thickness. The structure is a 7-layer, composite, orthotropic, rectangular strip that is simply supported on its supports on A and B and subjected to a line load of 10 N/mm at C onto the top surface (Fig. 12). The layers are oriented at 01/901 consecutively with respect to the strip edges. Each layer has the material properties given by Table 3 below. Our numerical results are compared with the reference solution given by NAFEMS [28] and a second verification is done using a comparison with the results obtained using the SC8R solid-shell Á element of ABAQUS software. Table 4 summarizes the results by comparing the normal displacement at the center of the strip (point E) to the reference Á solution and to those obtained using ABAQUS by means of the SC8R solid-shell FE model. The displacement field obtained using the present SCH8γ7 element is of a very good accuracy. We can see from Table 4, that even with only 1 element per edge AE, the obtained deflection (−0.807 mm) presents an error less than 23.9% compared to NAFEMS reference solution, while the obtained Á solution using ABAQUS solid-shell element, gave an error of 67.5%. The error vanishes as the number of element increases (only through the edge, while keeping only 1 element in the

thickness), for only 9 elements per edge AE we obtained a deflection (−1.059 mm) with an error less than 0.1%. 4.3. Multilayered composite plate with ply drop-offs This example demonstrates the applicability of the present solid-shell formulation to analyze a composite structure with ply drop-offs; it consists on a composite multilayered plate (with 6 layers) with ply drop-offs (Fig. 13). In this example, each layer is made of unidirectional fiberreinforced material [40], with the fiber directions aligned at þ 45=−45= þ 45=−45= þ 45=−451 with respect to the length direction. The structure, with length L¼ 12 m, width b ¼6 m and thickness h¼ 0.1 m is clamped on its thicker side and the free thinner end is subjected to a transverse normal load distribution uniformly along the free edge equivalent to a concentrated force of F¼600 N. The location of the ply drop-offs are at x ¼4 m and x¼ 8 m with the top two layers removed after each drop-off. The layer material properties are given in Table 5. Table 6 shows the computed solution obtained using the developed SCH8γ7 solid-shell model with 12  6 elements (12 along the length and 6 elements along the width). One can observe that the present model can predict accurately the global deflection even for a very thin structure (L=h ¼ 12=0:004 ¼ 3000). Comparison of the maximal deflection obtained with the present model with the Á one obtained using ABAQUS SC8R solid-shell element with the same mesh (Fig. 13), gave a small error of only 2.4%, which shows that the present model is free of shear locking. 4.4. Optimization of a warped thick cylinder under pressure In this application a warped thick cylinder submitted to an internal pressure is studied [28]. The continuum shell model employs one single element with a composite section based on two single-layers through the thickness. The structure is

Fig. 12. Geometry and loading of the laminated strip.

Table 3 Material properties of the laminated strip. Material property

Value

Units

E11 E22 ¼ E33 G12 G13 ¼ G23 ν13 ¼ ν23 ν12

100,000 5000 3000 2000 0.3 0.4

MPa MPa MPa MPa

Fig. 13. Multilayered composite plate with ply dropoffs: undeformed mesh.

Table 5 Material properties of the multilayered composite plate.

Table 4 Convergence of the laminated strip deflection. No. of elements (edge AE)

uE (mm) Present model

ABAQUS

1 3 6 9 12 15 18 30 60

−0.807 −1.034 −1.055 −1.059 −1.061 −1.061 −1.062 −1.062 −1.062

−1.776 −1.062 −1.052 −1.056 −1.057 −1.057 −1.058 −1.058 −1.058

NAFEMS

−1.060

−1.060

uE (mm) Á

Material property

Value

Units

E11 E22 ¼ E33 G12 ¼ G13 ¼ G23 ν12 ¼ ν13 ¼ ν23

25,000 1000 1000 0.2

MPa MPa MPa

Table 6 Deflection of the multilayered composite plate. Thickness

Deflection (m) ABAQUS

Deflection (m) Present model

0.0039 2.6933 40.9174

0.0029 2.5685 39.9480

Á

0.1 0.01 0.004

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

9

Fig. 14. Geometry and FE mesh of a warped thick cylinder.

composed of an inner isotropic cylinder E ¼210 GPa, ν ¼ 0:3 and an outer orthotropic circumferentially wound cylinder having the following material parameters, expressed in a cylindrical coordinate system: a circumferential modulus E11 ¼ 130 GPa, a longitudinal and a transverse moduli E22 ¼ E33 ¼ 5 GPa (as indicated in Fig. 14), ν12 ¼ ν13 ¼ 0:25, ν23 ¼ 0, G12 ¼ G13 ¼ 10 GPa, and G23 ¼ 5 GPa. The boundary conditions correspond to a displacement dz ¼ 0 at z¼0 and the cylinder is subjected to an internal pressure of 200 MPa (Fig. 14). At first and prior to the optimization of the material parameters, the cylinder is modeled using the initial material data and a mesh of 40 elements through its circumference, 10 elements through the length and only 1 element in thickness direction. The circumferential stress at z ¼0 is measured for 2 radius values (R ¼23 mm, R¼27 mm) and compared to the NAFEMS reference Á solution [28] and to the one obtained using the ABAQUS SC8R solid-shell FE model. For the inner flange we obtain S11 ¼ 1534:8 MPa, the SC8R Á ABAQUS gave S11 ¼ 1477 MPa, while the reference solution was 1565 MPa. For the outer flange we obtain S11 ¼ 892:8 MPa, the Á ABAQUS solution was S11 ¼ 900 MPa, while the reference solution was 875 MPa. These results confirm that our solid-shell model and the used mesh are good enough to carry out the optimization process. The optimization problem consists in finding optimal fiber orientation angle θ of the outer orthotropic circumferentially wound layer and the Young's modulus E of the inner isotropic cylinder while keeping constant the cylinder expansion for z ¼0 mm at a value of uR ¼ 0:5 mm. The objective function is based on the general Hill criterion [42]. nel

Jðθ,EÞ ¼ J ¼ ∑ J e

Fig. 15. Design of experiments: (a) central composite design and (b) Box-Behnken design.

Table 7 DOE plan using central composite design. Run

Factor x1

Factor x1

Response J (MPa2)

Disp. uR (mm)

1 2 3 4 5 6 7 8 9

−1 1 −1 1 −1 1 0 0 0

−1 −1 1 1 0 0 −1 1 0

268316.8 294069.1 282546.0 293785.5 278756.9 294282.9 280609.1 286957.9 285223.1

0.355 0.833 0.150 0.196 0.210 0.316 0.584 0.134 0.211

ð50Þ

e¼1

where e is the element index, nel is the total number of elements on the whole structure, the element OF is given by J e ¼ FðS22 −S33 Þ2 þ GðS33 −S11 Þ2 þ HðS11 −S22 Þ2 þ 2NS212 þ2LS223 þ2MS231 ð51Þ The design variables (x1 ¼ θ,x2 ¼ E) are constrained between 01≤θ≤901 and 60 GPa≤E≤300 GPa. Before performing the optimization procedure, a DOE based on composite design (Fig. 15a) built of 9 functions evaluations for the full cylinder, is carried out. In this case we used 3 groups of design points (4 two-level factorial design points, 4 axial points and a 1 center point). Then design variables are coded as (−1, þ 1) in order to facilitate the data treatment. The DOE plan with the OF values is given in Table 7. Fig. 16 shows the global quadratic response surface model based on DA which has been calculated and given explicitly by J~ ðθ,EÞ ¼ 291900 þ4107:30θ þ 1842:52E−7520:69θE −10230:38θ2 −261:44E2 The response approximation function for the constraint (uR ¼0.5) on the radial displacement of the cylinder at z ¼0 is also

Fig. 16. Approximated response surface.

10

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

carried out using (DA) and given explicitly by (see Fig. 17). ~ uðθ,EÞ ¼ 0:49 þ 0:026θ þ 0:071E þ 0:13θE−0:13θ2 −0:098E2

~ The minimization of J~ ðθ,EÞ under constraint uðθ,EÞ ¼ 0:5 (see Fig. 18) has been done using the SQP algorithm based on the work of Powell [31], the optimal solution was obtained in 5 iterations leading to the optimal solution in coded form (0.830, 0.961) which corresponds to ðθn ,En Þ ¼ ð82:341,295:40 GPaÞ. Fig. 19 shows the hoop stress distribution on the cylinder before and after optimization. We can observe just a few amount of stress reduction is obtained after optimization, this is due to the ~ presence of the constraint uðθ,EÞ ¼ 0:5 which cannot lead to the minimum of the unconstrained objective function J~ ðθ,EÞ.

4.5. Optimization of a multilayered strip panel

Fig. 17. Response surface of constraint function.

Fig. 18. Optimal solution in contour plot.

The structure is a thick composite multilayered shell (Fig. 20), composed by two layers having the same thickness. The lower layer is made of an aluminum alloy 6066-T651 (E¼67 GPa, ν ¼ 0:34) and the upper layer is an orthotropic lightweight based aramid fibers with strong characteristics (E11 ¼ 130 GPa, E22 ¼ E33 ¼ 5 GPa, ν12 ¼ ν13 ¼ 0:25, ν23 ¼ 0, G23 ¼ 5 GPa and G12 ¼ G13 ¼ 10 GPa). The structure is submitted to a transverse load in point A of 4.17 N/mm, and clamped at its end BC. The structure was meshed using 1596 3D-shell elements corresponding to 3588 nodes (Fig. 20). Prior to performing the material optimization, a FE model validation is carried out using initial material parameters. A numerical comparison of our results to those obtained with Á ABAQUS software by means of the SC8R solid-shell model is performed. Fig. 21 shows the deflection of the strip panel under the maximal load. Comparison of the deflection value obtained using the present SCH8γ7 solid-shell model gave 13.747 mm which Á is in good agreement with 13.676 mm obtained using ABAQUS . The optimization problem consists in finding optimal fiber orientation angle θ of the upper orthotropic lightweight based aramid layer and the Young's modulus E of the lower aluminum alloy layer, while the maximal deflection under load is maintained equal to 15 mm. The objective function is based on the general Hill criterion [1] and defined by Eq. (51).

Fig. 20. Initial undeformed mesh and loading conditions of a multilayered strip panel.

Fig. 19. Hoop stress distribution.

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

11

Fig. 21. Comparison of the global deflection of a multilayered strip panel (SC8R vs. SCH8γ7). (a) ABQUS (SC8R). (b) Present model ðSCH8γ7Þ.

Fig. 22. Approximated response surface.

The optimization problem is stated us nel

Minimize

Jðθ,EÞ ¼ ∑ J e ,

Fig. 23. Optimal solution in contour plot.

ðθ,EÞ∈Rn obtained for the optimal structure is about un ¼ 15:11 mm which respects the imposed constraint u≥15 (Fig. 23).

e¼1

with

uðθ,EÞ−15 mm≥0,

−22:51≤θ≤ þ 22:51, 60 GPa≤E≤79 GPa:

ð52Þ

The response surface model based on DA is built up using the same DOE based on central composite design using 9 experiments. Fig. 22 shows the approximated RSM calculated which is given by the quadratic metamodel J~ ðθ,EÞ ¼ 939800 þ58864:99θ þ 29763:57E þ 11549:15θE þ 146500θ2 þ 1431:34E2

The response approximation for the constraint (uðxÞ−15 ¼ 0) has been also carried out using DA by means of a quadratic metamodel, it is given explicitly by ~ uðθ,EÞ ¼ 13:08 þ3θ−1:22E−0:17θE þ4:58θ2 þ0:14E2 The optimal solution obtained after un-normalizing the DV corresponds to ðθn ,En Þ ¼ ð3:401,60:12 GPaÞ, the displacement

4.6. Optimization of a composite leaf spring for automobile suspension In this application, a composite leaf spring for automobile suspension is studied. The structure shown in Fig. 24 with a 10 mm constant thickness is made of a composite orthotropic material whose mechanical characteristics are given in Table 8. The optimization problem here, consists in finding the best orientation of the material fibers which allows reducing the maximal stresses over the structure under ultimate loading conditions. The FE mesh using 1220 solid-shell elements and the boundary conditions and loading are shown in Fig. 24. The structure supports a concentrated force at its free end (applied on the bottom shell surface) of 500 N, the other edge being clamped. Prior to the optimization process, the FE model has been validated using initial material parameters with an orientation angle θ ¼ 01 with reference to the global x-axis (Fig. 24). A first

12

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

Fig. 24. Geometry and boundary conditions of a composite leaf spring suspension.

Table 8 Material properties of a composite leaf spring suspension. Material property

Value

Units

E11 E22 ¼ E33 G12 ¼ G13 ¼ G23 ν12 ¼ ν13 ¼ ν23

39,000 10,000 5000 0.3

MPa MPa MPa

Á

analysis is carried out using ABAQUS software and its SC8R solidshell element with the same mesh. Fig. 25 shows the stress distribution on the top flange of the structure in the 1st principal direction. Comparison of maximal values of S11 obtained using the present solid-shell model gave 185.79 MPa which is almost the same maximal value 185.41 MPa obtained with the SC8R element Á of ABAQUS . Unlike the stress S11, the distribution of the transverse stress S33 in the thickness direction, shows that a maximal value of Á 20 MPa is reached, while the ABAQUS solution indicates zero stress as indicated in Fig. 26 (This results from the fact that the Á SC8R solid-shell element of ABAQUS is based on plane stress condition). Analysis of other stress components shows small values compared to the principal stress S11 values and a comparÁ ison with ABAQUS results shows good agreement. This first study confirms that the present model is valid and therefore the optimization process can be performed. In order to minimize the number of design variables, a 1-dimensional parameterization is used through the curvilinear abscissa s where 5 DVs are defined as indicated in Fig. 27. The points represented by θ1 ,θ2 ,…,θ5 are the optimization poles of control, where the fibers orientation angle θi is defined. Inside each finite element, the fiber orientation angle θe is calculated using a linear interpolation between two successive poles where the element is located. For instance, if an element is located between two poles θi and θj , the material orientation angle inside the element is defined as

 sj −s s−si θðsÞ ¼ θi þ θ sj −si sj −si j Here again, the objective function is based on the Hill criterion [1] and defined by Eq. (51). Thus the optimization problem can be stated in terms of the DV vector θ ¼ fθ1 ; θ2 ; θ3 ; θ4 ; θ5 g as follows: nel

Minimize

JðθÞ ¼ ∑ J e , e¼1

with

−601≤θi ≤þ 601,

θ∈Rn i ¼ 1,5

ð53Þ

In order to reduce the number of evaluations of the OF, a DOE based on the Box-Behnken algorithm is used with three levels coded between {−1, 0, þ1} (Fig. 15b). The Box-Behnken DOE with 5 variables, allows a screening of only 46 experiments (or combination of variables) coded between −1 and 1. The quadratic response surface model based on DA is built up using the above mentioned DOE based on 46 experiments in the space of 5 coded DV. Fig. 28 shows the approximated RSM calculated and given by the following quadratic metamodel: J~ðθÞ ¼ 1874:04−1206:35θ1 −2210:64θ2 −1240:79θ3 −331:03θ4 −71:52θ5

þ 2146:90θ1 θ2 þ 0:42θ1 θ3 þ 1:50θ1 θ4 þ 0:30θ1 θ5 þ 1281:27θ2 θ3 þ 4:25θ2 θ4 þ 1:14θ2 θ5 þ 444:04θ3 θ4 þ 2:08θ3 θ5 þ64:00θ4 θ5 þ 1245:75θ21 þ 2475:70θ22 þ 1375:84θ23 þ269:35θ24 −31:33θ25 The minimization of the resulting metamodel was carried out with two different algorithms, namely the VF02AD algorithm using the SQP method [31] and a second optimizer based on the simplex algorithm. Different initial points were used to minimize the OF, and the best result obtained after minimization corresponds to θn ¼ f15:771; 15:451; 17:331; 15:191; 0:001g. After obtaining the optimal material fibers orientations, we carried out a last FE calculation using optimal solution θn found, by means of the developed SCH8γ7 solid-shell model. Review of the principal results is given in Table 9, where maximal stresses in all principal directions are compared to values obtained using the initial material data with θ ¼ 01. As we can observe, the optimal fibers orientation, allow a reduction of the maximal stress distribution in the structure, we found also that the obtained results if physically good since the orientation of the material fibers follows the general orientation of the structure geometry.

5. Conclusion In the present investigation, an efficient 8-node solid-shell element formulation for the analysis of multilayered composite shell is presented. While the ANS method has been used in order to remedy to shear locking, the enhancement of transverse normal strain is adopted, thus full 3D anisotropic constitutive model is incorporated without resorting to the plane-stress assumption. The present formulation can predict the through-thickness effects with a high degree of accuracy. In the second part of the paper we proposed a specific method based on RSM for the optimization of laminated structures, where

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

Fig. 25. Comparison of stress distribution on the top flange in the 1st principal direction. (a) ABQUS (SC8R). (b) Present model ðSCH8γ7Þ.

Fig. 26. Comparison of stress distribution on the top flange in the 3rd principal direction (thickness direction). (a) ABQUS (SC8R). (b) Present model ðSCH8γ7Þ.

Fig. 27. Parameterization model of design variables.

the design variables are the material fibers orientation. To solve the optimization problem, we proposed a technique which has the advantage of replacing a complex response model by an approximate model evaluated in a limited number of points obtained through an experimental design. In this study, we proposed a specific response surface method based on Diffuse Approximation involving pattern search optimization. The resulting response

Fig. 28. Approximated response surface in the space fθ1 ,θ2 g.

13

14

H. Naceur et al. / Finite Elements in Analysis and Design 70–71 (2013) 1–14

Table 9 Summary of principal results obtained before and after optimization. Maximal stress component

Initial solution using θ ¼ 01 (in [MPa])

Optimal solution using θn (in [MPa])

S11 S22 S12 S33

185.79 20.41 34.53 20.43

173.90 19.12 12.66 19.20

surface algorithm involves iterative improvement of the objective function employing locally supported nonlinear approximations. Numerous applications have been treated, they confirmed that the optimization method based on RSM and DA coupled to the developed SCH8γ7 FE solid-shell model is efficient and particularly suited for industrial problems in structural mechanics. References [1] S. Abid, J.L. Batoz, C. Knopf-Lenoir, P. Lardeur, H. Naceur, Thickness Optimization of Beams and Shells with Large Displacements, Integrated Design and Manufacturing in Mechanical Engineering, Editions Kluwer Academic Publishers173–182. [2] R.J. Alves de Sousa, J.W. Yoon, R.P.R. Cardoso, R.A. Fontes Valente, J.J. Gracio, On the use of a reduced enhanced solid–shell (RESS) element for sheet forming simulations, Int. J. Plasticity 23 (2007) 490–515. [3] P.M.A. Areias, J.M.A. César de Sá, C.A. Conceiç ao António, A.A. Fernandes, Analysis of 3D problems using a new enhanced strain hexahedral element, Int. J. Numer. Meth. Eng. 58 (2003) 1637–1682. [4] B. Bassa, F. Sabourin, M. Brunet, A new nine-node solid-shell finite element using complete 3D constitutive laws, Int. J. Numer. Meth. Eng. 92 (7) (2012) 589–636 16. [5] K.J. Bathe, E.N. Dvorkin, A formulation of general shell elements- the use of mixed interpolation of tensorial components, Int. J. Numer. Meth. Eng. 22 (1986) 697–722. [6] J.L. Batoz, G. Dhatt, Finite Element Modeling of Structures, 3rd Edition, vol. 3, Hermès 1992, Paris (in french). [7] T. Belytschko, Y.Y. Lu, L. Gu, Element free Galerkin method, Int. J. Numer. Meth. Eng. 37 (1994) 229–256. [8] P. Betsch, E. Stein, An assumed strain approach avoiding artificial thickness straining for a non-linear 4-node shell element, Commun. Numer. Meth. Eng. 11 (1995) 899–909. [9] M. Bischoff, E. Ramm, Shear deformable shell elements for large strains and rotations, Int. J. Numer. Meth. Eng. 40 (1997) 4427–4449. [10] N. Büchter, E. Ramm, Shell theory versus degeneration – a comparison in large rotation finite element analysis, Int. J. Numer. Meth. Eng. 37 (1992) 55–62. [11] R.P.R. Cardoso, J.W. Yoon, M. Mahardika, S. Choudry, R.J. Alves de Sousa, R.A. F. Valente, Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for one-point quadrature solid-shell elements, Int. J. Numer. Meth. Eng. 75 (2008) 156–187. [12] C. Cho, H.C. Park, S.W. Lee, Stability analysis using a geometrically nonlinear assumed strain solid shell element model, Finite Elem. Anal. Des. 29 (1998) 121–135. [13] D.R. Cox, N. Reid, The Theory of the Design of Experiments, Edition, Chapman Hall/CRC, 2000. [14] S. Doll, K. Schweizerhof, R. Hauptmann, C. Freischlager, On volumetric locking of low-order solid and solid-shell elements for finite elastoviscoplastic deformations and selective reduced integration, Eng. Comp. 17 (2000) 874–902. [15] E. Domissy, S. Bouabdallah, J.L. Batoz, Formulation and evaluation of a solid finite element model for the linear and nonlinear shell analysis, in: 12eme Congrès Français de Mécanique Strasbourg, France, September 1995 (in french). [16] D. Durban, D. Givoli, J.G. Simmonds, Advances in the Mechanics of Plates and Shells, Kluwer Academic Publishers, 2002.

[17] E.N. Dvorkin, K.J. Bathe, A continuum mechanics based four node shell element for general nonlinear analysis, Eng. Comput. (1984) 1–77. [18] M. Hannachi, H. Naceur, J.L. Batoz, Continuum based solid-shell element modeling for the optimization of composite multilayered structures, Int. Rev. Mech. Eng. 1 (4) (2007) 150–163. [19] M. Harnau, K. Schweizerhof, Artificial kinematics and simple stabilization of solid-shell elements occurring in highly constrained situations and applications in composite sheet forming simulation, Finite Elem. Anal. Des. 28 (3–4) (2006) 1097–1111. [20] R. Hauptmann, K. Schweizerhof, A systematic development of Solid-Shell element formulation for linear and non-linear analysis employing only displacement degrees of freedom, Int. J. Numer. Meth. Eng. 42 (1998) 49–69. [21] S. Klinkel, F. Grutmann, W. Wagner, A continuum based d-shell element for laminated structures, Comput. Struct. 71 (1999) 43–62. [22] R.L. Mason, R.F. Gunst, J.L. Hess, Statistical Design and Analysis of Experiments —with Applications to Engineering and Science, 2nd ed, John Wiley & Sons, 2003. [23] S.K. Mazumdar, Composites Manufacturing—Materials, Product, and Process Engineering, Edition, CRC Press LLC, Florida, 2002. [24] R.H. Meyers, D.C. Montgomery, Response Surface Methodology Process and Product Optimization Using Designed Experiments, 2nd Edition, Wiley, 2003. [25] R.A.S. Moreira, R.J. Alves de Sousa, R.A.F. Valente, A solid-shell layerwise finite element for non-linear geometric and material analysis, Compos. Struct. 92 (6) (2010) 1517–1523. [26] H. Naceur, A. Delameziere, J.L. Batoz, Y.Q. Guo, C. Knopf-Lenoir, Some improvements on the optimum process design in deep drawing using the Inverse Approach, J. Mater. Process. Technol. 146 (2) (2004) 250–262. [27] H. Naceur, Y.Q. Guo, S. Ben-Elechi, Response surface methodology for design of sheet forming parameters to control springback effects, Comput. Struct. 84 (2006) 1651–1663. [28] NAFEMS, National Agency for Finite Element Methods and Standards (NAFEMS), Test R0031/2 from NAFEMS publication R0031, Composites Benchmarks Issue 2, February 5, UK, 2001. [29] N.H. Nguyen, V.N. Pham, M. Hogge, J.P. Ponthot, An assumed natural strain technique for 2D solid-shell elements, in: M. Hogge et al. (Eds.), Fourth International Conference on Advanced Computational Methods in Engineering, University of Liège, Belgium, 26–28 May 2008. [30] J. Nocedal, S.J. Wright, in: P. Glynn, S.M. Robinson (Eds.), Numerical Optimization, Springer-Verlag, New York, 1999. [31] M.J.D. Powell, Extensions to subroutine VF02AD, in: R.F. Drenick, F. Kozin (Eds.), Systems Modeling and Optimization, Lecture Notes in Control and Information Sciences, vol. 38, Springer-Verlag, Berlin, 1982, pp. 529–538. [32] N.-D. Quy, A. Matzenmiller, A solid-shell element with enhanced assumed strains for higher order shear deformations in laminates, Tech. Mech. 28 (3–4) (2008) 334–355. [33] K. Rah, W. Van Paepegem, A.M. Habraken, J. Degrieck, A partial hybrid stress solid-shell element for the analysis of laminated composites, Comp. Meth. Appl. Mech. Eng. 200 (49–52) (2011) 3526–3539. [34] M. Schwarze, S. Reese, A reduced integration solid-shell finite element based on the EAS and the ANS concept—geometrically linear problems, Int. J. Numer. Meth. Eng. 80 (2009) 1322–1355. [35] S. Shiri, H. Naceur, J. Roelandt, Numerical modeling of stamping and crashworthiness of steel/polymer/steel structures using solid-shell element, in: International Conference on Computational Plasticity, COMPLAS X 2009, Barcelona, Spain, 2009. [36] R.L. Spilker, S. Verbiese, O. Orringer, S.E. French, E.A. Witmer, A. Harris, Use of the hybrid-stress finite-element model for the static and dynamic analysis of multilayer composite plates and shells, Report for the Army Materials and Mechanics Research, Center Watertown, MA, 1976. [37] K.Y. Sze, S. Yi, M.H. Tay, An explicit hybrid-stabilized eighteen-node solid element for thin shell analysis, Int. J. Numer. Meth. Eng. 40 (1997) 1839–1856. [38] K.Y. Sze, L.Q. Yao, A hybrid stress ANS solid–shell element and its generalization for smart structure modeling, Part I-solid–shell element formulation, Int. J. Numer. Meth. Eng. 48 (2000) 545–564. [39] J.R. Vinson, R.L. Sierakowski, The Behavior of Structures Composed of Composite Materials, 2nd Edition, Kluwer Academic Publishers, 2004. [40] L. Vu-Quoc, X.G. Tan, Optimal solid shells for non-linear analyses of multilayer composites. I. Statics, Comp. Meth. Appl. Mech. Eng. 192 (2003) 975–1016. [41] G. Wempner, D. Talaslidis, Mechanics of Solids and Shells—Theories and Approximations, Edition, CRC PRESS LLC, New York, 2003. [42] X. Xiao, C. Hsiung, Z. Zhao, Analysis and modeling of flexural deformation of laminated steel, Int. J. Mech. Sci. 50 (2008) 69–82.