- Email: [email protected]
An efficient computational strategy for composite laminates assemblies including variability
International Journal of Solids and Structures 50 (2013) 2749–2757
Contents lists available at SciVerse ScienceDirect
International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr
An efficient computational strategy for composite laminates assemblies including variability V. Roulet, P.-A. Boucard ⇑, L. Champaney LMT-Cachan (ENS Cachan/CNRS/UPMC/PRES UniverSud Paris), 61 av. du Président Wilson, F-94230 Cachan, France
a r t i c l e
i n f o
Article history: Received 6 September 2012 Received in revised form 29 March 2013 Available online 10 May 2013 Keywords: Multiparametric strategy LATIN solver Parallel computation Laminate composites Damage
a b s t r a c t The aim of this work is to present an efficient numerical strategy for studying the influence of the material parameters on problems involving 3D assemblies of composite parts with contact and friction. This approach is based on the multiscale LATIN method with domain decomposition and its ability to reuse the solution of one problem (for a particular set of parameters) to solve many similar problems, each associated with another set of parameters. The presentation of the proposed calculation method will be followed by two sample applications, one with variabilities in the cohesive interface laws and the other with variabilities in the damageable behavior laws, which characterize the composite material. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction In order to reduce design and sizing costs, manufacturing industries are seeking to limit the number of actual tests and replace them by numerical simulations. However, this approach may not be entirely compatible with the recent adoption by aeronautical industries of composite materials, whose behavior is complex and often imperfectly known. Taking these material uncertainties into account leads to the introduction of variabilities into the constitutive model. Starting with the hypothesis that the material properties of a composite material can be modeled within the framework of the probability theory, a very rich mathematical background is available to completely characterize the probabilistic behavior and the evolution of the structure under an external perturbation. When dealing with complex systems, approximate methods and numerical integration are necessary. If the uncertainties are small, the perturbation theory is a very valuable tool for analyzing their effect (see Kleiber and Hien, 1992). If they are large, then the perturbation theory is not applicable. One obvious way to solve Stochastic Partial Differential Equation (SPDEs) is the Monte Carlo method. This method is expensive, especially if higher-order accuracy for mean values, standard variations, queues of distributions, etc. is sought. Various numerical methods for solving SPDEs have been proposed in the literature in Ghanem and Kruger (1996).
⇑ Corresponding author. E-mail addresses: [email protected] (V. Roulet), boucard@lmt. ens-cachan.fr, [email protected] (P.-A. Boucard), champa [email protected] (L. Champaney). 0020-7683/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijsolstr.2013.04.028
Ghanem and Spanos (1991), and Ghanem (1999) propose a hybrid finite element/spectral approach called the Spectral Stochastic Finite Element Method (SSFEM). An extended review of these methods is available in Kaminski (2005). We propose here an alternative approach to reduce the computational costs. Furthermore, the simulation of assemblies of composite parts raises two difficulties: the interactions among the components, which are often complex and nonlinear with frictional contact problems, and the material’s behavior, especially its degradation. Contact problems are characterized by constraints such as nonpenetration conditions, and an active area of contact – that is, an area where contact effectively occurs – that is unknown a priori. For these reasons, these problems lead to stiff non-linear system of equations. Several approaches exist for solving static contact problems (Kikuchi and Oden, 1988; Zhong and Mackerle, 1992; Wriggers, 1995). In most of them, the numerical methods that are employed for enforcing the contact constraints can be grouped into Lagrange multiplier and penalty methods (Wriggers et al., 1985). The penalty methods, used in Kikuchi (1982) and Armero and Petcz (1998), are closely related to the regularization of the contact constraints. They are usually formulated in terms of the displacement variables and therefore are primal methods. They allow treating contact as a material behaviour, as exemplified by the method of joint finite elements in Alart and Curnier (1991). Penalty methods can experience various numerical difficulties, especially ill-conditioning, when a too large or too small penalty parameter is introduced. Lagrange multiplier methods are dual methods where the multipliers, which represent the contact reaction forces, are introduced in order to enforce exactly the non-penetration conditions. Augmented Lagrange multiplier methods, used in Arora et
2750
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
al. (1991), Klarbring (1992), Radi et al. (1998), Chabrand et al. (1998), result in mixed formulations involving both displacement and force unknowns. The numerical solution schemes underlying both the Lagrange multiplier and augmented Lagrange multiplier methods are often related to the Uzawa algorithm (see Simo and Laursen, 1992; Champaney et al., 1999). Consequently, it is necessary to develop an efficient numerical strategy dedicated to the resolution of nonlinear problems involving very large numbers of degrees of freedom. Variabilities may be introduced either at the level of the interactions among parts (friction coefficients, preloads, etc.) or in the constitutive material model itself. Thus, the solution of the problem is no longer the result of a deterministic calculation, but must be viewed as a response function of the input parameters. In the rest of this work, we will assume that the variation range of each parameter is divided into a finite number of values. Thus, each set of the parameters defines an associated problem. Then, a parametric study consists in solving as many problems as there are parameter sets. Therefore, one must be able to solve a very large number of ‘‘similar’’ problems. In order to do that, we will use a multiparametric strategy, which was already presented in Boucard and Champaney (2003), Champaney et al. (2008) and Champaney (2011) for assemblies with contact and friction, but only in the case of linear and elastic materials. The main contribution of this work is an extension of this strategy to material nonlinearities. In particular, we will show that special care must be taken when using the multiparametric strategy in order to take into account damageable laws. This extension is mandatory if one wishes to model the behavior of laminated composites. Indeed, many approaches can be found in the literature for modeling the behavior and simulating the damage state of such materials. These approaches can be divided into three categories depending on the observation scale. Historically, the first models to be used in this context were macroscopic models. These models introduce a posteriori fracture criteria such as those presented by Tsai and Wu (1971) and Azzi and Tsai (1965) or, more recently, in the context of the World Wide Failure Exercise in Hinton et al. (2004). The second category of approaches seeks to represent the degradation of the materials explicitly. There are many such models, known as microscopic models because of the physical observation scale, and several review papers have been dedicated to them such as Nairn and Hu (1994), Nairn (2000), Berthelot (2003). Finally, halfway between these two categories, there is a third approach in which the material is considered on an intermediate scale, called the mesoscale. While these methods no longer represent the degradations explicitly, they seek to take into account their effect on the behavior mainly through damage mechanics Kachanov (1958). In this paper, we will use LMTs (Laboratoire de Mécanique et Technologie) standard damage mesomodel developed in Ladevèze (1986), Ladevèze and Le Dantec (1992), Allix (1992), which belongs to that last category of approaches. This model is summarized in the next section but is not developed in this work. The choice of a mesomodel is a compromise between efficiency and accuracy. Nevertheless, modeling an assembly of composite parts inevitably leads to the definition of a very large nonlinear problem, which cannot be addressed by direct solvers. The resolution of such a problem requires a dedicated strategy. A first category of approaches consists in the homogeneization of micro problems in order to solve a macro one. The micro problems can be computed with a classical finite element method (FE2 (Finite Element square) method (Feyel and Chaboche, 2000; Feyel, 2003)) or other methods such as the Voronoï cell finite element method (Ghosh et al., 1995; Ghosh et al., 2001). Other approaches consists in dividing the problem into several smaller subproblems, which leads to an optimal use of parallel computing architectures.
These methods, known as domain decomposition methods, can be categorized according to the boundary unknowns, which are shared by the subdomains. Primal methods, such as the BDD (Balancing Domain Decomposition) method shown in Le Tallec et al. (1991) and Mandel (1993)), define the relation among the subdomains in terms of displacements. Dual methods, such as the FETI (Finite Element by Tearing and Interconnecting) method shown in Farhat and Roux (1992) and Farhat et al. (2001), seek to express the link in terms of forces. Finally, there is a category of methods, called mixed methods with few examples in Rixen and Farhat (1999), Series et al. (2003), Farhat et al. (2000) and Mandel and Tezaur (1996), which use both primal and dual quantities to connect the subdomains. The introduction of contact and friction into the relations lends itself naturally to the use of mixed methods, such as in Dostal et al. (1998) and Dureisseix and Farhat (2001). The LATIN (for LArge Time INcrement) method, which was introduced in Ladevèze (1985) and Ladevèze (1999), relies on a mixed domain decomposition method. The domain is divided into volume entities, called substructures, which are separated from one another by surface entities, called interfaces. The scalability of the method, which is crucial if a large number of substructures and interfaces is to be considered, is provided by the definition of a macroscopic problem described by Ladevèze and Dureisseix (2000) and developed in the case of cohesive interfaces and elastic plies by Kerfriden et al. (2009). The main addition of this work is to take into account the material non-linearities in the plies, such as damage and plasticity, with a generalization of the LATIN method in the case of non-linear material behavior. After a rapid presentation of the damage mesomodel used in our work, we will give the details of our calculation strategy with non-linear materials and introduce the principle of the multiparametric approach. Finally, we will study two examples in order to illustrate the performance of our method in the case of variable interface or material behavior.
2. Description of the material’s constitutive model For the purpose of this paper, the material’s behavior is modeled using the standard version of the damage mesomodel developed in Ladevèze (1986), Ladevèze and Le Dantec (1992) and Allix (1992)). This model is chosen for its capabilities to reproduce complex damage phenomena such as delamination between plies (Allix et al., 1995; Bordeu et al., 2011) and distribution of damage in the plies of laminate (Ladevèze et al., 1998; Lubineau and Ladevèze, 2008). In these works, comparisons have been carried out with experimental tests that show perfectly consistent results in terms of global behaviour (load–displacement curves) or local behaviour (X-ray radiography on laminates). This model relies on the central idea that the composite consists of two constituents: the ply, which is damageable and inelastic, and the interface, which is elastic and damageable. Each mesoconstituent is defined independently of the others and possesses its own unknowns and its own constitutive relation (see Fig. 1). 2.1. Behavior of the interface The interface is assumed to be subjected to an orthotropic damageable behavior controlled by a single damage variable d. The associated thermodynamic force is defined by:
8 " #a " #a !1=a 9 = <1 2 X hr33 ðsÞi2þ hr3i ðsÞi2þ YðtÞ ¼ sup þ ci 0 0 2 2 ; s6t :2 kn ð1 dðsÞÞ kt ð1 dðsÞÞ i¼1 ð1Þ
2751
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
where indices 1; 2 and 3 denote respectively the direction of the fibers, the direction perpendicular to the fibers, and the normal to the ply. hh:ii and h:iþ designate respectively the mean value through the ply’s thickness and the positive part of the variable. The mixed damage forces are defined by:
Y d ðtÞ ¼ supðY d ðsÞ þ bY d0 ðsÞÞ
ð4Þ
Y F ðtÞ ¼ supðY F ðsÞÞ
ð5Þ
s6t s6t
where b is the coupling between the transverse and shear energies. For small damage rates, the damage variables are governed by:
Fig. 1. Decomposition of the composite into plies and interfaces.
Table 1 The material coefficients of the cohesive interface. Parameter
Value
Variability
0
7120 N mm1
Yes
kt
0
5900 N mm1
No
c1 ¼ c2 a
0.4 1.0 0.0 MPa 0.18 MPa 1 0.01 s
No No No Yes No No
kn
Y0 YC a
sc
where c1 , c2 represent the coupling between the tension and shear stresses in the two transverse directions, and a denotes the brittleness of the material. The damage evolution law depends on the material function x defined by:
"
x¼
n hY Y 0 iþ n þ 1 YC Y0
#n ð2Þ
Then, the evolution law, including what is called a delay effect, can be expressed as a function of x:
(
1
d_ ¼ s0 1 expða0 hx diþ Þ
if d < 1
c
d¼1
ð3Þ
otherwise
The values of the material coefficients, which define the behavior of the interface, are those of the material type studied in Allix (1992). They are given in Table 1. These values will be the same for the two examples presented at the end of the paper.
8 pffiffiffiffiffiffiffiffi pffiffiffiffi0 < x ¼ pYffiffiffiffi d ðtÞ pffiffiffiffiY if d < 1 d d¼ YC Y0 : 1 otherwise ( 0 0 xd ¼ bd if d < 1 and d0 < 1 0 d ¼ 1 otherwise ( 0 if Y F ðtÞ < Y CF dF ¼ 1 otherwise
ð6Þ
ð7Þ ð8Þ
The expressions for xd and xd0 introduce two material parameters: the damage threshold Y 0 and a critical value Y C . The influence of the latter will be studied in Section 5. The extension to large damage rates is carried out traditionally by introducing evolution laws with 0 delay effects for d and d . The last degradation phenomenon of the material’s structure is sliding of the matrix along the fibers with friction. This leads to inelastic strains, which are described in detail in Ladevèze and Le Dantec (1992). In order to introduce these irreversible strains into the model, we seek to describe the behavior in terms of ‘‘effective’’ e and e quantities: r e_ p represent respectively the effective stress and ee the effective strain, which satisfy Trðre_ p Þ ¼ Trð r e_ p Þ with:
r r r re 11 ¼ r11 ; re 12 ¼ 12 ; re 23 ¼ 23 ; re 13 ¼ 13 ; 1d 1d 1d hr i hr i re 22 ¼ 22 þ0 hr22 iþ ; re 33 ¼ 33 þ0 hr33 iþ 1d
1d
Classically, the yield function f satisfies the condition:
2 2 1=2 e ¼ r e ; RÞ e 12 þ r e 223 þ r e 213 þ c2 r e 22 þ r e 233 f ðr R R0 6 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} er eq where c is a material coupling parameter, p is the accumulated plastic strain, R0 is the initial yield stress and R is the strain hardening function defined by:
R ¼ kp pm
ð9Þ
e ¼ p , Relation (9) can be considered With the change of variable p to be linear. Therefore, the hardening function can be considered to be linear with respect to e p and written as: m
2.2. Description of the ply’s behavior Due to the very structure of the ply, its behavior is highly directional and, therefore, anisotropic. Therefore, the damage state of the ply is defined by the three damage variables dF (associated 0 with fiber breakage), d and d . These variables are assumed to be constant throughout the thickness of each ply. The thermodynamic damage forces associated with each variable are:
YF ¼ Yd ¼
1
**
r
2 11 E01
ð1 dF Þ2 ** 1 r212
2 X 3 X i¼1 j>i
r223
mij
r231
þ þ G012 G023 G031 ** ++ hr22 i2þ hr33 i2þ 1 þ Y d0 ¼ 0 E02 E03 ð1 d Þ2 ð1 dÞ2
mji
þ E0i E0j ++
!
++
rii rjj
R ¼ kp e p Finally, the plasticity evolution law must satisfy the following conditions:
@f e e_ p ¼ p_ e @r
and p_ ¼
8 < e @R 1
r eq
:
0
@p
if f ¼ 0 otherwise
In this paper, we will use the material coefficients of Table 2 for the plies, which are identified on the material studied in Ladevèze and Le Dantec (1992). In order to simplify the notations for the remainder of the T pa0 per, let us introduce the ‘‘meta’’ variables X ¼ e and p ; d; d ; dF Y ¼ ðR; Y d ; Y d0 ; Y F ÞT , and rewrite the ply’s equations in the following generic form:
2752
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
Table 2 The material coefficients of the ply (with 1 representing the fiber direction). Parameter
Value
E01 E02 E03
181.5 GPa
No
9.9 GPa
No
9.9 GPa
No
m12 m13 m23
0.34 0.34 0.49 6.16 GPa 6.16 GPa 3.08 GPa 0 MPa
No No No No No No No
0.5 MPa
Yes
40 MPa 0.2 1 0.05 s
No No No No
G12 G13 G23 Y0 YC YF b a
sc
Variability
Fig. 2. The LATIN domain decomposition.
Constitutive relation:
r ¼ Kðd; d0 ; dF Þee
ð10Þ
Y ¼ AX
ð11Þ
Fig. 3. Definition of the interface unknowns.
Evolution laws:
"
e_ p _ X
#
¼B
r
Y
ð12Þ
3. The computational strategy In this section, we will describe the LATIN method in the case of nonlinear material. The computational strategy we chose relies on the multiscale LATIN method developed in Ladevèze (1999) and Ladevèze et al. (2001), which consists in: dividing the global domain into substructures and interfaces in order to define a unique problem for each element resulting from the decomposition; dividing the equations of the problem into two groups: the linear equations and the local equations; solving the problem by means of an iterative calculation scheme. This method has already been applied in Kerfriden et al. (2009) for composites made of orthotropic elastic plies and damageable (cohesive) interfaces. In that context, the decomposition of the domain was natural: the composite’s plies were modeled by the (linear) substructures issued from the decomposition and the cohesive interfaces were identical to the interfaces of the decomposition. We want to preserve this rather intuitive decomposition, but extend the method’s domain of application to damageable materials, which follow the constitutive model of Section 2. Therefore, the nonlinearities will not be concentrated at the interface level, but will also apply to the substructures. 3.1. Decomposition of the domain Assuming small displacements and quasi-static isothermal evolution, the problem is defined for a structure X subjected to volume forces f d , to forces F d over a part @ X2 of its boundary, and to prescribed displacements U d over the complementary part @ X1 of @ X. Here, we consider that the problem has already been decomposed and that each substructure is separated from the adjacent substructures by interfaces (with cohesive or frictional contact behavior), as shown in Fig. 2.
Fig. 4. Graphical representation of the iterative scheme of the LATIN method.
Let us consider a substructure XE resulting from that decomposition. The state of this substructure is defined by a Cauchy stress field rE , an inelastic strain field eEp , and a set of values of the internal variables XE (damage variables and accumulated plastic strain) and their conjugate quantities YE . 0 Let us also consider an interface CEE between two substructures 0
XE and XE . The state of this interface is defined by a pair of force 0
0 _ E; W _ E Þ (see Fig. 3) along fields ðF E ; F E Þ, a pair of velocity fields ðW with the damage variables d in the case of a cohesive interface, and the conjugate quantities Y d . The global state of the structure, denoted s, is defined by the state of each subdomain (interface or substructure).
3.2. The LATIN method The LATIN method used here is based on two key points. The first key point is the separation of the equations of the problem into two groups: the linear equations, whose solutions (which can be local) are denoted Ad , and the local equations, whose solutions (which can be nonlinear) are denoted C. Thus, the exact solution of the global problem is the intersection of Ad and C. The second key point is the introduction of an iterative scheme in order to obtain this exact solution by building successive approximations of the solution s belonging to Ad and to C alternatively. Switching from one solution space to the other requires the introduction of
2753
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
new equations called the search directions (Eþ and E ). The role of the search directions is illustrated in Fig. 4, which represents the evolution of the algorithm in the space of the solutions. Thus, each LATIN iteration consists of two steps: the local step: starting from sn 2 Ad , find bs n such that bs n 2 C and bs n sn 2 Eþ the linear step: starting from bs n 2 C, find snþ1 such that snþ1 2 Ad et snþ1 bs n 2 E The search directions are linear equations and are parameters of the methods, which affect only the convergence rate of the algorithm, but not the final solution of the problem. They involve both the interface unknowns (Eqs. (13, 15)) and the substructure unknowns (Eqs. (14) and (16)). Their expressions are:
8 E
E _ > E E b _ c > ð13Þ W F F ¼ k W > 0 n n n > < n þ ( ) ( ) bs n sn 2 E () rEn > rb En > > ð14Þ > : b E ¼ YE Yn n
UðtÞ ¼ ½NfuðtÞg and eðUðtÞÞ ¼ ½BfuðtÞg On the interfaces, a compatible discretization is applied to the velocity fields and the forces:
_ E b_ E _ E ¼ ½Nfw c _ En g and W W n n ¼ ½Nf w n g E F En ¼ ½NffnE g and b F n ¼ ½Nfbf En g
nþ1
3.3. Determination of the interface problem The interface problem is classically solved during the local step 0
of the LATIN method. Let us consider an interface CEE between two 0
0 _ E; W _ E ; d; Y d Þ substructures XE and XE . At each iteration, ðF E ; F E ; W must satisfy: 0
0
_ E Þ and ðF E0 ; W _ E Þ; the search direction (13) for the pairs ðF E ; W 0
E E 0 the interface equilibrium: 8M 2 CEE ; b F n ðMÞ þ b F n ðMÞ ¼ 0;
the (usually nonlinear) evolution laws, which, if one defines the 0
velocity jump by EE0
8M 2 C ;
_ E _ E c c ¼W n W n , may take the form:
8 9 E
: Y EE ; d;n 0
( ¼R
_ EE0 ðMÞ W n 0 d_ EE
3.5. About discretization A standard finite element discretization is used for the displacement field in the substructures. So the displacement U and the strain e fields have the following form:
8
E > _ E > _ E W c > ð15Þ F Enþ1 b F n ¼ k0 W > nþ1 n > < 8 9 ( ) E snþ1 bs n 2 E () E < b e p;n = > > e_ p;nþ1 ¼ > ð16Þ > > :X : XE bE ; nþ1
_ EE0 W nþ1=2
or linear. The constitutive relation, which depends on internal variables (damage), which are a priori unknown, is nonlinear. The search direction (16), which is part of the constitutive relation (11), enables the calculations to be carried out for a damage value, which is fixed and known (because it was obtained in the local step). Thus, once the damage is known, the constitutive relation becomes a linear relation. Once this difficulty has been eliminated, the separation of the equations is a simple task: the solutions of the kinematic admissibility, equilibrium and constitutive relations (10, 11) belong to Ad , while the equations of the evolution laws (12) belong to C, as shown in Table 3.
The resolution of the linear step considers that the complete state of the interfaces (defined by forces fields, velocity fields and eventually damage variables) is computed during the last local step. The parameter k0 , defined in the Eq. (15), can be considered such as a stiffness between the interface and the substructures. Consequently, a substructure is subjected to mixed boundary conditions. Moreover, during this last local step, the damage state of each substructure is also computed. This leads to a material behavior, which can be considered elastic and linear because the damage state is given. The linear substructure problem takes into account the kinematic admissibility, the linear hebavior at prescribed damage state and the search direction. It becomes, after discretization:
b_ E ðtÞg _ k0 ½HE fuðtÞg þ ½K E fuðtÞg ¼ ½HE bf En1 ðtÞ þ k0 f w n where:
)! ð17Þ
n
The search directions at the interface (13,15) introduce a parameter k0 whose optimal value with respect to the convergence rate of the method was studied in Ladevèze et al. (2001). 3.4. Determination of the substructure problem The mechanical problem to be solved in each substructure XE consists in finding for each time step ðrE ; eEp ; XE YE Þ the evolution, which verifies: the kinematic admissibility; the equilibrium; the damageable behavior of the mesomodel’s plies (Eqs. (10)–(12)) A major difficulty related to the introduction of the nonlinear behavior of the substructures resides in the handling of the constitutive relation (10) in the LATIN algorithm. Indeed, the method is based on the separation of the equations, which are either local
½HE ¼
Z @ XE
½Nt ½NdSE and ½K E ¼
Z XE
½Bt Kd ½BdXE
½K E is the classical finite element stiffness matrix of substructure XE at given damage state and ½HE is the boundary term related to the interfaces.
Table 3 Separation of the substructure’s constitutive relations. Step
Equations to be solved
Linear
In the substructures: Kinematic admissibility Equilibrium r ¼ Kd ee Y ¼ AX Search directions (15, 16)
Local
On the substructures: " # e_ p r ¼B _ Y X At the interfaces: Behavior Search directions (13, 14)
2754
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
For the resolution of this differential equation, an Euler implicit time integration scheme is used. During the local step, the thermodynamical forces associated to the damage state of the cohesive interfaces (defining by the resolution of the Eqs. (1)–(3)) and of the plies (defining by the resolution of the Eqs. (4)–(8)) are given. The use of an Euler implicite time integration scheme leads to a non-linear problem, defined at each Gauss point, which is solved with a Newton–Raphson scheme. 3.6. The multiparametric strategy 3.6.1. Principle The objective of the multiparametric strategy is to solve a large number of mechanically ‘‘similar’’ problems, each associated with a set of parameters. The approach followed in this paper relies on the possibility to initialize the LATIN method using any solution, usually an elastic solution. We will take advantage of this possibility by initializing the resolution of a problem using the result of a previously solved problem whose material parameters are considered to be relatively similar. Thus, a first, mechanically relevant approximation of the solution of the problem is available at no additional computation cost. A graphical representation of this strategy is given by Fig. 5. Following the calculation of a solution associated with problem AAd \ CA , a new set of parameters is considered. In order to simplify the figure, since the variabilities usually concern the nonlinear equations, let us assume that AAd ABd . Thus, the difference between problem A and problem B is contained in the solutions of the local equations. Assuming that the parameters are only slightly different, the space defined by CB is relatively close to that defined by CA . Therefore, to initialize the resolution of problem B, it is usually more relevant to use the solution of problem A than to start from an elastic solution. 3.6.2. Sources of variabilities From a pragmatic standpoint, for designing an assembly of parts made of composite materials, the parameters of the problem can be categorized according to their range of variation: Elastic parameters: the elastic moduli and Poisson coefficients are usually well-known with few identification tests. Therefore, we will consider that a parametric study of their influence is irrelevant. Nonlinear interface parameters: even though the friction coefficients and the parameters involved in the cohesive laws often have a major influence on the solution, they are difficult to identify. Therefore, a parametric study of these parameters is particularly appropriate. Nonlinear substructure parameters: the influence of these parameters on the stiffness matrix is significant, but often very local, which makes them difficult to identify. Therefore, it may be very relevant to resort to a parametric study to identify them. We performed two parametric studies related to the last two groups of parameters. The first example concerns the case in which the interface parameters vary; the second example addresses the case in which the parameters affect the laws governing the substructures. 4. Application to a cohesive interface problem 4.1. Definition of the problem This first example assumes that the equations of the substructures are deterministic and that the variabilities concern the inter-
Fig. 5. Graphical interpretation of the multiparametric method.
faces alone. The objective of this section is to apply the multiparametric strategy to a simple case including the cohesive interfaces of the standard damage mesomodel. The example concerns the opening of a laminate subjected to a normal compression load. The laminate consisted of 5 tridimensional plies, oriented in the 0° and 90° directions alternatively, along with titanium blocks whose behavior was assumed to be linear, elastic and isotropic (see Fig. 6). 0 The multiparametric study was performed on parameters kn and Y 0 , which correspond to the initial stiffness of the interface and to the damage threshold and were assumed to be equal at each interface. Indeed, during a identification test, the value of these parameters are difficult to determinate and cannot be well known. The range of variation of each parameter was discretized into 26 arbitrarily values, leading to 676 different problems (see Table 4).
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
2755
Fig. 6. A view of the laminate – the behavior of a cohesive interface.
The efficiency of the strategy is expressed by the gain G defined by:
Table 4 0 The values of parameters Y C and kn . Parameter
Min. value
Max. value
Step
YC
0.1 MPa 103 N/mm
1.5 MPa 107 N/mm
0.056 MPa
0
kn a
G¼
s625 625 s1 ’ l625 l625
a
For all these non-linear analysis, the gain is computed considering that:
Denotes a constant step in the logarithmic space.
for the resolution of the first problem, we use the LATIN method to solve the non-linear problem (associated to CPU time s1 ), instead of an usual Finite Element software, the LATIN method is used to obtain the CPU time reference s1 ; however they are equivalent as showed in Roulet et al. (2011) with a comparison between the commercial software ABAQUS and the LATIN method, the CPU time associated to the 625 resolutions is estimated, assuming that the CPU time is constant for each non-linear analysis.
Fig. 7. Influence of the parameters on the opening of the laminate.
Table 5 CPU time and gain achieved with the multiparametric strategy. Number of processors
s1
l625
Gain
1 6
32 min 10 min
1098 min 330 min
18 19 0
Let us point out that on the figure the variation range of kn is expressed in logarithmic scale, and that it is on that scale that the points were regularly spaced.
The gain shown in Table 5 is about 18 for both parallel and sequential computations. Thus, despite its ease of implementation, the multiparametric strategy leads to a significant improvement in the computation time required for a parametric study. In this example, the nonlinearities were concentrated only at the interfaces. The substructures were considered to have linear elastic behavior. Therefore, the state of the complete structure was defined by the values of the fields at the interfaces, and the values of the fields within the substructures could be obtained directly through an elastic calculation with prescribed boundaries. Thus, only the interface data needed to be stored for each problem and transferred from one problem to another. The introduction of nonlinearities on the substructure level precludes the use of this trick and raises the crucial question of the storage of the solutions, which is addressed by the second example.
4.2. Results 5. Application to damageable plies The parameters of the cohesive laws play an important role in the opening of each ply. Indeed, Fig. 7 shows the opening of the laminate (defined as the difference in the end displacements of the upper and lower plies) as a function of the parameters. Let us designate by:
s1 , the CPU time required for the resolution of the first problem; s625 , the CPU time required for the 625 resolutions, here estimated to be 625 s1 ; l625 , the CPU time required for the 625 resolutions using the
multiparametric strategy.
For this example, the parameters involved in the evolution of the damage were allowed to vary during the parametric study. As explained above, the introduction of nonlinearities in the substructures leads to the question of the storage of the solutions. Indeed, while in the elastic case the fields at the interfaces were sufficient, now this is no longer the case. Nevertheless, we will try to perform a calculation initialized by the previously calculated interface fields alone. This calculation will be compared to a calculation whose initialization took into account the damage states of all the substructures.
2756
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
5.1. Description of the problem In the problem considered here, the damage state of the material influences the global behavior of the assembly drastically. The objective was to carry out a 3D simulation of a 4-point bending test (see Fig. 8). The three plates were made of a damageable laminate and assembled by elastic bolts. The loading was applied in two steps: a preloading step: a prescribed relative displacement U i between the head and the thread of each bolt; a bending step: a prescribed normal displacement at the head of each bolt to induce the bending of the assembly.
5.2. Parametric study
Table 6 Variations of the parameters. Parameter
Min. value
Max. value
step
Ui
0.01 mm 0.1 MPa
0.1 mm 1 MPa
0.01 mm 0.1 MPa
YC
Table 7 The memory storage requirements for the solution of one problem. Initialization from . . .
Memory requirement
the interface fields alone the interface fields and the internal variables of the substructures Stiffness matrices (for one time step)
1069.7 kB 2661.4 kB 32,944.6 kB
A parametric study was carried out for two parameters: the preload in the bolts U i and the critical value Y C . These parameters have a non-negligible influence on the solution and are often unknown. Each parameter was allowed 11 arbitrarily different values (see Table 6), leading to 121 problems. With regard to the richness of the initialization, which conditions the storage space required for each solution, we considered two approaches: an initialization from the interface fields alone, thus assuming no initial damage; an initialization from the interface fields and the damage state of the substructures. Table 7 gives an idea of the memory required to store a solution corresponding to the mesh used for this example. For each Gauss point of each interface element, 7 floating-point numbers are necessary (3 for the displacements, 3 for the loads, and one for the damage variable); for each Gauss point in the substructures, 3 0 numbers are necessary (d; d and dF ). In other words, for the discretization considered here (linear elements, 5 elements through the thickness of each ply), taking into account the damage state in each Gauss point of each substructure amounts to doubling the storage space required. While storing the matrices is not necessary in order to use the multiparametric strategy, this enables one to compare the amount of RAM required for the calculation and the disk space required for storing the results.
5.3. Results The relevant quantity of interest for an easy comparison of the solutions is necessarily a function of the damage state of the structure. Therefore, we chose to represent, in Fig. 9, the reaction forces over the surface where the normal displacements were prescribed. The study of the relevance of the two types of initialization is summarized in Table 8. One can see that for the multiparametric strategy to be relevant one must reinstate the damage state of the previous calculation right from the initialization. The informa-
Fig. 9. Influence of the parameters on the reaction force.
Table 8 The gain observed in the parametric study. Initialization
1 Resolution
121 Resolutions
Gain
Interfaces Interfaces & SST
810 s 810 s
95,590 s 7159 s
1.02 13.7
tion on the interface level is insufficiently rich, so the gain is practically inexistent. 6. Conclusion We presented an efficient computational strategy for the resolution of multiple nonlinear problems, especially when these problems involve damage. This strategy relies on a domain decomposition, which makes parallelization easy, and on the LATIN iterative solver. A relevant initialization of the calculation by reusing previously computed solutions leads to a significant reduction in the number of LATIN iterations, which are necessary, and, therefore, in the computation time. In the case of variabilities in the interface laws, only the interface fields are stored. In the case of nonlinear substructures with variable behavior, storing the state of the internal variables is necessary in order to preserve significant gains in computation time. The two examples given illustrate the efficiency of the multiparametric strategy as long as these initialization conditions are respected. Acknowledgement
Fig. 8. The geometry of the 4-point bending problem.
The research leading to these results was funded by the European Community’s Seventh Framework Program FP7/2007–2013 under Grant Agreement No. 213371 (http://www.maaximus.eu/).
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
References Alart, P., Curnier, A., 1991. A mixed formulation for frictional contact problems prone to newton like solution methods. Computer Methods in Applied Mechanics and Engineering 92 (3), 353–375. Allix, O., 1992. Damage analysis of delamination around a hole. New Advances in Computational Structural Mechanics, pp. 411–421. Allix, O., Ladevèze, P., Corigliano, A., 1995. Damage analysis of interlaminar fracture specimens. Composite Structures 31 (1), 61–74. Armero, F., Petcz, E., 1998. Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems. Computer Methods in Applied Mechanics and Engineering 158 (3–4), 269–300. Arora, J.S., Chahande, A.I., Paeng, J.K., 1991. Multiplier methods for engineering optimization. International Journal for Numerical Methods in Engineering 32 (7), 1485–1525. Azzi, V.D., Tsai, S.W., 1965. Anisotropic strength of composites. Experimental Mechanics 5 (9), 283–288. Berthelot, J.-M., 2003. Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: static and fatigue loading. Applied Mechanics Reviews 56 (1), 111–147. Bordeu, F., Boucard, P., Lubineau, G., 2011. A mesoscale model for damage, cracking and delamination prediction in composite materials. Composite Materials 17 (4), 271–282. Boucard, P.A., Champaney, L., 2003. A suitable computational strategy for the parametric analysis of problems with multiple contact. International Journal for Numerical Methods in Engineering 57 (9), 1259–1281. Chabrand, P., Dubois, F., Raous, M., 1998. Various numerical methods for solving unilateral contact problems with friction. Mathematical and Computer Modelling 28 (4–8), 97–108, recent Advances in Contact Mechanics. Champaney, L., 2011. A domain decomposition method for studying the effects of missing fasteners on the behavior of structural assemblies with contact and friction. Computer Methods in Applied Mechanics and Engineering (205/208), 121–129. http://dx.doi.org/10.1016/j.cma.2011.04.008. Champaney, L., Boucard, P.A., Guinard, S., 2008. Adaptive multi-analysis strategy for contact problems with friction. Computational Mechanics 42 (2), 305–315. Champaney, L., Cognard, J.Y., Ladevèze, P., 1999. Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions. Computers & Structures 73 (1–5), 249–266. Dostal, Z., Friedlander, A., Santos, S.A., 1998. Solution of coercive and semicoercive contact problems by feti domain decomposition. Contemporary Mathematics 218, 82–93. Dureisseix, D., Farhat, C., 2001. A numerically scalable domain decomposition method for the solution of frictionless contact problems. International Journal for Numerical Methods in Engineering 50 (12), 2643–2666. Farhat, C., Lesoinne, M., LeTallec, P., Pierson, K., Rixen, D., 2001. Feti-dp: a dualprimal unified feti method-part i: a faster alternative to the two-level feti method. International Journal for Numerical Methods in Engineering 50 (7), 1523–1544. Farhat, C., Lesoinne, M., Pierson, K., 2000. A scalable dual-primal domain decomposition method. Numerical linear algebra with applications 7 (7–8), 687–714. Farhat, C., Roux, F.X., 1992. An unconventional domain decomposition method for an efficient parallel solution of large-scale finite element systems. SIAM Journal on Scientific and Statistical Computing 13, 379. Feyel, F., 2003. A multilevel finite element method (FE2) to describe the response of highly non-linear structures using generalized continua. Computer Methods in Applied Mechanics and Engineering 192 (28–30), 3233–3244. Feyel, F., Chaboche, J.-L., 2000. Fe2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials. Computer Methods in Applied Mechanics and Engineering 183 (3–4), 309–330. Ghanem, R.G., 1999. Ingredients for a general purpose stochastic finite elements implementation. Computer Methods in Applied Mechanics and Engineering 168, 19–34. Ghanem, R.G., Kruger, R.M., 1996. Numerical solution of spectral stochastic finite element systems. Computer Methods in Applied Mechanics and Engineering 129 (3), 289–303. Ghanem, R.G., Spanos, P.D., 1991. Stochastic Finite Elements: A Spectral Approach. Springer, Berlin. Ghosh, S., Lee, K., Moorthy, S., 1995. Multiple scale analysis of heterogeneous elastic structures using homogenization theory and voronoi cell finite element method. International Journal of Solids and Structures 32 (1), 27–62. Ghosh, S., Lee, K., Raghavan, P., 2001. A multi-level computational model for multiscale damage analysis in composite and porous materials. International Journal of Solids and Structures 38 (14), 2335–2385.
2757
Hinton, M.J., Kaddour, A.S., Soden, P.D., 2004. Failure criteria in fibre reinforced polymer composites: the world-wide failure exercise. Elsevier Science. Kachanov, L.M., 1958. On creep rupture time. Izv. Acad. Nauk SSSR, Otd. Techn. Nauk 8, 26–31. Kaminski, M., 2005. Computational Mechanics of Composite Materials: Sensitivity, Randomness and Multiscale Behaviour. Springer. Kerfriden, P., Allix, O., Gosselet, P., 2009. A three-scale domain decomposition method for the 3D analysis of debonding in laminates. Computational Mechanics 44 (3), 343–362. Kikuchi, N., 1982. Penalty/finite element approximations of a class of unilateral contact problems. Penalty method and finite element method. ASME, New York. Kikuchi, N., Oden, J.T., 1988. Contact problems in elasticity: a study of variational inequalities and finite element methods. Society for Industrial Mathematics. Klarbring, A., 1992. Mathematical programming and augmented Lagrangian methods for frictional contact problems. In: Proceedings of Contact Mechanics International Symposium, PPUR, vol. 369, p. 390. Kleiber, M., Hien, T.D., 1992. The Stochastic Finite Element Method. John Wiley & Sons, New York. Ladevèze, P., 1985. Sur une famille d’algorithmes en mécanique des structures. Compte rendu de l’académie des Sciences 300 (2), 41–44. Ladevèze, P., 1986. Sur la mécanique de l’endommagement des composites. C. Bathias et D. Menkes, éds: Comptes-rendus des JNC5, Paris. Ladevèze, P., 1999. Nonlinear Computational Structural Mechanics: New Approaches and Non-incremental Methods of Calculation. Springer Verlag. Ladevèze, P., Allix, O., Gornet, L., Lévèque, D., Perret, L., 1998. A computational damage mechanics approach for laminates: identification and comparison with experimental results. In: George Z. Voyiadjis, J.-W.W.J., Chaboche, J.-L. (Eds.), Damage Mechanics in Engineering Materials, vol. 46 of Studies in Applied Mechanics. Elsevier, pp. 481–500. Ladevèze, P., Dureisseix, D., 2000. A micro/macro approach for parallel computing of heterogeneous structures. International Journal for Computational Civil and Structural Engineering 1 (1), 18–28. Ladevèze, P., Le Dantec, E., 1992. Damage modelling of the elementary ply for laminated composites. Composites Science and Technology 43 (3), 257–267. Ladevèze, P., Loiseau, O., Dureisseix, D., 2001. A micro-macro and parallel computational strategy for highly heterogeneous structures. International Journal for Numerical Methods in Engineering 52 (1–2), 121–138. Le Tallec, P., Roeck, Y.H., Vidrascu, M., 1991. Domain decomposition methods for large linearly elliptic three-dimensional problems. Journal of Computational and Applied Mathematics 34 (1), 93–117. Lubineau, G., Ladevèze, P., 2008. Construction of a micromechanics-based intralaminar mesomodel, and illustrations in abaqus/standard. Computational Materials Science 43 (1), 137–145. Mandel, J., 1993. Balancing domain decomposition. Communications in Numerical Methods in Engineering 9 (3), 233–241. Mandel, J., Tezaur, R., 1996. Convergence of a substructuring method with lagrange multipliers. Numerische Mathematik 73 (4), 473–487. Nairn, J.A., 2000. Matrix microcracking in composites. Polymer Matrix Composites 2, 403–432. Nairn, J.A., Hu, S., 1994. Matrix microcracking. Composite Materials Series, pp. 187– 243. Radi, B., Baba, O.A., Gelin, J.C., 1998. Treatment of the frictional contact via a Lagrangian formulation. Mathematical and Computer Modelling 28 (4–8), 407– 412, recent Advances in Contact Mechanics. Rixen, D.J., Farhat, C., 1999. A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems. International Journal for Numerical Methods in Engineering 44 (4), 489–516. Roulet, V., Champaney, L., Boucard, P.A., 2011. A parallel strategy for the multiparametric analysis of structures with large contact and friction surfaces. Advances in Engineering Software 42 (6), 347–358. Series, L., Feyel, F., Roux, F.X., 2003. Une méthode de décomposition de domaine avec deux multiplicateurs de lagrange. Simo, J.C., Laursen, T.A., 1992. An augmented lagrangian treatment of contact problems involving friction. Computers & Structures 42 (1), 97–116. Tsai, S.W., Wu, E.M., 1971. A general theory of strength for anisotropic materials. Journal of Composite Materials 5 (1), 58. Wriggers, P., 1995. Finite element algorithms for contact problems. Archives of Computational Methods in Engineering 2 (4), 1–49. Wriggers, P., Simo, J.C., Taylor, R.L., 1985. Penalty and augmented lagrangian formulations for contact problems 85, 97–106. Zhong, Z., Mackerle, J., 1992. Static contact problems – a review. Engineering Computations 9 (3–37).
Contents lists available at SciVerse ScienceDirect
International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr
An efficient computational strategy for composite laminates assemblies including variability V. Roulet, P.-A. Boucard ⇑, L. Champaney LMT-Cachan (ENS Cachan/CNRS/UPMC/PRES UniverSud Paris), 61 av. du Président Wilson, F-94230 Cachan, France
a r t i c l e
i n f o
Article history: Received 6 September 2012 Received in revised form 29 March 2013 Available online 10 May 2013 Keywords: Multiparametric strategy LATIN solver Parallel computation Laminate composites Damage
a b s t r a c t The aim of this work is to present an efficient numerical strategy for studying the influence of the material parameters on problems involving 3D assemblies of composite parts with contact and friction. This approach is based on the multiscale LATIN method with domain decomposition and its ability to reuse the solution of one problem (for a particular set of parameters) to solve many similar problems, each associated with another set of parameters. The presentation of the proposed calculation method will be followed by two sample applications, one with variabilities in the cohesive interface laws and the other with variabilities in the damageable behavior laws, which characterize the composite material. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction In order to reduce design and sizing costs, manufacturing industries are seeking to limit the number of actual tests and replace them by numerical simulations. However, this approach may not be entirely compatible with the recent adoption by aeronautical industries of composite materials, whose behavior is complex and often imperfectly known. Taking these material uncertainties into account leads to the introduction of variabilities into the constitutive model. Starting with the hypothesis that the material properties of a composite material can be modeled within the framework of the probability theory, a very rich mathematical background is available to completely characterize the probabilistic behavior and the evolution of the structure under an external perturbation. When dealing with complex systems, approximate methods and numerical integration are necessary. If the uncertainties are small, the perturbation theory is a very valuable tool for analyzing their effect (see Kleiber and Hien, 1992). If they are large, then the perturbation theory is not applicable. One obvious way to solve Stochastic Partial Differential Equation (SPDEs) is the Monte Carlo method. This method is expensive, especially if higher-order accuracy for mean values, standard variations, queues of distributions, etc. is sought. Various numerical methods for solving SPDEs have been proposed in the literature in Ghanem and Kruger (1996).
⇑ Corresponding author. E-mail addresses: [email protected] (V. Roulet), boucard@lmt. ens-cachan.fr, [email protected] (P.-A. Boucard), champa [email protected] (L. Champaney). 0020-7683/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijsolstr.2013.04.028
Ghanem and Spanos (1991), and Ghanem (1999) propose a hybrid finite element/spectral approach called the Spectral Stochastic Finite Element Method (SSFEM). An extended review of these methods is available in Kaminski (2005). We propose here an alternative approach to reduce the computational costs. Furthermore, the simulation of assemblies of composite parts raises two difficulties: the interactions among the components, which are often complex and nonlinear with frictional contact problems, and the material’s behavior, especially its degradation. Contact problems are characterized by constraints such as nonpenetration conditions, and an active area of contact – that is, an area where contact effectively occurs – that is unknown a priori. For these reasons, these problems lead to stiff non-linear system of equations. Several approaches exist for solving static contact problems (Kikuchi and Oden, 1988; Zhong and Mackerle, 1992; Wriggers, 1995). In most of them, the numerical methods that are employed for enforcing the contact constraints can be grouped into Lagrange multiplier and penalty methods (Wriggers et al., 1985). The penalty methods, used in Kikuchi (1982) and Armero and Petcz (1998), are closely related to the regularization of the contact constraints. They are usually formulated in terms of the displacement variables and therefore are primal methods. They allow treating contact as a material behaviour, as exemplified by the method of joint finite elements in Alart and Curnier (1991). Penalty methods can experience various numerical difficulties, especially ill-conditioning, when a too large or too small penalty parameter is introduced. Lagrange multiplier methods are dual methods where the multipliers, which represent the contact reaction forces, are introduced in order to enforce exactly the non-penetration conditions. Augmented Lagrange multiplier methods, used in Arora et
2750
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
al. (1991), Klarbring (1992), Radi et al. (1998), Chabrand et al. (1998), result in mixed formulations involving both displacement and force unknowns. The numerical solution schemes underlying both the Lagrange multiplier and augmented Lagrange multiplier methods are often related to the Uzawa algorithm (see Simo and Laursen, 1992; Champaney et al., 1999). Consequently, it is necessary to develop an efficient numerical strategy dedicated to the resolution of nonlinear problems involving very large numbers of degrees of freedom. Variabilities may be introduced either at the level of the interactions among parts (friction coefficients, preloads, etc.) or in the constitutive material model itself. Thus, the solution of the problem is no longer the result of a deterministic calculation, but must be viewed as a response function of the input parameters. In the rest of this work, we will assume that the variation range of each parameter is divided into a finite number of values. Thus, each set of the parameters defines an associated problem. Then, a parametric study consists in solving as many problems as there are parameter sets. Therefore, one must be able to solve a very large number of ‘‘similar’’ problems. In order to do that, we will use a multiparametric strategy, which was already presented in Boucard and Champaney (2003), Champaney et al. (2008) and Champaney (2011) for assemblies with contact and friction, but only in the case of linear and elastic materials. The main contribution of this work is an extension of this strategy to material nonlinearities. In particular, we will show that special care must be taken when using the multiparametric strategy in order to take into account damageable laws. This extension is mandatory if one wishes to model the behavior of laminated composites. Indeed, many approaches can be found in the literature for modeling the behavior and simulating the damage state of such materials. These approaches can be divided into three categories depending on the observation scale. Historically, the first models to be used in this context were macroscopic models. These models introduce a posteriori fracture criteria such as those presented by Tsai and Wu (1971) and Azzi and Tsai (1965) or, more recently, in the context of the World Wide Failure Exercise in Hinton et al. (2004). The second category of approaches seeks to represent the degradation of the materials explicitly. There are many such models, known as microscopic models because of the physical observation scale, and several review papers have been dedicated to them such as Nairn and Hu (1994), Nairn (2000), Berthelot (2003). Finally, halfway between these two categories, there is a third approach in which the material is considered on an intermediate scale, called the mesoscale. While these methods no longer represent the degradations explicitly, they seek to take into account their effect on the behavior mainly through damage mechanics Kachanov (1958). In this paper, we will use LMTs (Laboratoire de Mécanique et Technologie) standard damage mesomodel developed in Ladevèze (1986), Ladevèze and Le Dantec (1992), Allix (1992), which belongs to that last category of approaches. This model is summarized in the next section but is not developed in this work. The choice of a mesomodel is a compromise between efficiency and accuracy. Nevertheless, modeling an assembly of composite parts inevitably leads to the definition of a very large nonlinear problem, which cannot be addressed by direct solvers. The resolution of such a problem requires a dedicated strategy. A first category of approaches consists in the homogeneization of micro problems in order to solve a macro one. The micro problems can be computed with a classical finite element method (FE2 (Finite Element square) method (Feyel and Chaboche, 2000; Feyel, 2003)) or other methods such as the Voronoï cell finite element method (Ghosh et al., 1995; Ghosh et al., 2001). Other approaches consists in dividing the problem into several smaller subproblems, which leads to an optimal use of parallel computing architectures.
These methods, known as domain decomposition methods, can be categorized according to the boundary unknowns, which are shared by the subdomains. Primal methods, such as the BDD (Balancing Domain Decomposition) method shown in Le Tallec et al. (1991) and Mandel (1993)), define the relation among the subdomains in terms of displacements. Dual methods, such as the FETI (Finite Element by Tearing and Interconnecting) method shown in Farhat and Roux (1992) and Farhat et al. (2001), seek to express the link in terms of forces. Finally, there is a category of methods, called mixed methods with few examples in Rixen and Farhat (1999), Series et al. (2003), Farhat et al. (2000) and Mandel and Tezaur (1996), which use both primal and dual quantities to connect the subdomains. The introduction of contact and friction into the relations lends itself naturally to the use of mixed methods, such as in Dostal et al. (1998) and Dureisseix and Farhat (2001). The LATIN (for LArge Time INcrement) method, which was introduced in Ladevèze (1985) and Ladevèze (1999), relies on a mixed domain decomposition method. The domain is divided into volume entities, called substructures, which are separated from one another by surface entities, called interfaces. The scalability of the method, which is crucial if a large number of substructures and interfaces is to be considered, is provided by the definition of a macroscopic problem described by Ladevèze and Dureisseix (2000) and developed in the case of cohesive interfaces and elastic plies by Kerfriden et al. (2009). The main addition of this work is to take into account the material non-linearities in the plies, such as damage and plasticity, with a generalization of the LATIN method in the case of non-linear material behavior. After a rapid presentation of the damage mesomodel used in our work, we will give the details of our calculation strategy with non-linear materials and introduce the principle of the multiparametric approach. Finally, we will study two examples in order to illustrate the performance of our method in the case of variable interface or material behavior.
2. Description of the material’s constitutive model For the purpose of this paper, the material’s behavior is modeled using the standard version of the damage mesomodel developed in Ladevèze (1986), Ladevèze and Le Dantec (1992) and Allix (1992)). This model is chosen for its capabilities to reproduce complex damage phenomena such as delamination between plies (Allix et al., 1995; Bordeu et al., 2011) and distribution of damage in the plies of laminate (Ladevèze et al., 1998; Lubineau and Ladevèze, 2008). In these works, comparisons have been carried out with experimental tests that show perfectly consistent results in terms of global behaviour (load–displacement curves) or local behaviour (X-ray radiography on laminates). This model relies on the central idea that the composite consists of two constituents: the ply, which is damageable and inelastic, and the interface, which is elastic and damageable. Each mesoconstituent is defined independently of the others and possesses its own unknowns and its own constitutive relation (see Fig. 1). 2.1. Behavior of the interface The interface is assumed to be subjected to an orthotropic damageable behavior controlled by a single damage variable d. The associated thermodynamic force is defined by:
8 " #a " #a !1=a 9 = <1 2 X hr33 ðsÞi2þ hr3i ðsÞi2þ YðtÞ ¼ sup þ ci 0 0 2 2 ; s6t :2 kn ð1 dðsÞÞ kt ð1 dðsÞÞ i¼1 ð1Þ
2751
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
where indices 1; 2 and 3 denote respectively the direction of the fibers, the direction perpendicular to the fibers, and the normal to the ply. hh:ii and h:iþ designate respectively the mean value through the ply’s thickness and the positive part of the variable. The mixed damage forces are defined by:
Y d ðtÞ ¼ supðY d ðsÞ þ bY d0 ðsÞÞ
ð4Þ
Y F ðtÞ ¼ supðY F ðsÞÞ
ð5Þ
s6t s6t
where b is the coupling between the transverse and shear energies. For small damage rates, the damage variables are governed by:
Fig. 1. Decomposition of the composite into plies and interfaces.
Table 1 The material coefficients of the cohesive interface. Parameter
Value
Variability
0
7120 N mm1
Yes
kt
0
5900 N mm1
No
c1 ¼ c2 a
0.4 1.0 0.0 MPa 0.18 MPa 1 0.01 s
No No No Yes No No
kn
Y0 YC a
sc
where c1 , c2 represent the coupling between the tension and shear stresses in the two transverse directions, and a denotes the brittleness of the material. The damage evolution law depends on the material function x defined by:
"
x¼
n hY Y 0 iþ n þ 1 YC Y0
#n ð2Þ
Then, the evolution law, including what is called a delay effect, can be expressed as a function of x:
(
1
d_ ¼ s0 1 expða0 hx diþ Þ
if d < 1
c
d¼1
ð3Þ
otherwise
The values of the material coefficients, which define the behavior of the interface, are those of the material type studied in Allix (1992). They are given in Table 1. These values will be the same for the two examples presented at the end of the paper.
8 pffiffiffiffiffiffiffiffi pffiffiffiffi0 < x ¼ pYffiffiffiffi d ðtÞ pffiffiffiffiY if d < 1 d d¼ YC Y0 : 1 otherwise ( 0 0 xd ¼ bd if d < 1 and d0 < 1 0 d ¼ 1 otherwise ( 0 if Y F ðtÞ < Y CF dF ¼ 1 otherwise
ð6Þ
ð7Þ ð8Þ
The expressions for xd and xd0 introduce two material parameters: the damage threshold Y 0 and a critical value Y C . The influence of the latter will be studied in Section 5. The extension to large damage rates is carried out traditionally by introducing evolution laws with 0 delay effects for d and d . The last degradation phenomenon of the material’s structure is sliding of the matrix along the fibers with friction. This leads to inelastic strains, which are described in detail in Ladevèze and Le Dantec (1992). In order to introduce these irreversible strains into the model, we seek to describe the behavior in terms of ‘‘effective’’ e and e quantities: r e_ p represent respectively the effective stress and ee the effective strain, which satisfy Trðre_ p Þ ¼ Trð r e_ p Þ with:
r r r re 11 ¼ r11 ; re 12 ¼ 12 ; re 23 ¼ 23 ; re 13 ¼ 13 ; 1d 1d 1d hr i hr i re 22 ¼ 22 þ0 hr22 iþ ; re 33 ¼ 33 þ0 hr33 iþ 1d
1d
Classically, the yield function f satisfies the condition:
2 2 1=2 e ¼ r e ; RÞ e 12 þ r e 223 þ r e 213 þ c2 r e 22 þ r e 233 f ðr R R0 6 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} er eq where c is a material coupling parameter, p is the accumulated plastic strain, R0 is the initial yield stress and R is the strain hardening function defined by:
R ¼ kp pm
ð9Þ
e ¼ p , Relation (9) can be considered With the change of variable p to be linear. Therefore, the hardening function can be considered to be linear with respect to e p and written as: m
2.2. Description of the ply’s behavior Due to the very structure of the ply, its behavior is highly directional and, therefore, anisotropic. Therefore, the damage state of the ply is defined by the three damage variables dF (associated 0 with fiber breakage), d and d . These variables are assumed to be constant throughout the thickness of each ply. The thermodynamic damage forces associated with each variable are:
YF ¼ Yd ¼
1
**
r
2 11 E01
ð1 dF Þ2 ** 1 r212
2 X 3 X i¼1 j>i
r223
mij
r231
þ þ G012 G023 G031 ** ++ hr22 i2þ hr33 i2þ 1 þ Y d0 ¼ 0 E02 E03 ð1 d Þ2 ð1 dÞ2
mji
þ E0i E0j ++
!
++
rii rjj
R ¼ kp e p Finally, the plasticity evolution law must satisfy the following conditions:
@f e e_ p ¼ p_ e @r
and p_ ¼
8 < e @R 1
r eq
:
0
@p
if f ¼ 0 otherwise
In this paper, we will use the material coefficients of Table 2 for the plies, which are identified on the material studied in Ladevèze and Le Dantec (1992). In order to simplify the notations for the remainder of the T pa0 per, let us introduce the ‘‘meta’’ variables X ¼ e and p ; d; d ; dF Y ¼ ðR; Y d ; Y d0 ; Y F ÞT , and rewrite the ply’s equations in the following generic form:
2752
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
Table 2 The material coefficients of the ply (with 1 representing the fiber direction). Parameter
Value
E01 E02 E03
181.5 GPa
No
9.9 GPa
No
9.9 GPa
No
m12 m13 m23
0.34 0.34 0.49 6.16 GPa 6.16 GPa 3.08 GPa 0 MPa
No No No No No No No
0.5 MPa
Yes
40 MPa 0.2 1 0.05 s
No No No No
G12 G13 G23 Y0 YC YF b a
sc
Variability
Fig. 2. The LATIN domain decomposition.
Constitutive relation:
r ¼ Kðd; d0 ; dF Þee
ð10Þ
Y ¼ AX
ð11Þ
Fig. 3. Definition of the interface unknowns.
Evolution laws:
"
e_ p _ X
#
¼B
r
Y
ð12Þ
3. The computational strategy In this section, we will describe the LATIN method in the case of nonlinear material. The computational strategy we chose relies on the multiscale LATIN method developed in Ladevèze (1999) and Ladevèze et al. (2001), which consists in: dividing the global domain into substructures and interfaces in order to define a unique problem for each element resulting from the decomposition; dividing the equations of the problem into two groups: the linear equations and the local equations; solving the problem by means of an iterative calculation scheme. This method has already been applied in Kerfriden et al. (2009) for composites made of orthotropic elastic plies and damageable (cohesive) interfaces. In that context, the decomposition of the domain was natural: the composite’s plies were modeled by the (linear) substructures issued from the decomposition and the cohesive interfaces were identical to the interfaces of the decomposition. We want to preserve this rather intuitive decomposition, but extend the method’s domain of application to damageable materials, which follow the constitutive model of Section 2. Therefore, the nonlinearities will not be concentrated at the interface level, but will also apply to the substructures. 3.1. Decomposition of the domain Assuming small displacements and quasi-static isothermal evolution, the problem is defined for a structure X subjected to volume forces f d , to forces F d over a part @ X2 of its boundary, and to prescribed displacements U d over the complementary part @ X1 of @ X. Here, we consider that the problem has already been decomposed and that each substructure is separated from the adjacent substructures by interfaces (with cohesive or frictional contact behavior), as shown in Fig. 2.
Fig. 4. Graphical representation of the iterative scheme of the LATIN method.
Let us consider a substructure XE resulting from that decomposition. The state of this substructure is defined by a Cauchy stress field rE , an inelastic strain field eEp , and a set of values of the internal variables XE (damage variables and accumulated plastic strain) and their conjugate quantities YE . 0 Let us also consider an interface CEE between two substructures 0
XE and XE . The state of this interface is defined by a pair of force 0
0 _ E; W _ E Þ (see Fig. 3) along fields ðF E ; F E Þ, a pair of velocity fields ðW with the damage variables d in the case of a cohesive interface, and the conjugate quantities Y d . The global state of the structure, denoted s, is defined by the state of each subdomain (interface or substructure).
3.2. The LATIN method The LATIN method used here is based on two key points. The first key point is the separation of the equations of the problem into two groups: the linear equations, whose solutions (which can be local) are denoted Ad , and the local equations, whose solutions (which can be nonlinear) are denoted C. Thus, the exact solution of the global problem is the intersection of Ad and C. The second key point is the introduction of an iterative scheme in order to obtain this exact solution by building successive approximations of the solution s belonging to Ad and to C alternatively. Switching from one solution space to the other requires the introduction of
2753
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
new equations called the search directions (Eþ and E ). The role of the search directions is illustrated in Fig. 4, which represents the evolution of the algorithm in the space of the solutions. Thus, each LATIN iteration consists of two steps: the local step: starting from sn 2 Ad , find bs n such that bs n 2 C and bs n sn 2 Eþ the linear step: starting from bs n 2 C, find snþ1 such that snþ1 2 Ad et snþ1 bs n 2 E The search directions are linear equations and are parameters of the methods, which affect only the convergence rate of the algorithm, but not the final solution of the problem. They involve both the interface unknowns (Eqs. (13, 15)) and the substructure unknowns (Eqs. (14) and (16)). Their expressions are:
8 E
E _ > E E b _ c > ð13Þ W F F ¼ k W > 0 n n n > < n þ ( ) ( ) bs n sn 2 E () rEn > rb En > > ð14Þ > : b E ¼ YE Yn n
UðtÞ ¼ ½NfuðtÞg and eðUðtÞÞ ¼ ½BfuðtÞg On the interfaces, a compatible discretization is applied to the velocity fields and the forces:
_ E b_ E _ E ¼ ½Nfw c _ En g and W W n n ¼ ½Nf w n g E F En ¼ ½NffnE g and b F n ¼ ½Nfbf En g
nþ1
3.3. Determination of the interface problem The interface problem is classically solved during the local step 0
of the LATIN method. Let us consider an interface CEE between two 0
0 _ E; W _ E ; d; Y d Þ substructures XE and XE . At each iteration, ðF E ; F E ; W must satisfy: 0
0
_ E Þ and ðF E0 ; W _ E Þ; the search direction (13) for the pairs ðF E ; W 0
E E 0 the interface equilibrium: 8M 2 CEE ; b F n ðMÞ þ b F n ðMÞ ¼ 0;
the (usually nonlinear) evolution laws, which, if one defines the 0
velocity jump by EE0
8M 2 C ;
_ E _ E c c ¼W n W n , may take the form:
8 9 E
: Y EE ; d;n 0
( ¼R
_ EE0 ðMÞ W n 0 d_ EE
3.5. About discretization A standard finite element discretization is used for the displacement field in the substructures. So the displacement U and the strain e fields have the following form:
8
E > _ E > _ E W c > ð15Þ F Enþ1 b F n ¼ k0 W > nþ1 n > < 8 9 ( ) E snþ1 bs n 2 E () E < b e p;n = > > e_ p;nþ1 ¼ > ð16Þ > > :X : XE bE ; nþ1
_ EE0 W nþ1=2
or linear. The constitutive relation, which depends on internal variables (damage), which are a priori unknown, is nonlinear. The search direction (16), which is part of the constitutive relation (11), enables the calculations to be carried out for a damage value, which is fixed and known (because it was obtained in the local step). Thus, once the damage is known, the constitutive relation becomes a linear relation. Once this difficulty has been eliminated, the separation of the equations is a simple task: the solutions of the kinematic admissibility, equilibrium and constitutive relations (10, 11) belong to Ad , while the equations of the evolution laws (12) belong to C, as shown in Table 3.
The resolution of the linear step considers that the complete state of the interfaces (defined by forces fields, velocity fields and eventually damage variables) is computed during the last local step. The parameter k0 , defined in the Eq. (15), can be considered such as a stiffness between the interface and the substructures. Consequently, a substructure is subjected to mixed boundary conditions. Moreover, during this last local step, the damage state of each substructure is also computed. This leads to a material behavior, which can be considered elastic and linear because the damage state is given. The linear substructure problem takes into account the kinematic admissibility, the linear hebavior at prescribed damage state and the search direction. It becomes, after discretization:
b_ E ðtÞg _ k0 ½HE fuðtÞg þ ½K E fuðtÞg ¼ ½HE bf En1 ðtÞ þ k0 f w n where:
)! ð17Þ
n
The search directions at the interface (13,15) introduce a parameter k0 whose optimal value with respect to the convergence rate of the method was studied in Ladevèze et al. (2001). 3.4. Determination of the substructure problem The mechanical problem to be solved in each substructure XE consists in finding for each time step ðrE ; eEp ; XE YE Þ the evolution, which verifies: the kinematic admissibility; the equilibrium; the damageable behavior of the mesomodel’s plies (Eqs. (10)–(12)) A major difficulty related to the introduction of the nonlinear behavior of the substructures resides in the handling of the constitutive relation (10) in the LATIN algorithm. Indeed, the method is based on the separation of the equations, which are either local
½HE ¼
Z @ XE
½Nt ½NdSE and ½K E ¼
Z XE
½Bt Kd ½BdXE
½K E is the classical finite element stiffness matrix of substructure XE at given damage state and ½HE is the boundary term related to the interfaces.
Table 3 Separation of the substructure’s constitutive relations. Step
Equations to be solved
Linear
In the substructures: Kinematic admissibility Equilibrium r ¼ Kd ee Y ¼ AX Search directions (15, 16)
Local
On the substructures: " # e_ p r ¼B _ Y X At the interfaces: Behavior Search directions (13, 14)
2754
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
For the resolution of this differential equation, an Euler implicit time integration scheme is used. During the local step, the thermodynamical forces associated to the damage state of the cohesive interfaces (defining by the resolution of the Eqs. (1)–(3)) and of the plies (defining by the resolution of the Eqs. (4)–(8)) are given. The use of an Euler implicite time integration scheme leads to a non-linear problem, defined at each Gauss point, which is solved with a Newton–Raphson scheme. 3.6. The multiparametric strategy 3.6.1. Principle The objective of the multiparametric strategy is to solve a large number of mechanically ‘‘similar’’ problems, each associated with a set of parameters. The approach followed in this paper relies on the possibility to initialize the LATIN method using any solution, usually an elastic solution. We will take advantage of this possibility by initializing the resolution of a problem using the result of a previously solved problem whose material parameters are considered to be relatively similar. Thus, a first, mechanically relevant approximation of the solution of the problem is available at no additional computation cost. A graphical representation of this strategy is given by Fig. 5. Following the calculation of a solution associated with problem AAd \ CA , a new set of parameters is considered. In order to simplify the figure, since the variabilities usually concern the nonlinear equations, let us assume that AAd ABd . Thus, the difference between problem A and problem B is contained in the solutions of the local equations. Assuming that the parameters are only slightly different, the space defined by CB is relatively close to that defined by CA . Therefore, to initialize the resolution of problem B, it is usually more relevant to use the solution of problem A than to start from an elastic solution. 3.6.2. Sources of variabilities From a pragmatic standpoint, for designing an assembly of parts made of composite materials, the parameters of the problem can be categorized according to their range of variation: Elastic parameters: the elastic moduli and Poisson coefficients are usually well-known with few identification tests. Therefore, we will consider that a parametric study of their influence is irrelevant. Nonlinear interface parameters: even though the friction coefficients and the parameters involved in the cohesive laws often have a major influence on the solution, they are difficult to identify. Therefore, a parametric study of these parameters is particularly appropriate. Nonlinear substructure parameters: the influence of these parameters on the stiffness matrix is significant, but often very local, which makes them difficult to identify. Therefore, it may be very relevant to resort to a parametric study to identify them. We performed two parametric studies related to the last two groups of parameters. The first example concerns the case in which the interface parameters vary; the second example addresses the case in which the parameters affect the laws governing the substructures. 4. Application to a cohesive interface problem 4.1. Definition of the problem This first example assumes that the equations of the substructures are deterministic and that the variabilities concern the inter-
Fig. 5. Graphical interpretation of the multiparametric method.
faces alone. The objective of this section is to apply the multiparametric strategy to a simple case including the cohesive interfaces of the standard damage mesomodel. The example concerns the opening of a laminate subjected to a normal compression load. The laminate consisted of 5 tridimensional plies, oriented in the 0° and 90° directions alternatively, along with titanium blocks whose behavior was assumed to be linear, elastic and isotropic (see Fig. 6). 0 The multiparametric study was performed on parameters kn and Y 0 , which correspond to the initial stiffness of the interface and to the damage threshold and were assumed to be equal at each interface. Indeed, during a identification test, the value of these parameters are difficult to determinate and cannot be well known. The range of variation of each parameter was discretized into 26 arbitrarily values, leading to 676 different problems (see Table 4).
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
2755
Fig. 6. A view of the laminate – the behavior of a cohesive interface.
The efficiency of the strategy is expressed by the gain G defined by:
Table 4 0 The values of parameters Y C and kn . Parameter
Min. value
Max. value
Step
YC
0.1 MPa 103 N/mm
1.5 MPa 107 N/mm
0.056 MPa
0
kn a
G¼
s625 625 s1 ’ l625 l625
a
For all these non-linear analysis, the gain is computed considering that:
Denotes a constant step in the logarithmic space.
for the resolution of the first problem, we use the LATIN method to solve the non-linear problem (associated to CPU time s1 ), instead of an usual Finite Element software, the LATIN method is used to obtain the CPU time reference s1 ; however they are equivalent as showed in Roulet et al. (2011) with a comparison between the commercial software ABAQUS and the LATIN method, the CPU time associated to the 625 resolutions is estimated, assuming that the CPU time is constant for each non-linear analysis.
Fig. 7. Influence of the parameters on the opening of the laminate.
Table 5 CPU time and gain achieved with the multiparametric strategy. Number of processors
s1
l625
Gain
1 6
32 min 10 min
1098 min 330 min
18 19 0
Let us point out that on the figure the variation range of kn is expressed in logarithmic scale, and that it is on that scale that the points were regularly spaced.
The gain shown in Table 5 is about 18 for both parallel and sequential computations. Thus, despite its ease of implementation, the multiparametric strategy leads to a significant improvement in the computation time required for a parametric study. In this example, the nonlinearities were concentrated only at the interfaces. The substructures were considered to have linear elastic behavior. Therefore, the state of the complete structure was defined by the values of the fields at the interfaces, and the values of the fields within the substructures could be obtained directly through an elastic calculation with prescribed boundaries. Thus, only the interface data needed to be stored for each problem and transferred from one problem to another. The introduction of nonlinearities on the substructure level precludes the use of this trick and raises the crucial question of the storage of the solutions, which is addressed by the second example.
4.2. Results 5. Application to damageable plies The parameters of the cohesive laws play an important role in the opening of each ply. Indeed, Fig. 7 shows the opening of the laminate (defined as the difference in the end displacements of the upper and lower plies) as a function of the parameters. Let us designate by:
s1 , the CPU time required for the resolution of the first problem; s625 , the CPU time required for the 625 resolutions, here estimated to be 625 s1 ; l625 , the CPU time required for the 625 resolutions using the
multiparametric strategy.
For this example, the parameters involved in the evolution of the damage were allowed to vary during the parametric study. As explained above, the introduction of nonlinearities in the substructures leads to the question of the storage of the solutions. Indeed, while in the elastic case the fields at the interfaces were sufficient, now this is no longer the case. Nevertheless, we will try to perform a calculation initialized by the previously calculated interface fields alone. This calculation will be compared to a calculation whose initialization took into account the damage states of all the substructures.
2756
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
5.1. Description of the problem In the problem considered here, the damage state of the material influences the global behavior of the assembly drastically. The objective was to carry out a 3D simulation of a 4-point bending test (see Fig. 8). The three plates were made of a damageable laminate and assembled by elastic bolts. The loading was applied in two steps: a preloading step: a prescribed relative displacement U i between the head and the thread of each bolt; a bending step: a prescribed normal displacement at the head of each bolt to induce the bending of the assembly.
5.2. Parametric study
Table 6 Variations of the parameters. Parameter
Min. value
Max. value
step
Ui
0.01 mm 0.1 MPa
0.1 mm 1 MPa
0.01 mm 0.1 MPa
YC
Table 7 The memory storage requirements for the solution of one problem. Initialization from . . .
Memory requirement
the interface fields alone the interface fields and the internal variables of the substructures Stiffness matrices (for one time step)
1069.7 kB 2661.4 kB 32,944.6 kB
A parametric study was carried out for two parameters: the preload in the bolts U i and the critical value Y C . These parameters have a non-negligible influence on the solution and are often unknown. Each parameter was allowed 11 arbitrarily different values (see Table 6), leading to 121 problems. With regard to the richness of the initialization, which conditions the storage space required for each solution, we considered two approaches: an initialization from the interface fields alone, thus assuming no initial damage; an initialization from the interface fields and the damage state of the substructures. Table 7 gives an idea of the memory required to store a solution corresponding to the mesh used for this example. For each Gauss point of each interface element, 7 floating-point numbers are necessary (3 for the displacements, 3 for the loads, and one for the damage variable); for each Gauss point in the substructures, 3 0 numbers are necessary (d; d and dF ). In other words, for the discretization considered here (linear elements, 5 elements through the thickness of each ply), taking into account the damage state in each Gauss point of each substructure amounts to doubling the storage space required. While storing the matrices is not necessary in order to use the multiparametric strategy, this enables one to compare the amount of RAM required for the calculation and the disk space required for storing the results.
5.3. Results The relevant quantity of interest for an easy comparison of the solutions is necessarily a function of the damage state of the structure. Therefore, we chose to represent, in Fig. 9, the reaction forces over the surface where the normal displacements were prescribed. The study of the relevance of the two types of initialization is summarized in Table 8. One can see that for the multiparametric strategy to be relevant one must reinstate the damage state of the previous calculation right from the initialization. The informa-
Fig. 9. Influence of the parameters on the reaction force.
Table 8 The gain observed in the parametric study. Initialization
1 Resolution
121 Resolutions
Gain
Interfaces Interfaces & SST
810 s 810 s
95,590 s 7159 s
1.02 13.7
tion on the interface level is insufficiently rich, so the gain is practically inexistent. 6. Conclusion We presented an efficient computational strategy for the resolution of multiple nonlinear problems, especially when these problems involve damage. This strategy relies on a domain decomposition, which makes parallelization easy, and on the LATIN iterative solver. A relevant initialization of the calculation by reusing previously computed solutions leads to a significant reduction in the number of LATIN iterations, which are necessary, and, therefore, in the computation time. In the case of variabilities in the interface laws, only the interface fields are stored. In the case of nonlinear substructures with variable behavior, storing the state of the internal variables is necessary in order to preserve significant gains in computation time. The two examples given illustrate the efficiency of the multiparametric strategy as long as these initialization conditions are respected. Acknowledgement
Fig. 8. The geometry of the 4-point bending problem.
The research leading to these results was funded by the European Community’s Seventh Framework Program FP7/2007–2013 under Grant Agreement No. 213371 (http://www.maaximus.eu/).
V. Roulet et al. / International Journal of Solids and Structures 50 (2013) 2749–2757
References Alart, P., Curnier, A., 1991. A mixed formulation for frictional contact problems prone to newton like solution methods. Computer Methods in Applied Mechanics and Engineering 92 (3), 353–375. Allix, O., 1992. Damage analysis of delamination around a hole. New Advances in Computational Structural Mechanics, pp. 411–421. Allix, O., Ladevèze, P., Corigliano, A., 1995. Damage analysis of interlaminar fracture specimens. Composite Structures 31 (1), 61–74. Armero, F., Petcz, E., 1998. Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems. Computer Methods in Applied Mechanics and Engineering 158 (3–4), 269–300. Arora, J.S., Chahande, A.I., Paeng, J.K., 1991. Multiplier methods for engineering optimization. International Journal for Numerical Methods in Engineering 32 (7), 1485–1525. Azzi, V.D., Tsai, S.W., 1965. Anisotropic strength of composites. Experimental Mechanics 5 (9), 283–288. Berthelot, J.-M., 2003. Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: static and fatigue loading. Applied Mechanics Reviews 56 (1), 111–147. Bordeu, F., Boucard, P., Lubineau, G., 2011. A mesoscale model for damage, cracking and delamination prediction in composite materials. Composite Materials 17 (4), 271–282. Boucard, P.A., Champaney, L., 2003. A suitable computational strategy for the parametric analysis of problems with multiple contact. International Journal for Numerical Methods in Engineering 57 (9), 1259–1281. Chabrand, P., Dubois, F., Raous, M., 1998. Various numerical methods for solving unilateral contact problems with friction. Mathematical and Computer Modelling 28 (4–8), 97–108, recent Advances in Contact Mechanics. Champaney, L., 2011. A domain decomposition method for studying the effects of missing fasteners on the behavior of structural assemblies with contact and friction. Computer Methods in Applied Mechanics and Engineering (205/208), 121–129. http://dx.doi.org/10.1016/j.cma.2011.04.008. Champaney, L., Boucard, P.A., Guinard, S., 2008. Adaptive multi-analysis strategy for contact problems with friction. Computational Mechanics 42 (2), 305–315. Champaney, L., Cognard, J.Y., Ladevèze, P., 1999. Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions. Computers & Structures 73 (1–5), 249–266. Dostal, Z., Friedlander, A., Santos, S.A., 1998. Solution of coercive and semicoercive contact problems by feti domain decomposition. Contemporary Mathematics 218, 82–93. Dureisseix, D., Farhat, C., 2001. A numerically scalable domain decomposition method for the solution of frictionless contact problems. International Journal for Numerical Methods in Engineering 50 (12), 2643–2666. Farhat, C., Lesoinne, M., LeTallec, P., Pierson, K., Rixen, D., 2001. Feti-dp: a dualprimal unified feti method-part i: a faster alternative to the two-level feti method. International Journal for Numerical Methods in Engineering 50 (7), 1523–1544. Farhat, C., Lesoinne, M., Pierson, K., 2000. A scalable dual-primal domain decomposition method. Numerical linear algebra with applications 7 (7–8), 687–714. Farhat, C., Roux, F.X., 1992. An unconventional domain decomposition method for an efficient parallel solution of large-scale finite element systems. SIAM Journal on Scientific and Statistical Computing 13, 379. Feyel, F., 2003. A multilevel finite element method (FE2) to describe the response of highly non-linear structures using generalized continua. Computer Methods in Applied Mechanics and Engineering 192 (28–30), 3233–3244. Feyel, F., Chaboche, J.-L., 2000. Fe2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials. Computer Methods in Applied Mechanics and Engineering 183 (3–4), 309–330. Ghanem, R.G., 1999. Ingredients for a general purpose stochastic finite elements implementation. Computer Methods in Applied Mechanics and Engineering 168, 19–34. Ghanem, R.G., Kruger, R.M., 1996. Numerical solution of spectral stochastic finite element systems. Computer Methods in Applied Mechanics and Engineering 129 (3), 289–303. Ghanem, R.G., Spanos, P.D., 1991. Stochastic Finite Elements: A Spectral Approach. Springer, Berlin. Ghosh, S., Lee, K., Moorthy, S., 1995. Multiple scale analysis of heterogeneous elastic structures using homogenization theory and voronoi cell finite element method. International Journal of Solids and Structures 32 (1), 27–62. Ghosh, S., Lee, K., Raghavan, P., 2001. A multi-level computational model for multiscale damage analysis in composite and porous materials. International Journal of Solids and Structures 38 (14), 2335–2385.
2757
Hinton, M.J., Kaddour, A.S., Soden, P.D., 2004. Failure criteria in fibre reinforced polymer composites: the world-wide failure exercise. Elsevier Science. Kachanov, L.M., 1958. On creep rupture time. Izv. Acad. Nauk SSSR, Otd. Techn. Nauk 8, 26–31. Kaminski, M., 2005. Computational Mechanics of Composite Materials: Sensitivity, Randomness and Multiscale Behaviour. Springer. Kerfriden, P., Allix, O., Gosselet, P., 2009. A three-scale domain decomposition method for the 3D analysis of debonding in laminates. Computational Mechanics 44 (3), 343–362. Kikuchi, N., 1982. Penalty/finite element approximations of a class of unilateral contact problems. Penalty method and finite element method. ASME, New York. Kikuchi, N., Oden, J.T., 1988. Contact problems in elasticity: a study of variational inequalities and finite element methods. Society for Industrial Mathematics. Klarbring, A., 1992. Mathematical programming and augmented Lagrangian methods for frictional contact problems. In: Proceedings of Contact Mechanics International Symposium, PPUR, vol. 369, p. 390. Kleiber, M., Hien, T.D., 1992. The Stochastic Finite Element Method. John Wiley & Sons, New York. Ladevèze, P., 1985. Sur une famille d’algorithmes en mécanique des structures. Compte rendu de l’académie des Sciences 300 (2), 41–44. Ladevèze, P., 1986. Sur la mécanique de l’endommagement des composites. C. Bathias et D. Menkes, éds: Comptes-rendus des JNC5, Paris. Ladevèze, P., 1999. Nonlinear Computational Structural Mechanics: New Approaches and Non-incremental Methods of Calculation. Springer Verlag. Ladevèze, P., Allix, O., Gornet, L., Lévèque, D., Perret, L., 1998. A computational damage mechanics approach for laminates: identification and comparison with experimental results. In: George Z. Voyiadjis, J.-W.W.J., Chaboche, J.-L. (Eds.), Damage Mechanics in Engineering Materials, vol. 46 of Studies in Applied Mechanics. Elsevier, pp. 481–500. Ladevèze, P., Dureisseix, D., 2000. A micro/macro approach for parallel computing of heterogeneous structures. International Journal for Computational Civil and Structural Engineering 1 (1), 18–28. Ladevèze, P., Le Dantec, E., 1992. Damage modelling of the elementary ply for laminated composites. Composites Science and Technology 43 (3), 257–267. Ladevèze, P., Loiseau, O., Dureisseix, D., 2001. A micro-macro and parallel computational strategy for highly heterogeneous structures. International Journal for Numerical Methods in Engineering 52 (1–2), 121–138. Le Tallec, P., Roeck, Y.H., Vidrascu, M., 1991. Domain decomposition methods for large linearly elliptic three-dimensional problems. Journal of Computational and Applied Mathematics 34 (1), 93–117. Lubineau, G., Ladevèze, P., 2008. Construction of a micromechanics-based intralaminar mesomodel, and illustrations in abaqus/standard. Computational Materials Science 43 (1), 137–145. Mandel, J., 1993. Balancing domain decomposition. Communications in Numerical Methods in Engineering 9 (3), 233–241. Mandel, J., Tezaur, R., 1996. Convergence of a substructuring method with lagrange multipliers. Numerische Mathematik 73 (4), 473–487. Nairn, J.A., 2000. Matrix microcracking in composites. Polymer Matrix Composites 2, 403–432. Nairn, J.A., Hu, S., 1994. Matrix microcracking. Composite Materials Series, pp. 187– 243. Radi, B., Baba, O.A., Gelin, J.C., 1998. Treatment of the frictional contact via a Lagrangian formulation. Mathematical and Computer Modelling 28 (4–8), 407– 412, recent Advances in Contact Mechanics. Rixen, D.J., Farhat, C., 1999. A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems. International Journal for Numerical Methods in Engineering 44 (4), 489–516. Roulet, V., Champaney, L., Boucard, P.A., 2011. A parallel strategy for the multiparametric analysis of structures with large contact and friction surfaces. Advances in Engineering Software 42 (6), 347–358. Series, L., Feyel, F., Roux, F.X., 2003. Une méthode de décomposition de domaine avec deux multiplicateurs de lagrange. Simo, J.C., Laursen, T.A., 1992. An augmented lagrangian treatment of contact problems involving friction. Computers & Structures 42 (1), 97–116. Tsai, S.W., Wu, E.M., 1971. A general theory of strength for anisotropic materials. Journal of Composite Materials 5 (1), 58. Wriggers, P., 1995. Finite element algorithms for contact problems. Archives of Computational Methods in Engineering 2 (4), 1–49. Wriggers, P., Simo, J.C., Taylor, R.L., 1985. Penalty and augmented lagrangian formulations for contact problems 85, 97–106. Zhong, Z., Mackerle, J., 1992. Static contact problems – a review. Engineering Computations 9 (3–37).