- Email: [email protected]
Iterative solution methods for damage progression analysis
Composite Structures 69 (2005) 21–33 www.elsevier.com/locate/compstruct
Iterative solution methods for damage progression analysis G. Kress *, M. Siau, P. Ermanni Structure Technologies, ETH Z€urich, LEO C2, Leonhardstr. 27, CH-8092 Z€urich, Switzerland Available online 2 June 2004
Abstract Damage progression is a sequence of damage events. A simulation of it solves a sequence of systems of equations which differ from each other only by the effects of the modelling of the incrementally increasing damage. An initial matrix triangulation, to solve the equations of the virgin structure, provides a preconditioner that can be used for solving the sequence of system equations by approximate or iterative solution techniques. In addition, obtained solutions serve as starting vectors for the next solutions. The numerical efficiency potential of this concept is evaluated by considering a method inspired by the sensitivity formula for finite element models, an iterative method inspired by Jacobi’s method, and the untransformed preconditioned conjugate gradient method. The ability to reproduce the damage patterns obtained by using exact solution is checked and the efficiency improvements are assessed. 2004 Elsevier Ltd. All rights reserved. Keywords: FEM; Numerical methods; Damage simulation; Composites
1. Introduction Typically, iterative methods for solving linear systems of equations can outperform direct solution methods when the number of unknowns is very high and the coefficient matrix is sparsely populated. Iterative solution methods may be competitive even for smaller problems with a banded coefficient matrix when the sought solution is closely guessed by a starting vector and the coefficient matrix is well conditioned. These circumstances can be provided at reasonable cost when variations of the same problem must be evaluated sequentially. Examples include the sequential refinement of design solutions towards some specified optimum or the progression of damage in composite materials in small steps. 1.1. Damage in composite materials Fig. 1 shows a state of damage in a notched test specimen made from a laminate with the stacking sequence [45,90,)45,0]s made from a T300/5208 graphite/ epoxy prepreg tape. The specimen suffered 104 cycles at a maximum stress of 85% of their mean quasi-static tensile strength. The radiograph shows that matrix *
Corresponding author. Fax: +41-1-633-1125. E-mail address: [email protected] (G. Kress).
0263-8223/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2004.04.043
cracking in all layers as well as delaminations exist. The delaminations appear in the radiograph in darker shades than the undamaged parts of the laminate. Characteristic features are the matrix cracks in the longitudinal direction tangent to the hole in the 0-degree plies, the matrix cracks emanating from the hole and extending almost to the straight edges of the specimen, and the delaminations seeming to be companions of the tangent matrix cracks. The fact that the diagonal matrix cracks extend much farther in the 45- than in the )45-direction has to do with the stacking sequence placing the 45degree layer at the surface. Please note that this is an effect that could only be captured by a fully threedimensional model, which also holds true for the delaminations. Obviously, it is not an easy task to simulate the damage development under fatigue loading. For ease, a very simple damage model is considered sufficient for the efficiency study although more sophisticated models exist. 1.2. Overview of the present study The present study chooses composite material damage progression analysis just as one of the various possible backdrops to verify approaches for taking advantage of iterative solution techniques for simulations involving a sequence of structural models which
22
G. Kress et al. / Composite Structures 69 (2005) 21–33
accuracy, and computing time. With the other laminates more complex damage patterns are simulated in order to compare the results of the direct and iterative solution methods. The conclusions of this study are found in Section 6.
2. Progressive damage analysis 2.1. Mechanical and residual stresses
Fig. 1. Damage in notched [45,90,)45,0]s -specimen [1].
differ from each other by small variations. Section 2 presents a very simply damage analysis approach conceptually including the effects of residual stresses due to temperature or moisture loads. Failure criteria are taken from [2] but the effect of shear stress on tensile fiber failure is omitted to achieve damage patterns resembling more closely those due to fatigue. The damage is modelled by drastically reducing stiffness parameters in directions perpendicular to the found failure lines. Simple as it is, the modelling has led to fair predictions of quasistatic strength of specimens with severe edge effects [3]. The alteration of the numerical system, initially representing the virgin structure, by damage modelling is considered in Section 3. From this, convenient formulations of the change-of-displacement field are presented and preconditioning with hindsight to iterative solution techniques is explained. Specifically, an iteration inspired by the Jacobi Method and a suited version of the Conjugate Gradient Method are considered. Attention is also given to a method derived from the well-known sensitivity formula for gradient calculation of the total potential energy. It could be used to hurry most quickly through the steps of an iterative damage progression. The sample problem structure is a notched test specimen and Section 4 explains the FEM modelling. It maps only a plane-stress state and assumes symmetric damage states which is clearly controversial as Fig. 1 testifies. The laminate stacking sequences and material properties used in the simulations are also given in this section. Numerical results are presented in Section 5. The simplest damage propagation pattern, here dubbed Unzip Action, is reproduced by all considered solution methods. It is used to investigate aspects of convergence,
The unidirectionally reinforced composite is modelled as a transversely isotropic homogeneous material with different properties in the fiber direction and those transverse to it. When temperature changes, the free thermal expansion of an individual layer is constrained by the other layers of the multidirectional laminate. During the autoclave process, the temperature of the laminate increases but the viscosity of the matrix material prevents temperature stresses. When the matrix material cures, internal residual stresses rres are initially negligible but increase with falling temperature. When the difference DT to the stress-free reference temperature becomes very large, laminates made from composites with brittle matrices may experience matrix cracking even before any external load is applied. In order to find the mechanical load causing damage, the designer must therefore take into account the presence of residual stresses before application of any boundary conditions. We assume a situation where the operating temperature and the mechanical load case are specified and the designer must answer the question at what multiple of the mechanical loads damage or failure of the part will occur. This leads to the partition of the stresses at any point in the laminated structure r ¼ krmec þ rres ;
ð1Þ
where the load factor k is obtained from substituting (1) into suitable failure criteria so that f ðrÞ ¼ 1:
ð2Þ
The theory of laminated plates describes the relationship between combined mechanical and thermal loading and the plate deformations with
N M
A ¼ B
B C
0 j
by using so-called fictitious line loads T N N N ¼ þ : M M MT
ð3Þ
ð4Þ
The stresses rres follow from a FEM analysis where the temperature- and moisture caused line loads NT are calculated on element level and assembled to the righthand-side vector. The essential boundary conditions are
G. Kress et al. / Composite Structures 69 (2005) 21–33
defined to give a statically determinate support to allow for the free deformation of the structure. The resulting deformation along the boundary Cu with fixed displacements must be taken into account when implementing the homogeneous geometric boundary conditions. They must not change the stress state of the free structure caused by the temperature- and moisture loads. Then, the stresses rmec follow from a FEM analysis where the geometric boundary conditions and only the mechanical external forces are applied. 2.2. Failure criteria As a matrix crack and a fiber rupture will differently change the load-carrying behavior of the laminate, a suitable failure criterion must be able to predict the failure mode. It should also be able to take into account the effect of combined stressing on the failure probability. The present study uses the four criteria by [2] for fiber and matrix failure under tensile and compressive loading under plane-stress loading. The strength parameters are obtained from single-stress tests: tensile strength in fiber direction compressive strength in fiber direction tensile strength in transverse direction compressive strength in transverse direction shear strength longitudinal shear strength transverse
XT XC YT YC SL ST
r1 XT
2
þ
s21 SL
ð5Þ
Fiber failure under compressive fiber stress (r1 < 0): r1 þ XC ¼ 0: Matrix failure under tensile stress (r2 > 0): 2 2 r2 s12 þ ¼ 1: YT SL Matrix failure under compressive stress (r2 < 0): # 2 " 2 2 r2 YC r2 s12 þ ¼ 1: 1 2ST 2ST YC SL
Critical factors on the mechanical loading are obtained from combining (1) and (2) with a suitable criterion. Quadratic failure criteria then lead to quadratic equations for the critical load k2 Q þ kL þ P ¼ 1:
ð9Þ
A value of P greater than 1 indicates that residual stresses be high enough to cause failure. Then, the damage resulting from this failure must be modelled but the factor on the mechanical load case is zero. Else, k is calculated by k¼
L þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi L2 þ 4Qð1 P Þ 2Q
ð10Þ
or, if the criterion is actually linear, Q ¼ 0, by k¼
1P : L
ð11Þ
Failure is recognized element-wise. The intralaminar damage modes of fiber fracture and matrix failure are crudely modelled by reducing the respective Young’s moduli in the respective layers to a very small fraction of their original values, responding to the lost capability of transferring load across the fractured surfaces:
2 ¼ 1:
2.3. Prediction of critical load factors
2.4. Damage modelling
Fiber failure under tensile fiber stress (r1 > 0):
23
ð6Þ
G12 ¼ G12 ;
m12 ¼ m21 ¼ 0;
c: fiber:
G12 ¼ G12 ;
m12 ¼ m21 ¼ 0;
t: matrix: E2 ¼ E1 ; c: matrix:
G12 ¼ G12 ; G12 ¼ G12 ;
m12 ¼ m21 ¼ 0; m12 ¼ m21 ¼ 0:
t: fiber:
E1 ¼ E1 ;
ð12Þ 2.5. Termination criteria
ð7Þ
ð8Þ
The simulations presented in Section 5 are based on one rather arbitrary alteration of Hashin’s criteria: the shear-stress term is removed from the criterion for fiber fracture under tensile stress (5). This leads to damage patterns more similar to those appearing under cyclic loading, see Fig. 1, where extensive damage of the matrix occurs before fibers fail.
The damage progression analysis is a repeatedly applied first-ply-failure analysis. After a failure event is detected, the structural model is updated to account for the damage and, due to the altered stress state, a new failure event will be detected at a different location. The iterations stop upon fulfillment of one of the following criteria: • the number of damage events reaches a preset value, • the work of the external forces exceeds a preset value. The latter criterion can also be used to recognize total failure where the structure ceases to be one part.
24
G. Kress et al. / Composite Structures 69 (2005) 21–33 1
3. Damage states and displacement fields
ðI þ Kð0Þ DKÞD~u ¼ D~u þ LD~u;
Each damage event reduces the apparent material stiffness of the region surrounding it. Considering a progressive failure analysis, where the behavior of the structural model is simulated by a finite-element model,
where the operation LD~u first assembles the vector DKD~u by matrix multiplications on element level and then applies the back-substitution process on it.
K~u ¼ r;
3.1. Functional and gradient
ð13Þ
the region surrounding a new damage event is presented by one finite element. Its stiffness matrix will change, representing the loss of material stiffness, and possibly internal forces will be set free in the presence of residual stresses. This will also change the displacements of the structure and with it the local material deformations so that each damage influences the spatial probability distribution of subsequent failure events. All of these changes are captured in the following equation: ðKð0Þ þ DKÞð~ uð0Þ þ D~ uÞ ¼ ðrð0Þ þ DrÞ:
ð14Þ
ð18Þ
Some iterative methods obtain the displacement solution of (16) by minimizing the quadratic functional f ðD~uÞ i h f ¼ D~uT 12ðKð0Þ þ DKÞD~u ðDr DK~uð0Þ Þ : ð19Þ This is the same as seeking the zero-point of its gradient rf ¼ Kð0Þ D~u þ DKð~uð0Þ þ D~uÞ Dr:
ð20Þ
3.2. Jacobi method inspired iterative solution method
ð0Þ
refers the respective entities to the The upper index undamaged state of the structure. Removing the system equations referring to the initial, or undamaged, state of the simulated structure
The Jacobi Method transforms a given system of Eq. (13), by partitioning the original coefficient matrix K into the diagonal matrix D and the remainder E,
Kð0Þ ~ uð0Þ ¼ rð0Þ
K~u ¼ ~r;
ð15Þ
~u ¼ D E~u þ D1~r;
reduces (14) to the form ðKð0Þ þ DKÞD~ u ¼ Dr DK~ uð0Þ :
ð16Þ
To evaluate (16), let it be that the system (15) has already been solved to obtain the displacement field ~uð0Þ . This implies that the system matrix Kð0Þ can be inverted because the appropriate essential boundary conditions are implemented. While (14) obtains the total displacements corresponding to a given instant of the damage history, (16) obtains the changes from the initial displacements found by (15). The changes D~ u due to the consecutive failure events incrementally modify the initial solution ~ uð0Þ obtained by (15). If the differences between two consecutive displacement fields are small, the numerical effort of iterative solution techniques may become less than that of direct ones because the previous displacements provide a good starting vector for the new ones. Furthermore, the triangulized form of Kð0Þ obtained from the direct solution of (15) can be used as a preconditioner for improving the convergence of iterative methods for approximating the solution of (16). Applying Kð0Þ to (16) shows1 that only the damaged elements cause deviation Kð0Þ DK from the ideal unity matrix: 1
1
ðI þ Kð0Þ DKÞD~ u ¼ Kð0Þ ðDr DK~ uð0Þ Þ:
D~u ¼ E~u þ ~r; 1
ð21Þ
so that an iteration rule results: ~uiþ1 ¼ B~ui þ ~z:
ð22Þ
The iterations converge to the exact solution because the convergence radius of B is smaller than unity or qðBÞ ¼ max jkj j < 1 due to the fact that the values of the entries of the main diagonal D of K are larger than those of the off-diagonal entries. The operation B~uj calculates the entries of the new vector ~ujþ1 sequentially, starting from its first entry ð~uðjþ1Þ Þ1 . The convergence is enhanced by replacing the entries of the previous vector ~uj with the already calculated entries of ~ujþ1 . The improving feature distinguishes the Gauss-Seidel Method from the Jacobi Method. Re-ordering the terms in (17) provides a form similar to that in the last line of (21) 1
1
D~u ¼ Kð0Þ DKD~u þ Kð0Þ ðDr DK~uð0Þ Þ; 1
ð23Þ
1
where Kð0Þ DK and Kð0Þ ðDr DK~uð0Þ Þ have similar meanings as B and z, respectively. The equation can be solved in iterations analogous to (22): 1
D~uiþ1 ¼ Kð0Þ Dr Lð~uð0Þ þ D~ui Þ
ð24Þ
or, if D~u0 ¼ 0,
ð17Þ 1
The deviation, however, is 1 generally not symmetric even though both matrices, Kð0Þ and DK, are symmetric. However, the left-hand side of (17) can also be evaluated as
~u ¼ B~u þ ~z
D~uiþ1 ¼ Kð0Þ Dr þ
i X
i
ð1Þ Li ~u0 :
ð25Þ
1
The sum of the displacements u~ð0Þ of the undamaged structure and the changes D~uðpÞ due to previous damage
G. Kress et al. / Composite Structures 69 (2005) 21–33
events provides a better starting vector. Then, the iterations yield 1
D~uiþ1 ¼ Kð0Þ Dr þ
i X i ð1Þ Li ~ uðpÞ : u0 Li D~
ð26Þ
1
The series converges linearly if it holds for any vector v that i
jL vj < jvj:
ð27Þ
The method presented here is not identical with the Jacobi Method: the preconditioner is the inverse of the complete initial coefficient matrix, and it is applied to the varied stiffness matrices of the damage progression simulation sequence. Therefore, it is here referred to as Jacobi-inspired method (JIM). The iterations are continued with i ¼ i þ 1 until a maximum number is reached or the termination criterion indicates that the chosen accuracy measure is fulfilled. 3.3. Method derived from sensitivity formula
1
solves the problem of using the preconditioner Kð0Þ to improve convergence to the solution of (16) and (17). The gradient rf0 at the starting point of the search is calculated directly from (20). The initial search direction follows from the operation 1
s0 ¼ Kð0Þ rf0
ð30Þ
performing back substitution with the existing triangulized stiffness matrix of the undamaged structure. The factor qnew is the inner product of gradient and search direction qnew ¼ rf0T s0
ð31Þ
and it is needed for both step-length and conjugate search direction calculations. The initial value is stored q0 ¼ qnew
ð32Þ 2
for use in the termination criterion qnew < q0 stopping the following iterative part of the algorithm. A vector q follows from applying the stiffness matrix K to the current search direction: qi ¼ Ksi :
If the damage modelling leads only to small stiffness changes, the changes of displacements are also small compared with the displacements of the undamaged structure. This may be the case for structures made from multidirectional laminates consisting of fiber-reinforced materials whose limit strain in the transverse direction is much smaller than in the fiber direction and when only matrix failure occurs. Then, the nonlinear term, requiring iterations, may be dropped so that the linear form results: 1
D~u ¼ Kð0Þ ðDr DK~ uð0Þ Þ:
ð28Þ
This can also be derived from the formula used to evaluate the sensitivity of the displacement field on some structural parameters x: 1 o~u or oK ð0Þ ~ ¼ Kð0Þ ; ð29Þ u oxi oxi oxi where the xi are structural stiffness parameters. From the reference point, i.e. the displacement field of the undamaged structure, a step along the direction of the gradient must be taken the length of which is given by the amount of stiffness change of the element where damage is modelled. This leads again to (28). Obviously, the method (SFM) gives the same result as the first iteration step of the method presented in Section 3.2.
25
ð33Þ
The step length q ai ¼ Tnew s i qi
ð34Þ
leads from the current reference point D~ui1 to the new point D~ui minimizing f (19) along the search direction si : D~ui ¼ D~uði1Þ þ ai si :
ð35Þ
The gradient is updated by: rfi ¼ rfi1 þ ai qi :
ð36Þ
A vector pi is formed by applying back substitution on the gradient rfi pi ¼ Kð0Þ1 rf ðiÞ :
ð37Þ
The factor q is updated by qold ¼ qnew ;
qnew ¼ rfiT pi :
The factor bðiÞ , q bðiÞ ¼ new qold
ð38Þ
ð39Þ
renders consecutive search directions si ¼ pi þ bsi1 ;
ð40Þ
K-conjugate.
4. Sample problem 3.4. Conjugate gradient method The preconditioned form (17) of (16) is not symmetric and is therefore not satisfied by the stationary point of its functional. The Untransformed Preconditioned Conjugate Gradient Method (UPCG) is explained by [4] and
The sample problem is a flat test specimen with straight edges and a center notch. The specimen is a laminate made from carbon fiber reinforced plastics (CFRP). The sketch in Fig. 2 shows the specimen region between the clamps.
26
G. Kress et al. / Composite Structures 69 (2005) 21–33
2R
2W
2L
Fig. 2. Geometry of notched test specimen: L ¼ 150 mm, w ¼ 50 mm, r ¼ 5 mm.
The specimen is loaded by moving the two clamps apart from each other. 4.1. FEM modelling Excluding non-balanced laminates renders the undamaged specimen doubly symmetric and allows to model only one quarter of the specimen without losing relevant information. It is strictly not justifiable to assume a double symmetry of the damage patterns evolving under load, as becomes obvious when studying radiographs, such as the one shown in Fig. 1, of fatigued symmetric and balanced laminates found in the literature [1]. However, as this study is not primarily concerned with damage mechanics but rather with aspects of using iterative solution techniques for improving the numerical efficiency of damage progression analyses, it is hoped that the use of the quarter model shown in Fig. 3 may be forgiven. The mesh is automatically generated with regard to two regions of the quarter model. One is the square area around the hole and the other the rectangle extending to the clamped edge. The mesh density is controlled by the three parameters specifying the number of elements Nx , Ny , and Nr along the • longitudinal direction of the outer rectangle, • transverse direction of the outer rectangle, and • radial direction of the inner square, respectively. The mesh of the inner square is designed for a smooth transition between its outer rectangular contour and the circular edge of the hole to keep distortion of the linear plane-stress finite elements minimal. These elements, or their underlying mechanical assumption of plane stress, confine the model to symmetric laminates with symmetric damage states about
Fig. 3. Mesh with Nx ¼ 25, Ny ¼ 25, and Nr ¼ 25.
the midplane. The mesh shown in Fig. 3 gives a numerical system of equations with 3962 degrees of freedom and a bandwidth of 106. Geometric boundary conditions are imposed to preserve the symmetry conditions justifying the quarter model and to simulate the loading by clamps. The symmetries require that the longitudinal and transverse displacement ux and uy must be zero along x ¼ 0 and y ¼ 0, respectively. The clamped edge of the model must remain straight and preserve its original length under load. Therefore, at x ¼ L, the displacements uy and ux are specified to be zero and unity, respectively. The actual straining of the specimen during the damage propagation simulation depends on the load factors of the various damage events. Modelling of the loading by inhomogeneous geometric boundary conditions renders the change-of-load term Dr of the gradient rf (20) zero. 4.2. Laminates This study considers the following unidirectional (UD) and multidirectional (MD) laminates: All laminates are symmetric and balanced and are made from unidirectionally reinforced prepreg made from T300 carbon fibers and an epoxy resin. Thus, the initially undamaged laminate is free from elastic couplings such as between extension and shear or extension and bending. An extension-shear coupling may be introduced by damage and is well mapped within the confinement of the quarter model which is not able to map global antisymmetric damage patterns such as shown in Fig. 1. The engineering constants Ei , Gij , and mij as well as the coefficients of thermal expansion (CTE) ai , at a fibervolume fraction of vf ¼ 60%, are listed in Table 2.
5. Numerical results Simulated damage patterns are shown in Figs. 4, 13, 17, 19, and 22. The information is color-coded. If the part of the laminate represented by a finite element has suffered matrix failure in one of its layers, the element appears green. It cannot be seen from the plots whether the laminate has suffered one or several matrix failures. Fiber failure is coded blue, and the presence of both fiber and matrix failure in one element is indicated by a red color. Elements with failure are plotted without distinct edges so that regions with contiguous damage appear without mesh lines. In black-and-white printouts, red gives a darker grey shade than green and blue appears darker than red. 1
1 For interpretation of color in Figs. 4, 13, 17, 19, and 22, the reader is referred to the web version of this article.
G. Kress et al. / Composite Structures 69 (2005) 21–33
5.1. Unzip action
27
0.08
1. the damage propagation pattern must be reproduced, 2. the failure load factors should be closely approximated. Meeting these criteria justifies addressing the relative numerical efficiency of the iterative solution methods. The sequence of first-ply-failure events shown in Fig. 4 are triggered by the load factors presented in Fig. 5. Here, the displacement fields after damage are calculated by a direct solver (DS). The curve is used as a Table 1 Laminates No.
Characterization
Stacking sequence
1 2 3 4 5
Transversely reinforced UD Longitudinally reinforced UD Symmetric cross-ply laminate Quasi-isotropic laminate Angle-ply laminate
[90] [0] [0,90]s [0,±45,90]s [±45]s
Failure Load Factors
0.07
The simplest damage pattern occurs with the UD laminate 1 specified in Table 1. Due to the normal stress concentration at the hole edge and x ¼ 0 and the low matrix strength, the first damage event is bound to be matrix failure in the element indicated in the mesh plot of Fig. 3. The propagating damage must then resemble an unzip action from the hole to the straight edge. Fig. 4 reveals that the damage modelling by reducing elastic constants in the respective principal material directions is crude because it does not well capture the realistic forming of a crack tip. However, its simplicity makes the unzip action suited for assessing the various iterative solution techniques by testing the criteria:
0.06 0.05 0.04 0.03 0.02 0.01 0.00 1
2
3
4
5
6 7
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Failure Events
Fig. 5. Failure load factors: unzip action and DS.
reference for the other curves resulting from iterative techniques for obtaining the damage-model induced changes of the displacement field. 5.1.1. Unzip action and Jacobi-inspired method JIM The sequence defined by (26) in Section 3.2 obtains the critical load factors shown in Fig. 6. The curves are obtained with numbers of iterations 1–5, and 24. The curves approximate the reference curve shown in Fig. 5 increasingly better with increasing number of iteration. The highest number of iterations requires approximately the same computing time as when using exact, or direct, solution. The damage propagation pattern labelled ’unzip action’ is always achieved, regardless of the numbers of iterations invested. A measure of how well the load factors for a sequence of failure events obtained by the iterative method fit those obtained with exact solution is provided simply by the sum of their square deviations: s¼
25 X
2
ðk kds Þ :
ð41Þ
1
The measure is plotted in Figs. 7, 9 and 16.
Property
Unit
1/23
2/13
3/12
Young’s modulus Shear stiffness Poisson’s ratio CTE
MPa MPa – K1
135000 3846 0.3 0:66
10000 5000 0.27 30:06
10000 5000 0.27 30:06
0.08 0.07
Failure Load Factors
Table 2 UD T300 properties (DORNIER SYSTEM)
0.06 0.05 0.04 0.03 0.02 0.01 0.00 1 2 3
4 5 6 7 8
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Failure Events
Fig. 4. Damage pattern unzip action after 15 failure events.
Fig. 6. Failure load factors: unzip action and JIM.
28
G. Kress et al. / Composite Structures 69 (2005) 21–33 1.E+01
1.E+02
1.E+01
1.E+00 1
2
3
4
5
6
7
8
9
10
Load Factor Deviations
Load Factor Deviation
1.E+00 1
2
3
4
5
6
7
8
9
10
1.E-01 1.E-02 1.E-03 1.E-04 1.E-05 1.E-06
1.E-01
Number of Iterations
Number of Iterations
Fig. 7. Number of JIM iterations and accuracy measure (41).
Fig. 9. Number of UPCG iterations and accuracy measure (41).
1.0
0.08
Relative Computing Times
0.9
Failure Load Factors
0.07 0.06 0.05 0.04 0.03 0.02
0.8 0.7 Exact JIM UPCG Linear (JIM) Linear (UPCG)
0.6 0.5 0.4 0.3 0.2 0.1
0.01
0.0 1
0.00 1
3
5
7
9
11
13
15
17
19
21
23
2
3
4
5
6
7
8
9
10
Number of Iterations
25
Failure Events
Fig. 10. Relative computing times for the 25/25/25 problem.
Fig. 8. Failure load factors: unzip action and UPCG.
5.1.2. Unzip action and sensitivity formula method SFM The results of the sensitivity formula method are the same as that of the first iteration step presented in Section 5.1.1. Fig. 4 obviates that the damage-induced changes of the structural behavior are much too large to be captured accurately by a linear exploitation of sensitivity. 5.1.3. Unzip action and UPCG method The Untransformed Preconditioned Conjugate Gradient (UPCG) algorithm given in Section 3.4 obtains the failure–load curves shown in Fig. 8. They represent the predicted failure loads for iteration numbers ranging from one to five. The failure loads converge rapidly against those calculated by DS, shown in Fig. 5. The curves obtained with iteration numbers higher than four cannot be distinguished from the reference curve in Fig. 8. Fig. 9 obtains an exponential decay of the deviations from the reference curve with increasing number of iterations. The figure indicates that every two iterations reduce the error by one order of magnitude. 5.1.4. Unzip action and computing times The computing times indicated in Fig. 10 refer to the mesh size 25-25-25 indicated in Fig. 4 and the unzip
action over 25 elements. They are normalized with respect to the time needed by the direct solution method. Not surprisingly, JIM requires less computational effort per iteration than UPCG. Extrapolating the curves, the computing time required by the exact solver is exceeded by UPCG after about 14 iterations and by JIM after about 24 iterations. Nevertheless, UPCG is much more efficient because of its higher convergence rate: comparing Figs. 7 and 9 with each other obtains that UPCG finds better accuracy with two iterations than JIM with ten for the here considered problem. Next, the dependence of computing times on problem size parameters is investigated. Table 3 provides data in terms number of elements Nel , number of degrees-offreedom Ndgf , half-bandwidth Nhbw , computing time with direct solver tds , and computing time with UPCG tupcg . The computation times of the various mesh sizes can be estimated by using the following model. It is based on the significant time consuming processes performed during the damage progression analysis: Tgs Tbs Tsp Tfpf
complete Gauss solution process time back-substitution processing time time required for forming inner products first-ply failure processing time
G. Kress et al. / Composite Structures 69 (2005) 21–33
29
Table 3 Problem sizes and computing times Nx
Ny
Nr
Nel
Ndgf
Nhbw
tds
tupcg
15 15 15 15 15 15 25 25 25 25 25 25 50 50 50 50 50 50
15 15 25 25 50 50 15 15 25 25 50 50 15 15 25 25 50 50
25 50 25 50 25 50 25 50 25 50 25 50 25 50 25 50 25 50
975 1725 1625 2875 3250 5750 1125 1875 1875 3125 3750 6250 1500 2250 2500 3750 5000 7500
2092 3642 3432 5982 6782 11832 2412 3962 3952 6502 7802 12852 3212 4762 5252 7802 10352 15402
66 66 106 106 206 206 66 66 106 106 206 206 66 66 106 106 206 206
611 1031 1510 2570 9170 16500 671 1030 1690 2771 10300 17640 831 1179 2150 3224 13200 20500
484 708 896 1414 2760 4743 532 730 1000 1550 3130 5082 642 848 1350 1795 4040 5970
The number of floating-point operations of the Gauss triangulation depends linearly on the number of degreesof-freedom Ndgf and quadratically on the half-bandwidth Nhbw : 2 þ 32Nhbw 1 : ð42Þ tgs ¼ aNdgf 2Nhbw
tds ¼ Nde ðtgs þ tfpf Þ:
The number of floating-point operations for the backsubstitution process depends linearly on the number Ndgf and linearly on the half-bandwidth Nhbw :
The model parameter values achieving the fit illustrated in Fig. 11 are listed in Table 4. The fitting of the factors reveals that the scalar products performed by the UPCG do not significantly contribute to the computation times. Fig. 11 shows the measured times obtained with both solution methods. The data points are represented with the open symbols. It also shows the predicted times for both methods and these points are indicated with the solid symbols.
ð43Þ
The number of floating-point operations for calculating the scalar products required in the UPCG algorithm depends linearly on the number of degrees-of-freedom Ndgf and the number of iterations invested for the convergence of the displacement fields: tsp ¼ cNdgf Niter :
tupcg ¼ tgs þ ðNde 1Þðtbs þ tsp þ tfpf Þ þ tfpf :
ð44Þ
The number of operations involved in the first-ply-failure analysis to identify the next damage event and the associated critical load depends linear on the number of finite elements Nel : tfpf ¼ dNel :
When the here proposed method is used where the displacement fields are approximated by UPCG, the computation time is estimated at
ð45Þ
The computation time for post processing, i.e. calculation of stress fields from the displacement fields, scales also with the number of finite elements and may therefore be considered a part of tfpf . The damage progression simulation is an iteration over the specified number Nde of damage events. Each step consists of two stages: solving the system of equations and identifying a new damage event. The total computation time required when the displacements fields are calculated by the direct solution method is now estimated at
ð47Þ
2.50E+04
Direct Solver Computing Times
tbs ¼ bNdgf Nhbw :
ð46Þ
2.00E+04
DS Time Data DS Time Model UPCG Time Data UPCG Time Model
1.50E+04
1.00E+04
5.00E+03
0.00E+00 1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16 17 18 19 20
Ordered Sample Problems
Fig. 11. Model fit to computing times.
Table 4 Fitting parameters for fit Fig. 11 a
b 6
0.58 · 10
4
0.49 · 10
c
d
0.0
0.9 · 102
30
G. Kress et al. / Composite Structures 69 (2005) 21–33
For presentation in Fig. 11, the data of Table 3 have been rearranged in ascending order of the times tds . The fitted models are used for an assessment of the model-size dependent computation time savings which are possible by using the UPCG approach (see Fig. 12). The time savings are given in terms of the ratio between tds and tupcg . Let us use a problem-size parameter N which equals the number of unknowns Ndgf ¼ N . Having the currently considered model with linear plane elements in mind, the number of elements may roughly be set equal to the number of nodal points and, since these each carry two degrees-of-freedom, we have Nel ¼ 0:5N . The half-bandwidth is assumed to grow by the square pffiffiffiffi root of the problem-size parameter, or Nhbw ¼ N , which is in approximate agreement with the data in Table 3. The number of damage events necessary to give complete information, on the load carrying capability or the damage tolerance of a structure, is somewhat arbitrarily assumed to grow also by pthe ffiffiffiffi square root of the problem size parameter, Nde ¼ N . Under these assumptions the models obtain that pffiffiffiffi tupcg 1 ð2a þ bÞN þ ð1:5a þ 0:5dÞ N a pffiffiffiffi : ð48Þ ¼ pffiffiffiffi tds N 2aN þ 1:5a N þ 0:5d a The computing times ratio is plotted in Fig. 12 over a range 103 6 N 6 107 . The computing times using UPCG may be larger than those using DS for rather small problem sizes. At very large problem sizes N the here proposed damage progression analysis with p UPCG gains efficiency over ffiffiffiffi that using DS bypaffiffiffiffifactor N . This remains true when, instead of Nde ¼ N , Nde ¼ N is assumed.
event is fiber failure found in the element on the corner at the edge of the hole. The damage pattern shown in Fig. 13 is obtained with the direct as well as with the preconditioned conjugate gradient solution methods. Instead of using a fixed number of iterations, the UPCG algorithm is terminated as soon as qnew < 2 q0 applies. This has the advantage that only as many iterations are performed as needed for a necessary accuracy. The number of iterations for each displacement-field update tends to increase with increasing accuracy. The five curves in Fig. 14 show, in ascending order, the number of iterations needed to achieve the terminating -values ¼ 101 , ¼ 102 , ¼ 103 , ¼ 104 , and ¼ 105 . Fig. 14 also reveals that numbers of iterations increase with increasing number of damage events which may be due to decreasing quality of the preconditioner which remains, symbolically, the inverse of the stiffness matrix related to the undamaged structure. The damage pattern is already correctly obtained with an -value as low as ¼ 101 . Fig. 15 shows the load factors after DS and UPCG for an -value range 101 P P 105 . In contrast to unzip-action, the load curve does not monotonically decrease. UPCG approximates the loads quite well and the mean square errors shown in Fig. 16 indicate that -values smaller than about ¼ 103 does not further improve them.
5.2. Complex laminate damage patterns 5.2.1. Unidirectional laminate [0] Early damage pattern development in this laminate is dominated by matrix shear failure roughly along the tangents at the hole. This is due to a combination of high shearing stresses along the tangents and low matrix shear strength of the material. Only the last damage
Fig. 13. Damage in [0]s -laminate after 25 failures.
25
1.0 1.0E+03
1.0E+04
1.0E+05
1.0E+06
0.1
1.0E+07
PCG Iterations per Update
PCG over DS solution times
10.0
20
15
10
5
0
1
0.0
Number of Degrees of Freedom
Fig. 12. Predicted computing times ratio tupcg =tds .
3
5
7
9
11
13
15
17
19
21
23
Failure Events
Fig. 14. Numbers of iterations to achieve residual norms.
25
G. Kress et al. / Composite Structures 69 (2005) 21–33 0.10
Failure Load Factors
0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 1
3
5
7
9
11
13
15
17
19
21
23
25
Failure Events
Fig. 15. UPCG load factors during [0] damage simulation.
1.E+00 -1
-2
-3
-4
-5
Load Factor Deviations
1.E-01 1.E-02 1.E-03 1.E-04 1.E-05 1.E-06 1.E-07
31
The two damage patterns are similar but not identical. Fig. 18 reveals that both solution methods calculate almost identical values for the first ten critical loads but afterwards they find different values. The damage development seems to diverge at this point. It is believed that the damage development is not robust in the sense that, when several critically stressed layers exist, small differences in the previously calculated stiffness matrices may play a deciding role in which of the critically stressed elements will be recognized to fail first. The outcome of such decision will naturally alter the subsequent damage development. This has nothing to do with the accuracy of the iterative solution method but rather with other details of the code here used for the sample calculations. 5.2.3. Quasi-isotropic laminate [0, 45,90]s The striking difference of the damage pattern evolved in Fig. 19 to the previous ones is that the region with matrix failure is more confined towards the corner formed by the hole and the symmetry line. This is due to the arresting effect of the fiber reinforcement which is effective in all directions. After much matrix damage is accumulated, fiber fracture occurs and threatens the integrity of the laminate similar to the Unzip Action. The
1.E-08
Termination criterion value epsilon
0.10 0.09
5.2.2. Cross-ply laminate [0,90]s Multi-layered laminates can suffer several damage events within one finite element. The number of damage events leading to the damage pattern of the cross-ply laminate shown in Fig. 17 is one hundred. Similarly to the longitudinally reinforced unidirectional laminate, early damage is dominated by matrix failure. The failures in the 0-layer are caused by shear und those in the 90-layer are caused by normal stresses transverse to the fibers. At the end, fiber fracture appears.
Fig. 17. Damage in [0,90]s -laminate after 100 failures.
Failure Load Factors
Fig. 16. UPCG failure load deviations (41) versus -values.
0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 1
3
5
7
9
11
13
15
17
19
21
23
25
Failure Events
Fig. 18. DS and UPCG load factors during [0,90]s damage simulations.
Fig. 19. Damage in [0,±45,90]s -laminate after 100 failures.
32
G. Kress et al. / Composite Structures 69 (2005) 21–33
load curves plotted in Fig. 20 reveal that the first 24 damage events are found identical by both solution methods. Again, UPCG converges quickly and a lowprecision termination with ¼ 101 already obtains accurate results, and the divergence between the damage patterns obtained by both methods is suspected to be due to some numerical deviation of the stiffness matrix update rather than lack of accuracy of UPCG. The quick convergence, even at weak accuracy demands on UPCG termination, is due to the fact that displacement field perturbations caused by matrix failure in layers of the quasi-isotropic laminate must be small because the fiber reinforcement restrains equally ‘‘in all’’ directions (hence quasi-isotropic). Therefore, it is interesting to consider the question of whether the displacements changes encountered here are small enough to fall within the application regime of the sensitivity-formula-inspired method (SF) derived in Section 3.3. It has been found that the critical loads evaluated by SF cannot be distinguished by the eye from the curves obtained by UPCG plotted in Fig. 20. Fig. 21 illustrates the computing times when using direct, the UPCG with only one iteration, and the SF solution techniques, on two different model sizes. The smaller
problem and the larger problems have the mesh sizes Nx ; Ny ; Nr ¼ 25; 25; 25 and Nx ; Ny ; Nr ¼ 50; 50; 50, respectively, and the corresponding numerical sizes can be looked up in Table 3. The possible time savings of SF with respect to UPCG restricted to only one iteration appear small which may be due to the similar preconditioning techniques. The UPCG appears to give slightly better accuracy even after the first iteration so that the case for SF is relatively weak. However, both methods save much time if compared with the direct solution and Fig. 21 illustrates how quickly the effect increases with increasing problem size. 5.2.4. Laminate [ 45]s The ½ 45 s -fiber reinforcement allows the matrix damage to spread as shown in Fig. 22. The damage zone extends from the notch to the free edge, reminding of the distinct ±45 cracks of the damage pattern shown in Fig. 1. It must be noted that the quarter model with its symmetry assumption prevents the anti-symmetric shape appearing in the figure. The two solution techniques produce similar damage patterns although the load plots in Fig. 23 reveal that the damage paths of the two solution methods separate quite early.
0.12
Failure Load Factors
0.10 0.08
0.06 0.04 0.02 0.00 1
3
5
7
9
11
13
15
17
19
21
23
25
Failure Events
Fig. 22. Damage in [±45]s -laminate after 200 failures.
Fig. 20. Load factors of [0,±45,90]s -laminate.
50000
0.16
45000 0.14
0.12
35000 30000
Exact Solution Conjugate Gradient Method "Sensitivity Formula"
25000 20000 15000 10000
Failure Load
Computing Time
40000
0.10
0.08 0.06 0.04
5000 0.02
0 25
50
Methods at Mesh Sizes
0.00 1
2
3
4
5
6
7
8
9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Number of Damage Event
Fig. 21. Computing times with DS, UPCG, and SF at two problem sizes.
Fig. 23. Load factors of [±45]s -laminate.
G. Kress et al. / Composite Structures 69 (2005) 21–33
6. Conclusion Approximative and iterative solution methods have been evaluated in context with damage progression simulation. The efficiency of the approximate solution techniques is enhanced by rearranging the system of equations so that it represents damage stages of the structure and uses as primary unknowns the damageinduced changes of the displacement field of the virgin structure. The solution methods have been assessed by performing damage progression simulations. The damage modelling, as such rather crude, suffices to serve the purpose. Several unidirectional as well as multidirectional laminates have been considered. A one-step approximation method inspired by the sensitivity formula for finite element models, often used in the context of gradient-based design optimization methods, gives good results when simulating the effects of matrix cracking in multidirectional laminates where the fiber arrangements prevent damage-induced stiffness variations to become large. The considered iterative solution methods, Untransformed Preconditioned Conjugate Gradient Method (UPCG) and a method inspired by the Jacobi Method (JIM), are applicable to all types of laminate. They can reproduce any damage pattern resulting from using exact solution of the state equations. Not surprisingly, UPCG converges much quicker than JIM. Thanks to the preconditioning and a good starting vector, very good accuracy is obtained with a very small number of iterations if compared to the number of unknowns. The computing times, resulting from using a frontal solver and from using UPCG with five iterations, have been modelled and the free parameters adjusted to fit computing times obtained on a laptop computer. The impact of solution technique on the numerical effort of
33
damage progression simulation has been assessed. Under some assumptions on the growth characteristics of typical problems, it has been found that UPCG reduces thepoverall numerical-effort growth scale by a ffiffiffiffi factor of N if compared to using direct solution techniques. The most striking result is that UPCG, embedded in damage progression simulation and enhanced by the here suggested precondition technique, performs significantly better than direct solution methods although the encountered matrices have a size and a banded structure where direct solution methods could otherwise not be challenged by any iterative method.
Acknowledgements Efficient solution techniques in the presence of small structural stiffness matrix variations have suggested themselves in the course of our research on design optimization with Evolutionary Algorithms, so that grant no. 2100-066879.01/1 of the Swiss National Science Foundation is again gratefully acknowledged.
References [1] Kress GR, Stinchcomb WW 1985. Fatigue response of notched graphite/epoxy laminates. In: Vinson JR, Taya M, editors. Recent Advances in Composites in the United States and Japan, ASTM STP 864, American Society for Testing and Materials, Philadelphia. p. 173–96. [2] Hashin Z. Failure criteria for unidirectional fiber composites. J Compos Mater 1980;7:329–34. [3] Kress G. Free-edge influence on CFRP-laminate strength. Int J Damage Mech 1994;3:192–211. [4] Shewchuk JR. An introduction to the gradient method without the agonizing pain. Pittsburgh PA, 15213: School of Computer Science, Carnegie Mellon University; 1994.
Iterative solution methods for damage progression analysis G. Kress *, M. Siau, P. Ermanni Structure Technologies, ETH Z€urich, LEO C2, Leonhardstr. 27, CH-8092 Z€urich, Switzerland Available online 2 June 2004
Abstract Damage progression is a sequence of damage events. A simulation of it solves a sequence of systems of equations which differ from each other only by the effects of the modelling of the incrementally increasing damage. An initial matrix triangulation, to solve the equations of the virgin structure, provides a preconditioner that can be used for solving the sequence of system equations by approximate or iterative solution techniques. In addition, obtained solutions serve as starting vectors for the next solutions. The numerical efficiency potential of this concept is evaluated by considering a method inspired by the sensitivity formula for finite element models, an iterative method inspired by Jacobi’s method, and the untransformed preconditioned conjugate gradient method. The ability to reproduce the damage patterns obtained by using exact solution is checked and the efficiency improvements are assessed. 2004 Elsevier Ltd. All rights reserved. Keywords: FEM; Numerical methods; Damage simulation; Composites
1. Introduction Typically, iterative methods for solving linear systems of equations can outperform direct solution methods when the number of unknowns is very high and the coefficient matrix is sparsely populated. Iterative solution methods may be competitive even for smaller problems with a banded coefficient matrix when the sought solution is closely guessed by a starting vector and the coefficient matrix is well conditioned. These circumstances can be provided at reasonable cost when variations of the same problem must be evaluated sequentially. Examples include the sequential refinement of design solutions towards some specified optimum or the progression of damage in composite materials in small steps. 1.1. Damage in composite materials Fig. 1 shows a state of damage in a notched test specimen made from a laminate with the stacking sequence [45,90,)45,0]s made from a T300/5208 graphite/ epoxy prepreg tape. The specimen suffered 104 cycles at a maximum stress of 85% of their mean quasi-static tensile strength. The radiograph shows that matrix *
Corresponding author. Fax: +41-1-633-1125. E-mail address: [email protected] (G. Kress).
0263-8223/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2004.04.043
cracking in all layers as well as delaminations exist. The delaminations appear in the radiograph in darker shades than the undamaged parts of the laminate. Characteristic features are the matrix cracks in the longitudinal direction tangent to the hole in the 0-degree plies, the matrix cracks emanating from the hole and extending almost to the straight edges of the specimen, and the delaminations seeming to be companions of the tangent matrix cracks. The fact that the diagonal matrix cracks extend much farther in the 45- than in the )45-direction has to do with the stacking sequence placing the 45degree layer at the surface. Please note that this is an effect that could only be captured by a fully threedimensional model, which also holds true for the delaminations. Obviously, it is not an easy task to simulate the damage development under fatigue loading. For ease, a very simple damage model is considered sufficient for the efficiency study although more sophisticated models exist. 1.2. Overview of the present study The present study chooses composite material damage progression analysis just as one of the various possible backdrops to verify approaches for taking advantage of iterative solution techniques for simulations involving a sequence of structural models which
22
G. Kress et al. / Composite Structures 69 (2005) 21–33
accuracy, and computing time. With the other laminates more complex damage patterns are simulated in order to compare the results of the direct and iterative solution methods. The conclusions of this study are found in Section 6.
2. Progressive damage analysis 2.1. Mechanical and residual stresses
Fig. 1. Damage in notched [45,90,)45,0]s -specimen [1].
differ from each other by small variations. Section 2 presents a very simply damage analysis approach conceptually including the effects of residual stresses due to temperature or moisture loads. Failure criteria are taken from [2] but the effect of shear stress on tensile fiber failure is omitted to achieve damage patterns resembling more closely those due to fatigue. The damage is modelled by drastically reducing stiffness parameters in directions perpendicular to the found failure lines. Simple as it is, the modelling has led to fair predictions of quasistatic strength of specimens with severe edge effects [3]. The alteration of the numerical system, initially representing the virgin structure, by damage modelling is considered in Section 3. From this, convenient formulations of the change-of-displacement field are presented and preconditioning with hindsight to iterative solution techniques is explained. Specifically, an iteration inspired by the Jacobi Method and a suited version of the Conjugate Gradient Method are considered. Attention is also given to a method derived from the well-known sensitivity formula for gradient calculation of the total potential energy. It could be used to hurry most quickly through the steps of an iterative damage progression. The sample problem structure is a notched test specimen and Section 4 explains the FEM modelling. It maps only a plane-stress state and assumes symmetric damage states which is clearly controversial as Fig. 1 testifies. The laminate stacking sequences and material properties used in the simulations are also given in this section. Numerical results are presented in Section 5. The simplest damage propagation pattern, here dubbed Unzip Action, is reproduced by all considered solution methods. It is used to investigate aspects of convergence,
The unidirectionally reinforced composite is modelled as a transversely isotropic homogeneous material with different properties in the fiber direction and those transverse to it. When temperature changes, the free thermal expansion of an individual layer is constrained by the other layers of the multidirectional laminate. During the autoclave process, the temperature of the laminate increases but the viscosity of the matrix material prevents temperature stresses. When the matrix material cures, internal residual stresses rres are initially negligible but increase with falling temperature. When the difference DT to the stress-free reference temperature becomes very large, laminates made from composites with brittle matrices may experience matrix cracking even before any external load is applied. In order to find the mechanical load causing damage, the designer must therefore take into account the presence of residual stresses before application of any boundary conditions. We assume a situation where the operating temperature and the mechanical load case are specified and the designer must answer the question at what multiple of the mechanical loads damage or failure of the part will occur. This leads to the partition of the stresses at any point in the laminated structure r ¼ krmec þ rres ;
ð1Þ
where the load factor k is obtained from substituting (1) into suitable failure criteria so that f ðrÞ ¼ 1:
ð2Þ
The theory of laminated plates describes the relationship between combined mechanical and thermal loading and the plate deformations with
N M
A ¼ B
B C
0 j
by using so-called fictitious line loads T N N N ¼ þ : M M MT
ð3Þ
ð4Þ
The stresses rres follow from a FEM analysis where the temperature- and moisture caused line loads NT are calculated on element level and assembled to the righthand-side vector. The essential boundary conditions are
G. Kress et al. / Composite Structures 69 (2005) 21–33
defined to give a statically determinate support to allow for the free deformation of the structure. The resulting deformation along the boundary Cu with fixed displacements must be taken into account when implementing the homogeneous geometric boundary conditions. They must not change the stress state of the free structure caused by the temperature- and moisture loads. Then, the stresses rmec follow from a FEM analysis where the geometric boundary conditions and only the mechanical external forces are applied. 2.2. Failure criteria As a matrix crack and a fiber rupture will differently change the load-carrying behavior of the laminate, a suitable failure criterion must be able to predict the failure mode. It should also be able to take into account the effect of combined stressing on the failure probability. The present study uses the four criteria by [2] for fiber and matrix failure under tensile and compressive loading under plane-stress loading. The strength parameters are obtained from single-stress tests: tensile strength in fiber direction compressive strength in fiber direction tensile strength in transverse direction compressive strength in transverse direction shear strength longitudinal shear strength transverse
XT XC YT YC SL ST
r1 XT
2
þ
s21 SL
ð5Þ
Fiber failure under compressive fiber stress (r1 < 0): r1 þ XC ¼ 0: Matrix failure under tensile stress (r2 > 0): 2 2 r2 s12 þ ¼ 1: YT SL Matrix failure under compressive stress (r2 < 0): # 2 " 2 2 r2 YC r2 s12 þ ¼ 1: 1 2ST 2ST YC SL
Critical factors on the mechanical loading are obtained from combining (1) and (2) with a suitable criterion. Quadratic failure criteria then lead to quadratic equations for the critical load k2 Q þ kL þ P ¼ 1:
ð9Þ
A value of P greater than 1 indicates that residual stresses be high enough to cause failure. Then, the damage resulting from this failure must be modelled but the factor on the mechanical load case is zero. Else, k is calculated by k¼
L þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi L2 þ 4Qð1 P Þ 2Q
ð10Þ
or, if the criterion is actually linear, Q ¼ 0, by k¼
1P : L
ð11Þ
Failure is recognized element-wise. The intralaminar damage modes of fiber fracture and matrix failure are crudely modelled by reducing the respective Young’s moduli in the respective layers to a very small fraction of their original values, responding to the lost capability of transferring load across the fractured surfaces:
2 ¼ 1:
2.3. Prediction of critical load factors
2.4. Damage modelling
Fiber failure under tensile fiber stress (r1 > 0):
23
ð6Þ
G12 ¼ G12 ;
m12 ¼ m21 ¼ 0;
c: fiber:
G12 ¼ G12 ;
m12 ¼ m21 ¼ 0;
t: matrix: E2 ¼ E1 ; c: matrix:
G12 ¼ G12 ; G12 ¼ G12 ;
m12 ¼ m21 ¼ 0; m12 ¼ m21 ¼ 0:
t: fiber:
E1 ¼ E1 ;
ð12Þ 2.5. Termination criteria
ð7Þ
ð8Þ
The simulations presented in Section 5 are based on one rather arbitrary alteration of Hashin’s criteria: the shear-stress term is removed from the criterion for fiber fracture under tensile stress (5). This leads to damage patterns more similar to those appearing under cyclic loading, see Fig. 1, where extensive damage of the matrix occurs before fibers fail.
The damage progression analysis is a repeatedly applied first-ply-failure analysis. After a failure event is detected, the structural model is updated to account for the damage and, due to the altered stress state, a new failure event will be detected at a different location. The iterations stop upon fulfillment of one of the following criteria: • the number of damage events reaches a preset value, • the work of the external forces exceeds a preset value. The latter criterion can also be used to recognize total failure where the structure ceases to be one part.
24
G. Kress et al. / Composite Structures 69 (2005) 21–33 1
3. Damage states and displacement fields
ðI þ Kð0Þ DKÞD~u ¼ D~u þ LD~u;
Each damage event reduces the apparent material stiffness of the region surrounding it. Considering a progressive failure analysis, where the behavior of the structural model is simulated by a finite-element model,
where the operation LD~u first assembles the vector DKD~u by matrix multiplications on element level and then applies the back-substitution process on it.
K~u ¼ r;
3.1. Functional and gradient
ð13Þ
the region surrounding a new damage event is presented by one finite element. Its stiffness matrix will change, representing the loss of material stiffness, and possibly internal forces will be set free in the presence of residual stresses. This will also change the displacements of the structure and with it the local material deformations so that each damage influences the spatial probability distribution of subsequent failure events. All of these changes are captured in the following equation: ðKð0Þ þ DKÞð~ uð0Þ þ D~ uÞ ¼ ðrð0Þ þ DrÞ:
ð14Þ
ð18Þ
Some iterative methods obtain the displacement solution of (16) by minimizing the quadratic functional f ðD~uÞ i h f ¼ D~uT 12ðKð0Þ þ DKÞD~u ðDr DK~uð0Þ Þ : ð19Þ This is the same as seeking the zero-point of its gradient rf ¼ Kð0Þ D~u þ DKð~uð0Þ þ D~uÞ Dr:
ð20Þ
3.2. Jacobi method inspired iterative solution method
ð0Þ
refers the respective entities to the The upper index undamaged state of the structure. Removing the system equations referring to the initial, or undamaged, state of the simulated structure
The Jacobi Method transforms a given system of Eq. (13), by partitioning the original coefficient matrix K into the diagonal matrix D and the remainder E,
Kð0Þ ~ uð0Þ ¼ rð0Þ
K~u ¼ ~r;
ð15Þ
~u ¼ D E~u þ D1~r;
reduces (14) to the form ðKð0Þ þ DKÞD~ u ¼ Dr DK~ uð0Þ :
ð16Þ
To evaluate (16), let it be that the system (15) has already been solved to obtain the displacement field ~uð0Þ . This implies that the system matrix Kð0Þ can be inverted because the appropriate essential boundary conditions are implemented. While (14) obtains the total displacements corresponding to a given instant of the damage history, (16) obtains the changes from the initial displacements found by (15). The changes D~ u due to the consecutive failure events incrementally modify the initial solution ~ uð0Þ obtained by (15). If the differences between two consecutive displacement fields are small, the numerical effort of iterative solution techniques may become less than that of direct ones because the previous displacements provide a good starting vector for the new ones. Furthermore, the triangulized form of Kð0Þ obtained from the direct solution of (15) can be used as a preconditioner for improving the convergence of iterative methods for approximating the solution of (16). Applying Kð0Þ to (16) shows1 that only the damaged elements cause deviation Kð0Þ DK from the ideal unity matrix: 1
1
ðI þ Kð0Þ DKÞD~ u ¼ Kð0Þ ðDr DK~ uð0Þ Þ:
D~u ¼ E~u þ ~r; 1
ð21Þ
so that an iteration rule results: ~uiþ1 ¼ B~ui þ ~z:
ð22Þ
The iterations converge to the exact solution because the convergence radius of B is smaller than unity or qðBÞ ¼ max jkj j < 1 due to the fact that the values of the entries of the main diagonal D of K are larger than those of the off-diagonal entries. The operation B~uj calculates the entries of the new vector ~ujþ1 sequentially, starting from its first entry ð~uðjþ1Þ Þ1 . The convergence is enhanced by replacing the entries of the previous vector ~uj with the already calculated entries of ~ujþ1 . The improving feature distinguishes the Gauss-Seidel Method from the Jacobi Method. Re-ordering the terms in (17) provides a form similar to that in the last line of (21) 1
1
D~u ¼ Kð0Þ DKD~u þ Kð0Þ ðDr DK~uð0Þ Þ; 1
ð23Þ
1
where Kð0Þ DK and Kð0Þ ðDr DK~uð0Þ Þ have similar meanings as B and z, respectively. The equation can be solved in iterations analogous to (22): 1
D~uiþ1 ¼ Kð0Þ Dr Lð~uð0Þ þ D~ui Þ
ð24Þ
or, if D~u0 ¼ 0,
ð17Þ 1
The deviation, however, is 1 generally not symmetric even though both matrices, Kð0Þ and DK, are symmetric. However, the left-hand side of (17) can also be evaluated as
~u ¼ B~u þ ~z
D~uiþ1 ¼ Kð0Þ Dr þ
i X
i
ð1Þ Li ~u0 :
ð25Þ
1
The sum of the displacements u~ð0Þ of the undamaged structure and the changes D~uðpÞ due to previous damage
G. Kress et al. / Composite Structures 69 (2005) 21–33
events provides a better starting vector. Then, the iterations yield 1
D~uiþ1 ¼ Kð0Þ Dr þ
i X i ð1Þ Li ~ uðpÞ : u0 Li D~
ð26Þ
1
The series converges linearly if it holds for any vector v that i
jL vj < jvj:
ð27Þ
The method presented here is not identical with the Jacobi Method: the preconditioner is the inverse of the complete initial coefficient matrix, and it is applied to the varied stiffness matrices of the damage progression simulation sequence. Therefore, it is here referred to as Jacobi-inspired method (JIM). The iterations are continued with i ¼ i þ 1 until a maximum number is reached or the termination criterion indicates that the chosen accuracy measure is fulfilled. 3.3. Method derived from sensitivity formula
1
solves the problem of using the preconditioner Kð0Þ to improve convergence to the solution of (16) and (17). The gradient rf0 at the starting point of the search is calculated directly from (20). The initial search direction follows from the operation 1
s0 ¼ Kð0Þ rf0
ð30Þ
performing back substitution with the existing triangulized stiffness matrix of the undamaged structure. The factor qnew is the inner product of gradient and search direction qnew ¼ rf0T s0
ð31Þ
and it is needed for both step-length and conjugate search direction calculations. The initial value is stored q0 ¼ qnew
ð32Þ 2
for use in the termination criterion qnew < q0 stopping the following iterative part of the algorithm. A vector q follows from applying the stiffness matrix K to the current search direction: qi ¼ Ksi :
If the damage modelling leads only to small stiffness changes, the changes of displacements are also small compared with the displacements of the undamaged structure. This may be the case for structures made from multidirectional laminates consisting of fiber-reinforced materials whose limit strain in the transverse direction is much smaller than in the fiber direction and when only matrix failure occurs. Then, the nonlinear term, requiring iterations, may be dropped so that the linear form results: 1
D~u ¼ Kð0Þ ðDr DK~ uð0Þ Þ:
ð28Þ
This can also be derived from the formula used to evaluate the sensitivity of the displacement field on some structural parameters x: 1 o~u or oK ð0Þ ~ ¼ Kð0Þ ; ð29Þ u oxi oxi oxi where the xi are structural stiffness parameters. From the reference point, i.e. the displacement field of the undamaged structure, a step along the direction of the gradient must be taken the length of which is given by the amount of stiffness change of the element where damage is modelled. This leads again to (28). Obviously, the method (SFM) gives the same result as the first iteration step of the method presented in Section 3.2.
25
ð33Þ
The step length q ai ¼ Tnew s i qi
ð34Þ
leads from the current reference point D~ui1 to the new point D~ui minimizing f (19) along the search direction si : D~ui ¼ D~uði1Þ þ ai si :
ð35Þ
The gradient is updated by: rfi ¼ rfi1 þ ai qi :
ð36Þ
A vector pi is formed by applying back substitution on the gradient rfi pi ¼ Kð0Þ1 rf ðiÞ :
ð37Þ
The factor q is updated by qold ¼ qnew ;
qnew ¼ rfiT pi :
The factor bðiÞ , q bðiÞ ¼ new qold
ð38Þ
ð39Þ
renders consecutive search directions si ¼ pi þ bsi1 ;
ð40Þ
K-conjugate.
4. Sample problem 3.4. Conjugate gradient method The preconditioned form (17) of (16) is not symmetric and is therefore not satisfied by the stationary point of its functional. The Untransformed Preconditioned Conjugate Gradient Method (UPCG) is explained by [4] and
The sample problem is a flat test specimen with straight edges and a center notch. The specimen is a laminate made from carbon fiber reinforced plastics (CFRP). The sketch in Fig. 2 shows the specimen region between the clamps.
26
G. Kress et al. / Composite Structures 69 (2005) 21–33
2R
2W
2L
Fig. 2. Geometry of notched test specimen: L ¼ 150 mm, w ¼ 50 mm, r ¼ 5 mm.
The specimen is loaded by moving the two clamps apart from each other. 4.1. FEM modelling Excluding non-balanced laminates renders the undamaged specimen doubly symmetric and allows to model only one quarter of the specimen without losing relevant information. It is strictly not justifiable to assume a double symmetry of the damage patterns evolving under load, as becomes obvious when studying radiographs, such as the one shown in Fig. 1, of fatigued symmetric and balanced laminates found in the literature [1]. However, as this study is not primarily concerned with damage mechanics but rather with aspects of using iterative solution techniques for improving the numerical efficiency of damage progression analyses, it is hoped that the use of the quarter model shown in Fig. 3 may be forgiven. The mesh is automatically generated with regard to two regions of the quarter model. One is the square area around the hole and the other the rectangle extending to the clamped edge. The mesh density is controlled by the three parameters specifying the number of elements Nx , Ny , and Nr along the • longitudinal direction of the outer rectangle, • transverse direction of the outer rectangle, and • radial direction of the inner square, respectively. The mesh of the inner square is designed for a smooth transition between its outer rectangular contour and the circular edge of the hole to keep distortion of the linear plane-stress finite elements minimal. These elements, or their underlying mechanical assumption of plane stress, confine the model to symmetric laminates with symmetric damage states about
Fig. 3. Mesh with Nx ¼ 25, Ny ¼ 25, and Nr ¼ 25.
the midplane. The mesh shown in Fig. 3 gives a numerical system of equations with 3962 degrees of freedom and a bandwidth of 106. Geometric boundary conditions are imposed to preserve the symmetry conditions justifying the quarter model and to simulate the loading by clamps. The symmetries require that the longitudinal and transverse displacement ux and uy must be zero along x ¼ 0 and y ¼ 0, respectively. The clamped edge of the model must remain straight and preserve its original length under load. Therefore, at x ¼ L, the displacements uy and ux are specified to be zero and unity, respectively. The actual straining of the specimen during the damage propagation simulation depends on the load factors of the various damage events. Modelling of the loading by inhomogeneous geometric boundary conditions renders the change-of-load term Dr of the gradient rf (20) zero. 4.2. Laminates This study considers the following unidirectional (UD) and multidirectional (MD) laminates: All laminates are symmetric and balanced and are made from unidirectionally reinforced prepreg made from T300 carbon fibers and an epoxy resin. Thus, the initially undamaged laminate is free from elastic couplings such as between extension and shear or extension and bending. An extension-shear coupling may be introduced by damage and is well mapped within the confinement of the quarter model which is not able to map global antisymmetric damage patterns such as shown in Fig. 1. The engineering constants Ei , Gij , and mij as well as the coefficients of thermal expansion (CTE) ai , at a fibervolume fraction of vf ¼ 60%, are listed in Table 2.
5. Numerical results Simulated damage patterns are shown in Figs. 4, 13, 17, 19, and 22. The information is color-coded. If the part of the laminate represented by a finite element has suffered matrix failure in one of its layers, the element appears green. It cannot be seen from the plots whether the laminate has suffered one or several matrix failures. Fiber failure is coded blue, and the presence of both fiber and matrix failure in one element is indicated by a red color. Elements with failure are plotted without distinct edges so that regions with contiguous damage appear without mesh lines. In black-and-white printouts, red gives a darker grey shade than green and blue appears darker than red. 1
1 For interpretation of color in Figs. 4, 13, 17, 19, and 22, the reader is referred to the web version of this article.
G. Kress et al. / Composite Structures 69 (2005) 21–33
5.1. Unzip action
27
0.08
1. the damage propagation pattern must be reproduced, 2. the failure load factors should be closely approximated. Meeting these criteria justifies addressing the relative numerical efficiency of the iterative solution methods. The sequence of first-ply-failure events shown in Fig. 4 are triggered by the load factors presented in Fig. 5. Here, the displacement fields after damage are calculated by a direct solver (DS). The curve is used as a Table 1 Laminates No.
Characterization
Stacking sequence
1 2 3 4 5
Transversely reinforced UD Longitudinally reinforced UD Symmetric cross-ply laminate Quasi-isotropic laminate Angle-ply laminate
[90] [0] [0,90]s [0,±45,90]s [±45]s
Failure Load Factors
0.07
The simplest damage pattern occurs with the UD laminate 1 specified in Table 1. Due to the normal stress concentration at the hole edge and x ¼ 0 and the low matrix strength, the first damage event is bound to be matrix failure in the element indicated in the mesh plot of Fig. 3. The propagating damage must then resemble an unzip action from the hole to the straight edge. Fig. 4 reveals that the damage modelling by reducing elastic constants in the respective principal material directions is crude because it does not well capture the realistic forming of a crack tip. However, its simplicity makes the unzip action suited for assessing the various iterative solution techniques by testing the criteria:
0.06 0.05 0.04 0.03 0.02 0.01 0.00 1
2
3
4
5
6 7
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Failure Events
Fig. 5. Failure load factors: unzip action and DS.
reference for the other curves resulting from iterative techniques for obtaining the damage-model induced changes of the displacement field. 5.1.1. Unzip action and Jacobi-inspired method JIM The sequence defined by (26) in Section 3.2 obtains the critical load factors shown in Fig. 6. The curves are obtained with numbers of iterations 1–5, and 24. The curves approximate the reference curve shown in Fig. 5 increasingly better with increasing number of iteration. The highest number of iterations requires approximately the same computing time as when using exact, or direct, solution. The damage propagation pattern labelled ’unzip action’ is always achieved, regardless of the numbers of iterations invested. A measure of how well the load factors for a sequence of failure events obtained by the iterative method fit those obtained with exact solution is provided simply by the sum of their square deviations: s¼
25 X
2
ðk kds Þ :
ð41Þ
1
The measure is plotted in Figs. 7, 9 and 16.
Property
Unit
1/23
2/13
3/12
Young’s modulus Shear stiffness Poisson’s ratio CTE
MPa MPa – K1
135000 3846 0.3 0:66
10000 5000 0.27 30:06
10000 5000 0.27 30:06
0.08 0.07
Failure Load Factors
Table 2 UD T300 properties (DORNIER SYSTEM)
0.06 0.05 0.04 0.03 0.02 0.01 0.00 1 2 3
4 5 6 7 8
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Failure Events
Fig. 4. Damage pattern unzip action after 15 failure events.
Fig. 6. Failure load factors: unzip action and JIM.
28
G. Kress et al. / Composite Structures 69 (2005) 21–33 1.E+01
1.E+02
1.E+01
1.E+00 1
2
3
4
5
6
7
8
9
10
Load Factor Deviations
Load Factor Deviation
1.E+00 1
2
3
4
5
6
7
8
9
10
1.E-01 1.E-02 1.E-03 1.E-04 1.E-05 1.E-06
1.E-01
Number of Iterations
Number of Iterations
Fig. 7. Number of JIM iterations and accuracy measure (41).
Fig. 9. Number of UPCG iterations and accuracy measure (41).
1.0
0.08
Relative Computing Times
0.9
Failure Load Factors
0.07 0.06 0.05 0.04 0.03 0.02
0.8 0.7 Exact JIM UPCG Linear (JIM) Linear (UPCG)
0.6 0.5 0.4 0.3 0.2 0.1
0.01
0.0 1
0.00 1
3
5
7
9
11
13
15
17
19
21
23
2
3
4
5
6
7
8
9
10
Number of Iterations
25
Failure Events
Fig. 10. Relative computing times for the 25/25/25 problem.
Fig. 8. Failure load factors: unzip action and UPCG.
5.1.2. Unzip action and sensitivity formula method SFM The results of the sensitivity formula method are the same as that of the first iteration step presented in Section 5.1.1. Fig. 4 obviates that the damage-induced changes of the structural behavior are much too large to be captured accurately by a linear exploitation of sensitivity. 5.1.3. Unzip action and UPCG method The Untransformed Preconditioned Conjugate Gradient (UPCG) algorithm given in Section 3.4 obtains the failure–load curves shown in Fig. 8. They represent the predicted failure loads for iteration numbers ranging from one to five. The failure loads converge rapidly against those calculated by DS, shown in Fig. 5. The curves obtained with iteration numbers higher than four cannot be distinguished from the reference curve in Fig. 8. Fig. 9 obtains an exponential decay of the deviations from the reference curve with increasing number of iterations. The figure indicates that every two iterations reduce the error by one order of magnitude. 5.1.4. Unzip action and computing times The computing times indicated in Fig. 10 refer to the mesh size 25-25-25 indicated in Fig. 4 and the unzip
action over 25 elements. They are normalized with respect to the time needed by the direct solution method. Not surprisingly, JIM requires less computational effort per iteration than UPCG. Extrapolating the curves, the computing time required by the exact solver is exceeded by UPCG after about 14 iterations and by JIM after about 24 iterations. Nevertheless, UPCG is much more efficient because of its higher convergence rate: comparing Figs. 7 and 9 with each other obtains that UPCG finds better accuracy with two iterations than JIM with ten for the here considered problem. Next, the dependence of computing times on problem size parameters is investigated. Table 3 provides data in terms number of elements Nel , number of degrees-offreedom Ndgf , half-bandwidth Nhbw , computing time with direct solver tds , and computing time with UPCG tupcg . The computation times of the various mesh sizes can be estimated by using the following model. It is based on the significant time consuming processes performed during the damage progression analysis: Tgs Tbs Tsp Tfpf
complete Gauss solution process time back-substitution processing time time required for forming inner products first-ply failure processing time
G. Kress et al. / Composite Structures 69 (2005) 21–33
29
Table 3 Problem sizes and computing times Nx
Ny
Nr
Nel
Ndgf
Nhbw
tds
tupcg
15 15 15 15 15 15 25 25 25 25 25 25 50 50 50 50 50 50
15 15 25 25 50 50 15 15 25 25 50 50 15 15 25 25 50 50
25 50 25 50 25 50 25 50 25 50 25 50 25 50 25 50 25 50
975 1725 1625 2875 3250 5750 1125 1875 1875 3125 3750 6250 1500 2250 2500 3750 5000 7500
2092 3642 3432 5982 6782 11832 2412 3962 3952 6502 7802 12852 3212 4762 5252 7802 10352 15402
66 66 106 106 206 206 66 66 106 106 206 206 66 66 106 106 206 206
611 1031 1510 2570 9170 16500 671 1030 1690 2771 10300 17640 831 1179 2150 3224 13200 20500
484 708 896 1414 2760 4743 532 730 1000 1550 3130 5082 642 848 1350 1795 4040 5970
The number of floating-point operations of the Gauss triangulation depends linearly on the number of degreesof-freedom Ndgf and quadratically on the half-bandwidth Nhbw : 2 þ 32Nhbw 1 : ð42Þ tgs ¼ aNdgf 2Nhbw
tds ¼ Nde ðtgs þ tfpf Þ:
The number of floating-point operations for the backsubstitution process depends linearly on the number Ndgf and linearly on the half-bandwidth Nhbw :
The model parameter values achieving the fit illustrated in Fig. 11 are listed in Table 4. The fitting of the factors reveals that the scalar products performed by the UPCG do not significantly contribute to the computation times. Fig. 11 shows the measured times obtained with both solution methods. The data points are represented with the open symbols. It also shows the predicted times for both methods and these points are indicated with the solid symbols.
ð43Þ
The number of floating-point operations for calculating the scalar products required in the UPCG algorithm depends linearly on the number of degrees-of-freedom Ndgf and the number of iterations invested for the convergence of the displacement fields: tsp ¼ cNdgf Niter :
tupcg ¼ tgs þ ðNde 1Þðtbs þ tsp þ tfpf Þ þ tfpf :
ð44Þ
The number of operations involved in the first-ply-failure analysis to identify the next damage event and the associated critical load depends linear on the number of finite elements Nel : tfpf ¼ dNel :
When the here proposed method is used where the displacement fields are approximated by UPCG, the computation time is estimated at
ð45Þ
The computation time for post processing, i.e. calculation of stress fields from the displacement fields, scales also with the number of finite elements and may therefore be considered a part of tfpf . The damage progression simulation is an iteration over the specified number Nde of damage events. Each step consists of two stages: solving the system of equations and identifying a new damage event. The total computation time required when the displacements fields are calculated by the direct solution method is now estimated at
ð47Þ
2.50E+04
Direct Solver Computing Times
tbs ¼ bNdgf Nhbw :
ð46Þ
2.00E+04
DS Time Data DS Time Model UPCG Time Data UPCG Time Model
1.50E+04
1.00E+04
5.00E+03
0.00E+00 1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16 17 18 19 20
Ordered Sample Problems
Fig. 11. Model fit to computing times.
Table 4 Fitting parameters for fit Fig. 11 a
b 6
0.58 · 10
4
0.49 · 10
c
d
0.0
0.9 · 102
30
G. Kress et al. / Composite Structures 69 (2005) 21–33
For presentation in Fig. 11, the data of Table 3 have been rearranged in ascending order of the times tds . The fitted models are used for an assessment of the model-size dependent computation time savings which are possible by using the UPCG approach (see Fig. 12). The time savings are given in terms of the ratio between tds and tupcg . Let us use a problem-size parameter N which equals the number of unknowns Ndgf ¼ N . Having the currently considered model with linear plane elements in mind, the number of elements may roughly be set equal to the number of nodal points and, since these each carry two degrees-of-freedom, we have Nel ¼ 0:5N . The half-bandwidth is assumed to grow by the square pffiffiffiffi root of the problem-size parameter, or Nhbw ¼ N , which is in approximate agreement with the data in Table 3. The number of damage events necessary to give complete information, on the load carrying capability or the damage tolerance of a structure, is somewhat arbitrarily assumed to grow also by pthe ffiffiffiffi square root of the problem size parameter, Nde ¼ N . Under these assumptions the models obtain that pffiffiffiffi tupcg 1 ð2a þ bÞN þ ð1:5a þ 0:5dÞ N a pffiffiffiffi : ð48Þ ¼ pffiffiffiffi tds N 2aN þ 1:5a N þ 0:5d a The computing times ratio is plotted in Fig. 12 over a range 103 6 N 6 107 . The computing times using UPCG may be larger than those using DS for rather small problem sizes. At very large problem sizes N the here proposed damage progression analysis with p UPCG gains efficiency over ffiffiffiffi that using DS bypaffiffiffiffifactor N . This remains true when, instead of Nde ¼ N , Nde ¼ N is assumed.
event is fiber failure found in the element on the corner at the edge of the hole. The damage pattern shown in Fig. 13 is obtained with the direct as well as with the preconditioned conjugate gradient solution methods. Instead of using a fixed number of iterations, the UPCG algorithm is terminated as soon as qnew < 2 q0 applies. This has the advantage that only as many iterations are performed as needed for a necessary accuracy. The number of iterations for each displacement-field update tends to increase with increasing accuracy. The five curves in Fig. 14 show, in ascending order, the number of iterations needed to achieve the terminating -values ¼ 101 , ¼ 102 , ¼ 103 , ¼ 104 , and ¼ 105 . Fig. 14 also reveals that numbers of iterations increase with increasing number of damage events which may be due to decreasing quality of the preconditioner which remains, symbolically, the inverse of the stiffness matrix related to the undamaged structure. The damage pattern is already correctly obtained with an -value as low as ¼ 101 . Fig. 15 shows the load factors after DS and UPCG for an -value range 101 P P 105 . In contrast to unzip-action, the load curve does not monotonically decrease. UPCG approximates the loads quite well and the mean square errors shown in Fig. 16 indicate that -values smaller than about ¼ 103 does not further improve them.
5.2. Complex laminate damage patterns 5.2.1. Unidirectional laminate [0] Early damage pattern development in this laminate is dominated by matrix shear failure roughly along the tangents at the hole. This is due to a combination of high shearing stresses along the tangents and low matrix shear strength of the material. Only the last damage
Fig. 13. Damage in [0]s -laminate after 25 failures.
25
1.0 1.0E+03
1.0E+04
1.0E+05
1.0E+06
0.1
1.0E+07
PCG Iterations per Update
PCG over DS solution times
10.0
20
15
10
5
0
1
0.0
Number of Degrees of Freedom
Fig. 12. Predicted computing times ratio tupcg =tds .
3
5
7
9
11
13
15
17
19
21
23
Failure Events
Fig. 14. Numbers of iterations to achieve residual norms.
25
G. Kress et al. / Composite Structures 69 (2005) 21–33 0.10
Failure Load Factors
0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 1
3
5
7
9
11
13
15
17
19
21
23
25
Failure Events
Fig. 15. UPCG load factors during [0] damage simulation.
1.E+00 -1
-2
-3
-4
-5
Load Factor Deviations
1.E-01 1.E-02 1.E-03 1.E-04 1.E-05 1.E-06 1.E-07
31
The two damage patterns are similar but not identical. Fig. 18 reveals that both solution methods calculate almost identical values for the first ten critical loads but afterwards they find different values. The damage development seems to diverge at this point. It is believed that the damage development is not robust in the sense that, when several critically stressed layers exist, small differences in the previously calculated stiffness matrices may play a deciding role in which of the critically stressed elements will be recognized to fail first. The outcome of such decision will naturally alter the subsequent damage development. This has nothing to do with the accuracy of the iterative solution method but rather with other details of the code here used for the sample calculations. 5.2.3. Quasi-isotropic laminate [0, 45,90]s The striking difference of the damage pattern evolved in Fig. 19 to the previous ones is that the region with matrix failure is more confined towards the corner formed by the hole and the symmetry line. This is due to the arresting effect of the fiber reinforcement which is effective in all directions. After much matrix damage is accumulated, fiber fracture occurs and threatens the integrity of the laminate similar to the Unzip Action. The
1.E-08
Termination criterion value epsilon
0.10 0.09
5.2.2. Cross-ply laminate [0,90]s Multi-layered laminates can suffer several damage events within one finite element. The number of damage events leading to the damage pattern of the cross-ply laminate shown in Fig. 17 is one hundred. Similarly to the longitudinally reinforced unidirectional laminate, early damage is dominated by matrix failure. The failures in the 0-layer are caused by shear und those in the 90-layer are caused by normal stresses transverse to the fibers. At the end, fiber fracture appears.
Fig. 17. Damage in [0,90]s -laminate after 100 failures.
Failure Load Factors
Fig. 16. UPCG failure load deviations (41) versus -values.
0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 1
3
5
7
9
11
13
15
17
19
21
23
25
Failure Events
Fig. 18. DS and UPCG load factors during [0,90]s damage simulations.
Fig. 19. Damage in [0,±45,90]s -laminate after 100 failures.
32
G. Kress et al. / Composite Structures 69 (2005) 21–33
load curves plotted in Fig. 20 reveal that the first 24 damage events are found identical by both solution methods. Again, UPCG converges quickly and a lowprecision termination with ¼ 101 already obtains accurate results, and the divergence between the damage patterns obtained by both methods is suspected to be due to some numerical deviation of the stiffness matrix update rather than lack of accuracy of UPCG. The quick convergence, even at weak accuracy demands on UPCG termination, is due to the fact that displacement field perturbations caused by matrix failure in layers of the quasi-isotropic laminate must be small because the fiber reinforcement restrains equally ‘‘in all’’ directions (hence quasi-isotropic). Therefore, it is interesting to consider the question of whether the displacements changes encountered here are small enough to fall within the application regime of the sensitivity-formula-inspired method (SF) derived in Section 3.3. It has been found that the critical loads evaluated by SF cannot be distinguished by the eye from the curves obtained by UPCG plotted in Fig. 20. Fig. 21 illustrates the computing times when using direct, the UPCG with only one iteration, and the SF solution techniques, on two different model sizes. The smaller
problem and the larger problems have the mesh sizes Nx ; Ny ; Nr ¼ 25; 25; 25 and Nx ; Ny ; Nr ¼ 50; 50; 50, respectively, and the corresponding numerical sizes can be looked up in Table 3. The possible time savings of SF with respect to UPCG restricted to only one iteration appear small which may be due to the similar preconditioning techniques. The UPCG appears to give slightly better accuracy even after the first iteration so that the case for SF is relatively weak. However, both methods save much time if compared with the direct solution and Fig. 21 illustrates how quickly the effect increases with increasing problem size. 5.2.4. Laminate [ 45]s The ½ 45 s -fiber reinforcement allows the matrix damage to spread as shown in Fig. 22. The damage zone extends from the notch to the free edge, reminding of the distinct ±45 cracks of the damage pattern shown in Fig. 1. It must be noted that the quarter model with its symmetry assumption prevents the anti-symmetric shape appearing in the figure. The two solution techniques produce similar damage patterns although the load plots in Fig. 23 reveal that the damage paths of the two solution methods separate quite early.
0.12
Failure Load Factors
0.10 0.08
0.06 0.04 0.02 0.00 1
3
5
7
9
11
13
15
17
19
21
23
25
Failure Events
Fig. 22. Damage in [±45]s -laminate after 200 failures.
Fig. 20. Load factors of [0,±45,90]s -laminate.
50000
0.16
45000 0.14
0.12
35000 30000
Exact Solution Conjugate Gradient Method "Sensitivity Formula"
25000 20000 15000 10000
Failure Load
Computing Time
40000
0.10
0.08 0.06 0.04
5000 0.02
0 25
50
Methods at Mesh Sizes
0.00 1
2
3
4
5
6
7
8
9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
Number of Damage Event
Fig. 21. Computing times with DS, UPCG, and SF at two problem sizes.
Fig. 23. Load factors of [±45]s -laminate.
G. Kress et al. / Composite Structures 69 (2005) 21–33
6. Conclusion Approximative and iterative solution methods have been evaluated in context with damage progression simulation. The efficiency of the approximate solution techniques is enhanced by rearranging the system of equations so that it represents damage stages of the structure and uses as primary unknowns the damageinduced changes of the displacement field of the virgin structure. The solution methods have been assessed by performing damage progression simulations. The damage modelling, as such rather crude, suffices to serve the purpose. Several unidirectional as well as multidirectional laminates have been considered. A one-step approximation method inspired by the sensitivity formula for finite element models, often used in the context of gradient-based design optimization methods, gives good results when simulating the effects of matrix cracking in multidirectional laminates where the fiber arrangements prevent damage-induced stiffness variations to become large. The considered iterative solution methods, Untransformed Preconditioned Conjugate Gradient Method (UPCG) and a method inspired by the Jacobi Method (JIM), are applicable to all types of laminate. They can reproduce any damage pattern resulting from using exact solution of the state equations. Not surprisingly, UPCG converges much quicker than JIM. Thanks to the preconditioning and a good starting vector, very good accuracy is obtained with a very small number of iterations if compared to the number of unknowns. The computing times, resulting from using a frontal solver and from using UPCG with five iterations, have been modelled and the free parameters adjusted to fit computing times obtained on a laptop computer. The impact of solution technique on the numerical effort of
33
damage progression simulation has been assessed. Under some assumptions on the growth characteristics of typical problems, it has been found that UPCG reduces thepoverall numerical-effort growth scale by a ffiffiffiffi factor of N if compared to using direct solution techniques. The most striking result is that UPCG, embedded in damage progression simulation and enhanced by the here suggested precondition technique, performs significantly better than direct solution methods although the encountered matrices have a size and a banded structure where direct solution methods could otherwise not be challenged by any iterative method.
Acknowledgements Efficient solution techniques in the presence of small structural stiffness matrix variations have suggested themselves in the course of our research on design optimization with Evolutionary Algorithms, so that grant no. 2100-066879.01/1 of the Swiss National Science Foundation is again gratefully acknowledged.
References [1] Kress GR, Stinchcomb WW 1985. Fatigue response of notched graphite/epoxy laminates. In: Vinson JR, Taya M, editors. Recent Advances in Composites in the United States and Japan, ASTM STP 864, American Society for Testing and Materials, Philadelphia. p. 173–96. [2] Hashin Z. Failure criteria for unidirectional fiber composites. J Compos Mater 1980;7:329–34. [3] Kress G. Free-edge influence on CFRP-laminate strength. Int J Damage Mech 1994;3:192–211. [4] Shewchuk JR. An introduction to the gradient method without the agonizing pain. Pittsburgh PA, 15213: School of Computer Science, Carnegie Mellon University; 1994.