- Email: [email protected]
Pre- and post-degradation analysis of composite materials with different moduli in tension and compression
134
S. Idelsohn et al., Degradation analysis of composite materials
pression. In fact, one of the important characteristics of composite materials is that they often exhibit different tensile and compressive moduli. (4) A good stiffness evaluation of the faminate must be achieved even when geometric nonlinearities become important. (5) The total number of degrees of freedom must be independent of the number of layers in order to achieve computational efficiency. This finite element consideration must be complemented with a correct nonlinear solution technique in order to eliminate numerical errors due to the geometric nonlinearities and the integration of the constitutive equations.
2. Finite element model It is generahy recognized that the classical thin plate theory is not su~~iently accurate for most composite laminates in which transverse shear defo~atio~ effects are important 131. In fact, for anisotropic plates and shells, the distorsion of the deformed normal due to transverse shear is dependent not only on the laminate thickness, but also on the orientation and degree of orthotropy of the individual layers. To establish a rational laminate analysis technique which can be applied to conventional composite st~cture, the effects of transverse shear deformation should be included. If one wants to derive a finite element that can be easily generalized to include geometric nonlinearities, the most straightforward solution consists in making a displacement assumption inside the element [4, 51. On the other hand, the equilibrium equations at the boundaries, if expressed in terms of surface tractions, remain always linear, and consequently the connectors between the elements can be arbitrarily chosen as being either of a displacement or an equilibrium nature. Different solutions may be adopted to analyse a multilayer laminate using displacement type elements. Firstly, each layer of the laminate can be separately discretized using a kinemati~ally admissible element. An example is shown in Fig. 1 for a three layer laminate, using the classical 8 node isoparametric axisymmetric displacement element as a basic element. In this case each layer has its own displacement field and severe cross-sectional warping is allowed. It
Fig. 1. One efement per Iayer discretization.
S. Idelsohn et al., Degradation analysis
of composite materials
135
Fig. 2. Muitilayer element,
leads to a very accurate solution but the procedure becomes rapidly intractable for practical structures due to the excessive number of degrees of freedom which result. Secondly, one can postulate a unique displacement field for the whole laminate. In this case, the stiffness matrix is obtained by summing the contribution of each layer using the same strain field for the whole element. The total number of degrees of freedom of the element is then independent of the number of layers and the corresponding solution remains quite economic (Fig. 2). However, it must be emphasized that the strain field is the same for the whole element and cross sectional warping cannot be represented with high accuracy. This last solution was adopted in the present element. The basic model is the classical 8 node isoparamet~c element. It can be used in plane stress, plane strain and axisymmetric structures. Orthotropic stress-strain relations are introduced for each layer. Several possibilities are presented and this is one of the main originalities of the paper. These stress-strain relations may be linear or nonlinear before degradation of the material matrix and an original model is presented for composite materials with different moduli in tension and compression. In addition, three different stress-strain relations are adopted for the post-degradation behavior.
3. The strain-stress
relation before degradation
Three types of material behavior are considered in the element before the matrix degradation: (a) linear elastic material; (b) nonlinear material; (c) material with different moduli in tension and compression. 3.1. Linear elastic material The strain-stress
relation in each layer for a given material is written as
with Et
=
[&II,
,522,
-533,
Y12,
Y-23,
Y1317
d = [UII,a;?z,
m3r
712,
.i-23, 7131
.
S. Idelsoh
136 HI is
the 6 x
6 s~~~~~~~c
et al., Degradation
analysis
of composite materials
compliance matrix and its inverse C = HI’ is the Hooke matrix.
The experimental stress-strain relations of a composite material show that when one of the principal stresses is nearly coincident with the fiber direction, the behavior is practically linear, However, for another stress state the matrix material pIays an important role and the behavior becomes nonlinear. An example is given in Fig. 3. Several models have been developed in order to describe mathematically the nonlinear behavior of the composite material [4]. An incremental model due to Sandhu [6] is adopted here which takes into account the multiaxial strain state. The mechanical properties are functions of equivalent strains defined as dG = ds&d&
~&L/d~1z)
= de&f
- ~~@m/d~&-
dE% = ddl-
dy; = dyij
~&hld~xJ
- ~~3~d~~~~d~~~~~ I r&d&dcz)) - v&h/dtr33))
, ,
for if j,
where dEii and dy+ are the strain variations in the orthotropic strain variations are related to the orthotropic stress variations
directions. These equivalent bi-univocally by the tangent
Go 50
Jo 20 10 0
3
ride
“25 5
_E
f
:
7
y,=amstant=.25
c 5 4 3 2 I .l
.2
-3
Fig. 3. Stress-strain
.L
relation for a fiber g&x-epoxy
composite (3M X p251S).
I.25
5
S. Kdelsohn et al., Degradation analysis of composite tnaterials
mechanical properties
of the mat&a1 E: and Gi, (not summed),
d&:t = dcFjJE::f&) , Combining
137
dr$ = drii~~~j(~~~.
(2) and (3) provides the incremental
(31
strain-stress
relations for a layer
dc:=H’du with the tangent compliance
(4) matrix W’ defined as 1
0
l/G&
_l
Cubic spline functions have been used to represent the experimental stress-strain curves from which the tangent moduli can be determined. Analytical results determined by the present theory correlate very well with experimental data [6].
An important characteristic of composite materials in that they ~equently exhibit different stiffnesses under tensile loading and under compressive loading. Some materials have tension moduli greater than compression mod&, others have the opposite behavior [7]. Although this behavior does not admit at the present time a general physical explanation, it must be taken into consideration in a finite element analysis in order to obtain realistic results. Several models have been presented during the last few years [8]. In all cases the uniaxial stress-strain relations are approximated by a bilinear representation (Fig. 4) and the different models differ in the manner by which the multiaxial constitutive relations are devised.
Fig. 4, Bilinear uniaxiaf stress-strain relations.
S. Idelsohn
138
et al., Degradation
analysis
of composite materials
Although a discontinuity of the slope at the origin is not physically meaningful, the model is a good approximation for all the materials that are normally represented by otherwise linear stress-strain relations. It should be noted that even if bilinear uniaxial stress-strain relations are used for each stress component, the multi-axial state is n-linear (piecewise linear) and then special care must be taken in the numerical integration of the constitutive relations. For ‘multimodulus’ materials the structure of the matrix Ii defined by (1) is not well defined when some stress components are in tension and others are in compression. Ambartsumyan [S] developed a model for multimodulus orthotropic material under plane stress. He gives the strain-stress relation in principal stress coordinates as Hjj =
Hij H;
if ull>O and cr22)0, if rrll
(6)
where subscript t or c denotes either the tensile matrix Hi, or the compressive matrix. The advantage of working in the principal stress axes is that the shear stresses are zero and thus all the corresponding terms in the matrix H need not be taken into account. Nevertheless as the principal strain directions are not aligned with the principal stress, ylz is not zero and the terms H& are present in this fo~ulation. An important deficiency of this model is that the symmetry of the matrix H is not preserved. Even if symmetry is enforced in (6) putting Hi2 = H&, the matrix will not be symmetric in any coordinate other than the principal stress directions. To overcome this difficulty, Jones [7] proposes to use weighting factors Ki for the terms HI2 and H,, of the matrix KZ, + &HE;? H21 = H12 = KlHC12+ K2Hi2
for cn>O, for all
ouCO, gu>O,
(7)
where
In the present paper another criterion for obtaining the strain-stress relation for a multimodulus orthotropic material is presented. Instead of working in the principa1 stress components, the matrix C is built in the orthotropic directions. In these axes the strain-stress matrix becomes
Then only 9 terms are di~erent from zero. Moreover, using the property that the shear moduli in the o~botrop~c frictions are the same in tension and compression [7], only 6 terms of the matrix C must be determined for a rn~lt~udu~us material. For the first 3 x 3 components, the foffowing criterion is proposed.
Instead of using the arbiter values Ki as proposed by Jones, it is preferred here to use constants Qij and pij that are easily changed for di~erent materials, In particular, if in each layer there is only one fiber direction coincident with the orthotropic direction 1, the authors recommend the following coefficients:
In this way, the preponderant elastic modulus of the fibers is taken into account in the coupling terms. This factor can be adjusted fur each material. Fur bi~direct~u~a~ laminates with the same fiber density in directions 1 and 2, the following coefficients may be used: a21 = OS, pzl = 0.5, a3%= 0.8, /3s1= 0.2, cxi32 = 0.8, & = 0.2 and for mat laminates all coefficients must be taken equal to 0.5. For bilinear material, in addition to G12, G13 and Gs, the values of El, E2, E3, v12, 293, v23in tension and in compression are introduced for each element layer. For nonlinear materials, spline functions are used for each stress component having a different shape for positive and for negative vafues.
Failure criteria must be i~truduced to model the n~n~iuear behavior due to partial d~~radati~~ and to represent pin-by”ply progressive failure in the real laminates. Therefore, one must specify for each layer a failure criterion and a behavior after degradation. Five different failure criteria which have been proposed in the literature [4, 9] have been adopted in the development of the element. In all cases the average stress of the layer is used to check whether the criterion is verified or not.
The Ts~~-~i~~ cri~~~i~~[lo] states that in the orthotropic degradation is reached when
directions of a layer, the matrix
140
S. Idelsohn
et al., Degradation
analysis
of composite materials
with f(a) = F(u~ - ~~3)’+ G(a, - al)’ + H(al - ~2)~+ H(al - CQ)*+ 2L& + 2M7:1+ 2N& and 2F=X;*+X;*-X;*,
2L= R;*,
2G=X;*+X;*-X;*,
2M = R;*,
2H=X;*+X;*-X;*,
2N= R;*,
where Xi are the longitudinal are the shear strengths.
tensile strength,
(12)
and the transverse
Tsai- Wu criterion The Tsai- Wu criterion [ll] takes into account the different
compression.
tensile strengths, while Ri
stress limits in tension and in
The function f(a) in (11) is replaced by
f(a)
= Eui + l+iaj
i, j = 1) . e e 36 7
,
(13)
where E = 1/x:
- 1/x;
)
&
I;l=o,
for i = 1,2,3,
Ei = l/X:X; =
for i = 4,5,6,
l/R:
for i#j; for if j;
i,j= 1,2,3, i, j = 4,5, 6.
(14) (15)
In these formulas, X: and XT represent the longitudinal tensile strength and compression strength of the i-component respectively, Ri the shear strength and Ui the shear strength in a test at 45” with respect to the orthotropy direction in the planes 1-2, l-3 and 2-3. Maximum stress criterion The maximum stress criterion [2] states that the matrix degradation
the orthotropic Uii
>
This criterion, criterion.
components Xi
or
is attained when one of the tensile strength or compression strength;
of stress overcomes
Uii < Xy,
/Tij( > Ri e
(16)
expressed as a function of the strain tensor, becomes the maximum strain
Sandhu criterion The Sandhu criterion [6] is based on the concept
transverse, and shear loadings are independent total strain energies. It can be written as aij
d&ij Kij = 1 >/ I
that strain energies under longitudinal, parameters. The failure criterion is based upon
(17)
where the & are the failure principal strain energies under uniaxiaf tension and compression, This faihre criterion is especiaily adapted to the nonlinear Sandhu model described in Section 3.2. This model uses the equivalent strains defined in (7). In this way it seems better to use also the equivalent strain in the failure criterion. The modified Sandhu criterion becomes
After the matrix degradation of a gayer, the transverse and shear stiffness decrease to a small vafue. Only the elastic modulus in the fiber direction is preserved and all the coupling terms in the stress-strain relations are relaxed. However, in filament wound composite materials the degraded elastic and shear moduli are not negligible. For this reason, after the matrix degradation, the compliance matrix has been changed into
where Ed and Gd are the transverse and shear moduli after degradation obtained by experimental tests. Three possibilities have been retained in this finite element. (a) The transverse and shear stresses in the degraded Iayer of the composite are assumed to be maintained at their value. That means
where t+* and E* are the stress and strain vectors at degradation. (b) The transverse and shear stresses in a degraded layer are relaxed: a =
I&l&
(21)
(c) The transverse and shear stresses in the degraded layer decrease Iinearly from the failure value to zero in such a way that the failure principal strain energies (i.e. Kij) keep the Same value (cf. Fig. 5 and [4]). When the l~n~tu~nal stress reaches a critical value (different for each failure criterion) the fibers break, aII the stresses are then relaxed and are carried by the neighbuuring layers. These total failure criteria are:
142
S. Idelsohn
et al., Degradation
analysis of composite
materials
MATRIX _LEGRAMTION ------._______________
Fig. 5. Third model of post-degradation
(a) Tsai-Hill
: 2>H 1 *
(au/X)
(b) Tsai-Wu Flail
behavior.
(22)
: + i+T:,
(c) Maximum
2 1.
(23)
stress or maximum
crI12XI
or
strain:
El1 2 Xf .
(24)
(d) Sandhu: (&I( 2 Ed (uniaxial)
6. Newton-Raphson
.
(25)
method in composite materials
The Newton-Raphson method is used to solve the nonlinear problem. Two special considerations must be taken when implementing this algorithm: the correct integration of the constitutive equations in the presence of multimodulus material or during the matrix degradation and the possibility of numerical degradations during geometrically nonlinear problems. 6.1. Geometrically The structural
linear problems discretization
of the virtual
statement
yields the nonlinear
problem (26)
r(q) = f(4) - g = 0 where g are the generalized external loads and f(q) the internal adding the internal force contribution of each element B’u dv I Li
forces that can be obtained
by
(27)
S. Idelsohn et ul., Degradation analysis of composite materials
in which
B
is the matrix connecting
the displacement
parameter
143
with the strain tensor
A.e=BAq.
e=Bq, The Newton-Raphson
(28)
method consists in adopting an estimate q to the solution
q=B+Aq,
(2%
a=$++~,
(30)
and assuming that in the neighbourhood
of q the linear solution
Ao, = C~~,~)AE is a good approximation. r(q)
A presumably
(31) From (27)-(32), (27) becomes
(B’CB)
=
dv Aq f
better approximation
KAq=g-1
v
I0
B’C(6,2)B
Thus, Newton’s method iteratively by (33) qk+l=
qk
-
B’S dv - g.
(32)
can be obtained by putting (32) to zero, which yields
B’&dv
where the tangent stiffness matrix K =
I”
(33) K
is defined as
dv .
takes an initial approximation
(34) go
to q and attempts
to improve
it
(35)
Kh’k.
Then the new strain and stress inurement values are calculated by
A-%= B(qic+l-
qk)
(36)
and (37) and a new tangent matrix is evaluated. The procedure is stopped when the residual vector (r) becomes less than a desired precision. It is in the integration of the constitutive relations (37) that the composite material has a particular behavior.
144 6.1.1.
S. Idelsohn
et al., Degradation
analysis
of composite
materials
In N-linear problems
If the matrix degradation is reached in any point of the structure, or if a change in the stress sign occurs for a multimodulus composite material, the stress increment is written
where N is the number of stress-strain relation changes during the iteration. To evaluate Aci, a factor ri is calculated at each step so that ri As = hei is the part of the increment necessary to reach a discontinuity point (degradation or stress sign change). The stresses are evaluated at this point and a new strain-stress matrix is calculated. The process is repeated until the total strain increment is obtained. 6.1.2. Nonlinear material problems
In this case the total strain increment is subdivided into m sufficiently small parts and the total stress increment is evaluated by Euler’s explicit integration. In this way it is easy to verify when the matrix degradation or a stress sign change occurs. 6.2. Geometrically nonlinear problems The use of the Newton-Raphson’s method to solve geometrically nonlinear problems involving composite materials may introduce some numerical difficulties. In fact, when geometrical nonlinearities occur, the displacement fields are over-evaluated at the first iterations. Thus, it is possible that a failure criterion is verified in a layer before the real stress reaches the collapse value. Owing to the irreversiblity of the phenomena the algorithm may converge to a non physical solution. In order to improve this situation the nonlinear geometrical iterations are carried out separately from the nonlinear material iterations. In this way the problem is first considered as a linear material problem and the geometrically nonlinear iterations are performed. When the equilibrium is attained the correct numerical integration of the constitutive relation is achieved and a new step with only geometrically nonlinear iterations is performed.
7. Numerical results 7.1. Axisymmetric cylinder
To test the behavior of the multilayer finite element developed here and to compare the different failure criteria used, a simple four layer cylinder submitted to a shear load at the end was examined [9]. The structure is modeled by two elements of two layers in the radial direction and only one element in the axial direction (Fig. 6). The fiber direction in each layer is of 90”, O”,45” and 75” with respect to the axial direction starting from the internal layer to the external one. The
S. Idelsohn et al., Degradation analysis
. 0 C. 0
of composite materials
145
TSAI-HILL TSAI-vvu SANDHU sigma max. “c
Fig. 6. Axisymmetric cylinder; maximal radial displacement for the different faifure criteria.
mechanical properties
of the glass epoxy composite material before degradation
E,=5500hb, zf12=
v23
=
E2=Eg= VI3
=
1800hb,
G,, = G13= GZ3= 900 hb ,
0.25 )
and after the matrix failure of a layer, the transverse reduced to Ed= lhb,
are
and shear moduli are considerably
G,=0.3hb.
The longitudinal strength is 105 hb; the transverse and shear strengths are only 7 hb and 5 hb respectively, and for the Tsai-Wu criterion, the shear strength at 45” is 7.5 hb. The axisymmetric shear load applied at the end of the cylinder was incremented by a first step of 1.5 hb and then by 7 steps of 0.5 hb each. The maximal radial displacement (point c) for the different failure criteria used is shown in Fig. 6. It can be noted that with a shear load of 1.5 hb no degradation is reached for any criteria and in the next increment in ail cases one layer is degraded. However, between a shear load of 2 hb and 4 hb the failure of the 2nd and 3rd layer is reached at different moments with each failure criterion used and then different solutions are obtained. When the maximal shear load (5 hb) is applied, neither layer of the cylinder is completely broken. 7.2. Spherical ring with different behavior in tension and compression The second example considered is that of a clamped spherical ring subjected to an initial compression load and then successive incremental rotations of the ring are carried out to introduce centrifugal load. Its geometric characteristics are summarized in Fig. 7. This spherical ring is built of 5 Iayers from a composite material made up of Kevlar 49 fibers and epoxy matrix. Their mechanical properties are described in Table 1. The spherical shell thickness is h = 2 mm and there is an initial layer of 0.4 h with meridia1 fiber direction. Above and below this layer there are respectively 2 layers of 0.15 h each with fiber directions oriented at +45” and -45” with respect to the internal layer.
146
S. Idelsohn
et al., Degradation
analysis
of composite materials
Fig. 18. Spherical ring discretization.
The structure is modeled by using 7 axisymmetric finite elements of 5 layers each. The density of the external element is considered 10 times larger than the composite density. This is to take into account some concentrated mass placed there. The maximum stress criterion is used assuming that after degradation transverse and shear stress are maintained at their degradation values. The stress-strain relations are considered bilinear after degradation and the geometrical nonlinear strain is neglected. To check the importance of the multimodulus behavior in structures where the sign of the stress state varied, three different tests are performed (1) With only tensile modulus. (2) With only compressive modulus. (3) With the multimodulus (model developed in Section 3.3). Fig. 8 shows the maximal displacement at the boundary of the ring for the three tests. In the first increment only the initial compression load exists. Then the rotations are incremented to 31.6, 400.0, 700, 1000, 1200 and 1400 rad/sec and produce then a tension state in the structure. Table 1 Material properties of the composite
layer
Longitudinal tensile modulus E: Transverse tensile modulus E; Longitudinal compressive modulus E’, Transverse compressive modulus ES Shear modulus G12 Tension Poisson’s ratio Compression Poisson’s ratio Density
8.500 560 6.000 700 210
kg/mm’ kg/mm’ kg/mm’ kg/mm’ kg/mm’
0.34 0.3 1.56 gm/cc
Longitudinal tensile strength Transverse tensile strength Longitudinal compressive strength Transverse compressive strength Shear strength
141 2.8 2.8 14.1 4.5
kg/mm’ kg/mm’ kg/mm2 kg/mm’ kg/mm’
S. Idelsohn et al., Degradation analysis of composite materials
___o_._
CCMPRESSIVE
-o-
TENSILE
+
MULTlMOWUJS
147
MODUWS
MODULUS
Fig. 8. Spherical
ring; maximal
displacement.
At a rotational velocity of 1000 rad/sec, 2 layers reach the degradation in the multimodulus material instead of 1 layer with the traction test and 5 layers in the compression one. At the final rotational velocity of 1400 rad/sec 12 layers reach the degradation instead of 18 layers using the compression modulus. The maximal displacement obtained in this case is 34% above the value corresponding to the multimodulus theory. In fact, a small difference between the tensile and compressive moduli introduces for a given load a difference in the number of degraded layers. When new increments are performed this difference produces a new break-up and so on. This behavior explains the large differences encountered between a multimodulus analysis and a single modulus one in the presence of matrix degradation. 7.3. Geometrically nonlinear cantilever beam To test the element behavior in geometrically nonlinear problems a very simple cantilever beam is discretized by a plane stress element made up of four layers each oriented at O”, 90”, 45” and -45”. Table 2 Geometrically
nonlinear
cantilever
beam
Number of geometrically nonlinear iterations Increment 1 2 3
Tsai-Hill 6 5 4+4
Sandhu 6+1 5+1 5+1+2+1
Number of material nonlinear iterations Tsai-Hill 1 1 2
Sandhu 2 2 4
Matrix degradation non non yes
148
S. Idelsohn
.75-
et al., Degradation
analysis
SIRESS-SlRAIN RELATIONSFORWWU
Fig. 9. Cantilever
of composite materials
CRIlERK)N
beam.
All the mechanical properties before and after degradation are the same as those used in the axisymmetric cylinder described in Section 7.1. Two different failure criteria are tested here: Tsai-Hill and Sandhu. In the Sandhu criterion, nonlinear stress-strain relations are used (Fig. 9). Nevertheless these nonlinearities are mild and the two failure criteria can be compared. The external load is a uniform pressure p distributed on the top of the beam. Three increments are applied to obtain respectively: 2.5, 5.0 and 7.5 hb. Table 2 summarizes the performances obtained using the two failure criteria. It can be seen that in both cases, 2 layers reach the matrix degradation at the 3rd increment. Six geometrically nonlinear iterations are necessary for the first increment and 5 for the second. The Sandhu criterion needs a second nonlinear material iteration for these 2 increments. The third increment shows clearly the uncoupled behavior between the geometrical and material iterations: when the Tsai-Hill failure criterion is used, 4 geometrically nonlinear iterations are required to obtain equilibrium, In the subsequent mechanical iteration two layers degrade and to achieve a new equilib~um, 4 geometrical iterations are again needed. A similar behavior is observed by using Sandhu’s failure criterion.
8. Conclusions Some difficulties inherent to the development of a finite element for the analysis of fiber reinforced composite materials have been presented. The proposed element has been conceived for use in a research laboratory in order to compare numerical solutions with experimental results. For these reasons a large number of possibilities and criteria are present in the model. The desire of the authors is to present a mathematical tool to aid a better unde~tanding of the complicated composite material behavior. Unfortunately, very few experimental results are available in the literature in order to confirm the analysis. Nevertheless, the examples show that the solution of a problem is
S. Idelsohn et al., Degradation analysis of composite materials
149
sensitive to the failure criteria used and that the differences between the tension and compression moduli must not be neglected when change in the stress signs occurs. Future research and numerical tests should be devised in order to adjust the coupling terms cy, /3 in the strain-stress relation of material with different moduli in tension and compression, the degraded elastic and shear moduli E, and Gd and the differences between tension and compression stiffness values. The effect of shear stresses in the multimodulus material should be explored to ascertain their influence in the coupling terms. The analysis of decreasing load and numerically decreasing stresses after the matrix degradation are actually unsolved problems.
References [l] G. Sander and C. Nyssen, Modelisation des materiaux composites dans les methodes d’elements finis, L.T.A.S. SF-90, presented at XIVth Internat. Aeronautical Congress, Paris, 1979. [2] R.M. Jones, Mechanics of Composite Materials (McGraw-Hill, New York, 1975). [3] E. Hinton, The flexural analysis of laminated composite using a parabolic isoparametric plate bending element, Internat. J. Numer. Meths. Engrg. 11 (1979) 174-179. [4] C. Nyssen, Modelisation par elements finis du comportement non lineaire de structures aerospatiales, Doctoral Thesis, Aerospace Laboratory, University of Liege, 1979. [5] C. Nyssen and P. Beckers, Finite element linear and nonlinear analysis of composite materials, Proc. Conf. Nonlinear Finite Element Analysis and ADINA, MIT, Cambridge, MA, 1979. [6] R.S. Sandhu, Ultimate strength analysis of symmetric laminates, Tech. Rept. AFFDL-TR-73-137, 1974. [7] R.M. Jones, Stress-strain relations for materials with different moduli in tension and compression, AIAA J. 15(l) (1977) 16-23. [8] R.M. Jones and D.A.R. Nelson, Material models for nonlinear deformation of graphite, AIAA J. 14(6) (1976) 709-717. [9] G. Laschet, S. Idelsohn and M. Hogge, Element multicouche pour l’analyse des r&sines armees, L.T.A.S. SF-91, Aerospace Laboratory, University of Liege, 1979. [lo] S.W. Tsai, Strength characteristics of composite materials, NASA-CR-224, NASA, 1965. [ll] S.W. Tsai and E. Wu, A general theory of strength of anisotropic materials, J. Comput. Materials 5 (1971) 58-80. [12] J.L. Walsch, J.H. Ahlberg and E.N. Wilson, Best approximation properties of the spline fit, J. Math. Mech. 2(2) (1962) 225-233.
S. Idelsohn et al., Degradation analysis of composite materials
pression. In fact, one of the important characteristics of composite materials is that they often exhibit different tensile and compressive moduli. (4) A good stiffness evaluation of the faminate must be achieved even when geometric nonlinearities become important. (5) The total number of degrees of freedom must be independent of the number of layers in order to achieve computational efficiency. This finite element consideration must be complemented with a correct nonlinear solution technique in order to eliminate numerical errors due to the geometric nonlinearities and the integration of the constitutive equations.
2. Finite element model It is generahy recognized that the classical thin plate theory is not su~~iently accurate for most composite laminates in which transverse shear defo~atio~ effects are important 131. In fact, for anisotropic plates and shells, the distorsion of the deformed normal due to transverse shear is dependent not only on the laminate thickness, but also on the orientation and degree of orthotropy of the individual layers. To establish a rational laminate analysis technique which can be applied to conventional composite st~cture, the effects of transverse shear deformation should be included. If one wants to derive a finite element that can be easily generalized to include geometric nonlinearities, the most straightforward solution consists in making a displacement assumption inside the element [4, 51. On the other hand, the equilibrium equations at the boundaries, if expressed in terms of surface tractions, remain always linear, and consequently the connectors between the elements can be arbitrarily chosen as being either of a displacement or an equilibrium nature. Different solutions may be adopted to analyse a multilayer laminate using displacement type elements. Firstly, each layer of the laminate can be separately discretized using a kinemati~ally admissible element. An example is shown in Fig. 1 for a three layer laminate, using the classical 8 node isoparametric axisymmetric displacement element as a basic element. In this case each layer has its own displacement field and severe cross-sectional warping is allowed. It
Fig. 1. One efement per Iayer discretization.
S. Idelsohn et al., Degradation analysis
of composite materials
135
Fig. 2. Muitilayer element,
leads to a very accurate solution but the procedure becomes rapidly intractable for practical structures due to the excessive number of degrees of freedom which result. Secondly, one can postulate a unique displacement field for the whole laminate. In this case, the stiffness matrix is obtained by summing the contribution of each layer using the same strain field for the whole element. The total number of degrees of freedom of the element is then independent of the number of layers and the corresponding solution remains quite economic (Fig. 2). However, it must be emphasized that the strain field is the same for the whole element and cross sectional warping cannot be represented with high accuracy. This last solution was adopted in the present element. The basic model is the classical 8 node isoparamet~c element. It can be used in plane stress, plane strain and axisymmetric structures. Orthotropic stress-strain relations are introduced for each layer. Several possibilities are presented and this is one of the main originalities of the paper. These stress-strain relations may be linear or nonlinear before degradation of the material matrix and an original model is presented for composite materials with different moduli in tension and compression. In addition, three different stress-strain relations are adopted for the post-degradation behavior.
3. The strain-stress
relation before degradation
Three types of material behavior are considered in the element before the matrix degradation: (a) linear elastic material; (b) nonlinear material; (c) material with different moduli in tension and compression. 3.1. Linear elastic material The strain-stress
relation in each layer for a given material is written as
with Et
=
[&II,
,522,
-533,
Y12,
Y-23,
Y1317
d = [UII,a;?z,
m3r
712,
.i-23, 7131
.
S. Idelsoh
136 HI is
the 6 x
6 s~~~~~~~c
et al., Degradation
analysis
of composite materials
compliance matrix and its inverse C = HI’ is the Hooke matrix.
The experimental stress-strain relations of a composite material show that when one of the principal stresses is nearly coincident with the fiber direction, the behavior is practically linear, However, for another stress state the matrix material pIays an important role and the behavior becomes nonlinear. An example is given in Fig. 3. Several models have been developed in order to describe mathematically the nonlinear behavior of the composite material [4]. An incremental model due to Sandhu [6] is adopted here which takes into account the multiaxial strain state. The mechanical properties are functions of equivalent strains defined as dG = ds&d&
~&L/d~1z)
= de&f
- ~~@m/d~&-
dE% = ddl-
dy; = dyij
~&hld~xJ
- ~~3~d~~~~d~~~~~ I r&d&dcz)) - v&h/dtr33))
, ,
for if j,
where dEii and dy+ are the strain variations in the orthotropic strain variations are related to the orthotropic stress variations
directions. These equivalent bi-univocally by the tangent
Go 50
Jo 20 10 0
3
ride
“25 5
_E
f
:
7
y,=amstant=.25
c 5 4 3 2 I .l
.2
-3
Fig. 3. Stress-strain
.L
relation for a fiber g&x-epoxy
composite (3M X p251S).
I.25
5
S. Kdelsohn et al., Degradation analysis of composite tnaterials
mechanical properties
of the mat&a1 E: and Gi, (not summed),
d&:t = dcFjJE::f&) , Combining
137
dr$ = drii~~~j(~~~.
(2) and (3) provides the incremental
(31
strain-stress
relations for a layer
dc:=H’du with the tangent compliance
(4) matrix W’ defined as 1
0
l/G&
_l
Cubic spline functions have been used to represent the experimental stress-strain curves from which the tangent moduli can be determined. Analytical results determined by the present theory correlate very well with experimental data [6].
An important characteristic of composite materials in that they ~equently exhibit different stiffnesses under tensile loading and under compressive loading. Some materials have tension moduli greater than compression mod&, others have the opposite behavior [7]. Although this behavior does not admit at the present time a general physical explanation, it must be taken into consideration in a finite element analysis in order to obtain realistic results. Several models have been presented during the last few years [8]. In all cases the uniaxial stress-strain relations are approximated by a bilinear representation (Fig. 4) and the different models differ in the manner by which the multiaxial constitutive relations are devised.
Fig. 4, Bilinear uniaxiaf stress-strain relations.
S. Idelsohn
138
et al., Degradation
analysis
of composite materials
Although a discontinuity of the slope at the origin is not physically meaningful, the model is a good approximation for all the materials that are normally represented by otherwise linear stress-strain relations. It should be noted that even if bilinear uniaxial stress-strain relations are used for each stress component, the multi-axial state is n-linear (piecewise linear) and then special care must be taken in the numerical integration of the constitutive relations. For ‘multimodulus’ materials the structure of the matrix Ii defined by (1) is not well defined when some stress components are in tension and others are in compression. Ambartsumyan [S] developed a model for multimodulus orthotropic material under plane stress. He gives the strain-stress relation in principal stress coordinates as Hjj =
Hij H;
if ull>O and cr22)0, if rrll
(6)
where subscript t or c denotes either the tensile matrix Hi, or the compressive matrix. The advantage of working in the principal stress axes is that the shear stresses are zero and thus all the corresponding terms in the matrix H need not be taken into account. Nevertheless as the principal strain directions are not aligned with the principal stress, ylz is not zero and the terms H& are present in this fo~ulation. An important deficiency of this model is that the symmetry of the matrix H is not preserved. Even if symmetry is enforced in (6) putting Hi2 = H&, the matrix will not be symmetric in any coordinate other than the principal stress directions. To overcome this difficulty, Jones [7] proposes to use weighting factors Ki for the terms HI2 and H,, of the matrix KZ, + &HE;? H21 = H12 = KlHC12+ K2Hi2
for cn>O, for all
ouCO, gu>O,
(7)
where
In the present paper another criterion for obtaining the strain-stress relation for a multimodulus orthotropic material is presented. Instead of working in the principa1 stress components, the matrix C is built in the orthotropic directions. In these axes the strain-stress matrix becomes
Then only 9 terms are di~erent from zero. Moreover, using the property that the shear moduli in the o~botrop~c frictions are the same in tension and compression [7], only 6 terms of the matrix C must be determined for a rn~lt~udu~us material. For the first 3 x 3 components, the foffowing criterion is proposed.
Instead of using the arbiter values Ki as proposed by Jones, it is preferred here to use constants Qij and pij that are easily changed for di~erent materials, In particular, if in each layer there is only one fiber direction coincident with the orthotropic direction 1, the authors recommend the following coefficients:
In this way, the preponderant elastic modulus of the fibers is taken into account in the coupling terms. This factor can be adjusted fur each material. Fur bi~direct~u~a~ laminates with the same fiber density in directions 1 and 2, the following coefficients may be used: a21 = OS, pzl = 0.5, a3%= 0.8, /3s1= 0.2, cxi32 = 0.8, & = 0.2 and for mat laminates all coefficients must be taken equal to 0.5. For bilinear material, in addition to G12, G13 and Gs, the values of El, E2, E3, v12, 293, v23in tension and in compression are introduced for each element layer. For nonlinear materials, spline functions are used for each stress component having a different shape for positive and for negative vafues.
Failure criteria must be i~truduced to model the n~n~iuear behavior due to partial d~~radati~~ and to represent pin-by”ply progressive failure in the real laminates. Therefore, one must specify for each layer a failure criterion and a behavior after degradation. Five different failure criteria which have been proposed in the literature [4, 9] have been adopted in the development of the element. In all cases the average stress of the layer is used to check whether the criterion is verified or not.
The Ts~~-~i~~ cri~~~i~~[lo] states that in the orthotropic degradation is reached when
directions of a layer, the matrix
140
S. Idelsohn
et al., Degradation
analysis
of composite materials
with f(a) = F(u~ - ~~3)’+ G(a, - al)’ + H(al - ~2)~+ H(al - CQ)*+ 2L& + 2M7:1+ 2N& and 2F=X;*+X;*-X;*,
2L= R;*,
2G=X;*+X;*-X;*,
2M = R;*,
2H=X;*+X;*-X;*,
2N= R;*,
where Xi are the longitudinal are the shear strengths.
tensile strength,
(12)
and the transverse
Tsai- Wu criterion The Tsai- Wu criterion [ll] takes into account the different
compression.
tensile strengths, while Ri
stress limits in tension and in
The function f(a) in (11) is replaced by
f(a)
= Eui + l+iaj
i, j = 1) . e e 36 7
,
(13)
where E = 1/x:
- 1/x;
)
&
I;l=o,
for i = 1,2,3,
Ei = l/X:X; =
for i = 4,5,6,
l/R:
for i#j; for if j;
i,j= 1,2,3, i, j = 4,5, 6.
(14) (15)
In these formulas, X: and XT represent the longitudinal tensile strength and compression strength of the i-component respectively, Ri the shear strength and Ui the shear strength in a test at 45” with respect to the orthotropy direction in the planes 1-2, l-3 and 2-3. Maximum stress criterion The maximum stress criterion [2] states that the matrix degradation
the orthotropic Uii
>
This criterion, criterion.
components Xi
or
is attained when one of the tensile strength or compression strength;
of stress overcomes
Uii < Xy,
/Tij( > Ri e
(16)
expressed as a function of the strain tensor, becomes the maximum strain
Sandhu criterion The Sandhu criterion [6] is based on the concept
transverse, and shear loadings are independent total strain energies. It can be written as aij
d&ij Kij = 1 >/ I
that strain energies under longitudinal, parameters. The failure criterion is based upon
(17)
where the & are the failure principal strain energies under uniaxiaf tension and compression, This faihre criterion is especiaily adapted to the nonlinear Sandhu model described in Section 3.2. This model uses the equivalent strains defined in (7). In this way it seems better to use also the equivalent strain in the failure criterion. The modified Sandhu criterion becomes
After the matrix degradation of a gayer, the transverse and shear stiffness decrease to a small vafue. Only the elastic modulus in the fiber direction is preserved and all the coupling terms in the stress-strain relations are relaxed. However, in filament wound composite materials the degraded elastic and shear moduli are not negligible. For this reason, after the matrix degradation, the compliance matrix has been changed into
where Ed and Gd are the transverse and shear moduli after degradation obtained by experimental tests. Three possibilities have been retained in this finite element. (a) The transverse and shear stresses in the degraded Iayer of the composite are assumed to be maintained at their value. That means
where t+* and E* are the stress and strain vectors at degradation. (b) The transverse and shear stresses in a degraded layer are relaxed: a =
I&l&
(21)
(c) The transverse and shear stresses in the degraded layer decrease Iinearly from the failure value to zero in such a way that the failure principal strain energies (i.e. Kij) keep the Same value (cf. Fig. 5 and [4]). When the l~n~tu~nal stress reaches a critical value (different for each failure criterion) the fibers break, aII the stresses are then relaxed and are carried by the neighbuuring layers. These total failure criteria are:
142
S. Idelsohn
et al., Degradation
analysis of composite
materials
MATRIX _LEGRAMTION ------._______________
Fig. 5. Third model of post-degradation
(a) Tsai-Hill
: 2>H 1 *
(au/X)
(b) Tsai-Wu Flail
behavior.
(22)
: + i+T:,
(c) Maximum
2 1.
(23)
stress or maximum
crI12XI
or
strain:
El1 2 Xf .
(24)
(d) Sandhu: (&I( 2 Ed (uniaxial)
6. Newton-Raphson
.
(25)
method in composite materials
The Newton-Raphson method is used to solve the nonlinear problem. Two special considerations must be taken when implementing this algorithm: the correct integration of the constitutive equations in the presence of multimodulus material or during the matrix degradation and the possibility of numerical degradations during geometrically nonlinear problems. 6.1. Geometrically The structural
linear problems discretization
of the virtual
statement
yields the nonlinear
problem (26)
r(q) = f(4) - g = 0 where g are the generalized external loads and f(q) the internal adding the internal force contribution of each element B’u dv I Li
forces that can be obtained
by
(27)
S. Idelsohn et ul., Degradation analysis of composite materials
in which
B
is the matrix connecting
the displacement
parameter
143
with the strain tensor
A.e=BAq.
e=Bq, The Newton-Raphson
(28)
method consists in adopting an estimate q to the solution
q=B+Aq,
(2%
a=$++~,
(30)
and assuming that in the neighbourhood
of q the linear solution
Ao, = C~~,~)AE is a good approximation. r(q)
A presumably
(31) From (27)-(32), (27) becomes
(B’CB)
=
dv Aq f
better approximation
KAq=g-1
v
I0
B’C(6,2)B
Thus, Newton’s method iteratively by (33) qk+l=
qk
-
B’S dv - g.
(32)
can be obtained by putting (32) to zero, which yields
B’&dv
where the tangent stiffness matrix K =
I”
(33) K
is defined as
dv .
takes an initial approximation
(34) go
to q and attempts
to improve
it
(35)
Kh’k.
Then the new strain and stress inurement values are calculated by
A-%= B(qic+l-
qk)
(36)
and (37) and a new tangent matrix is evaluated. The procedure is stopped when the residual vector (r) becomes less than a desired precision. It is in the integration of the constitutive relations (37) that the composite material has a particular behavior.
144 6.1.1.
S. Idelsohn
et al., Degradation
analysis
of composite
materials
In N-linear problems
If the matrix degradation is reached in any point of the structure, or if a change in the stress sign occurs for a multimodulus composite material, the stress increment is written
where N is the number of stress-strain relation changes during the iteration. To evaluate Aci, a factor ri is calculated at each step so that ri As = hei is the part of the increment necessary to reach a discontinuity point (degradation or stress sign change). The stresses are evaluated at this point and a new strain-stress matrix is calculated. The process is repeated until the total strain increment is obtained. 6.1.2. Nonlinear material problems
In this case the total strain increment is subdivided into m sufficiently small parts and the total stress increment is evaluated by Euler’s explicit integration. In this way it is easy to verify when the matrix degradation or a stress sign change occurs. 6.2. Geometrically nonlinear problems The use of the Newton-Raphson’s method to solve geometrically nonlinear problems involving composite materials may introduce some numerical difficulties. In fact, when geometrical nonlinearities occur, the displacement fields are over-evaluated at the first iterations. Thus, it is possible that a failure criterion is verified in a layer before the real stress reaches the collapse value. Owing to the irreversiblity of the phenomena the algorithm may converge to a non physical solution. In order to improve this situation the nonlinear geometrical iterations are carried out separately from the nonlinear material iterations. In this way the problem is first considered as a linear material problem and the geometrically nonlinear iterations are performed. When the equilibrium is attained the correct numerical integration of the constitutive relation is achieved and a new step with only geometrically nonlinear iterations is performed.
7. Numerical results 7.1. Axisymmetric cylinder
To test the behavior of the multilayer finite element developed here and to compare the different failure criteria used, a simple four layer cylinder submitted to a shear load at the end was examined [9]. The structure is modeled by two elements of two layers in the radial direction and only one element in the axial direction (Fig. 6). The fiber direction in each layer is of 90”, O”,45” and 75” with respect to the axial direction starting from the internal layer to the external one. The
S. Idelsohn et al., Degradation analysis
. 0 C. 0
of composite materials
145
TSAI-HILL TSAI-vvu SANDHU sigma max. “c
Fig. 6. Axisymmetric cylinder; maximal radial displacement for the different faifure criteria.
mechanical properties
of the glass epoxy composite material before degradation
E,=5500hb, zf12=
v23
=
E2=Eg= VI3
=
1800hb,
G,, = G13= GZ3= 900 hb ,
0.25 )
and after the matrix failure of a layer, the transverse reduced to Ed= lhb,
are
and shear moduli are considerably
G,=0.3hb.
The longitudinal strength is 105 hb; the transverse and shear strengths are only 7 hb and 5 hb respectively, and for the Tsai-Wu criterion, the shear strength at 45” is 7.5 hb. The axisymmetric shear load applied at the end of the cylinder was incremented by a first step of 1.5 hb and then by 7 steps of 0.5 hb each. The maximal radial displacement (point c) for the different failure criteria used is shown in Fig. 6. It can be noted that with a shear load of 1.5 hb no degradation is reached for any criteria and in the next increment in ail cases one layer is degraded. However, between a shear load of 2 hb and 4 hb the failure of the 2nd and 3rd layer is reached at different moments with each failure criterion used and then different solutions are obtained. When the maximal shear load (5 hb) is applied, neither layer of the cylinder is completely broken. 7.2. Spherical ring with different behavior in tension and compression The second example considered is that of a clamped spherical ring subjected to an initial compression load and then successive incremental rotations of the ring are carried out to introduce centrifugal load. Its geometric characteristics are summarized in Fig. 7. This spherical ring is built of 5 Iayers from a composite material made up of Kevlar 49 fibers and epoxy matrix. Their mechanical properties are described in Table 1. The spherical shell thickness is h = 2 mm and there is an initial layer of 0.4 h with meridia1 fiber direction. Above and below this layer there are respectively 2 layers of 0.15 h each with fiber directions oriented at +45” and -45” with respect to the internal layer.
146
S. Idelsohn
et al., Degradation
analysis
of composite materials
Fig. 18. Spherical ring discretization.
The structure is modeled by using 7 axisymmetric finite elements of 5 layers each. The density of the external element is considered 10 times larger than the composite density. This is to take into account some concentrated mass placed there. The maximum stress criterion is used assuming that after degradation transverse and shear stress are maintained at their degradation values. The stress-strain relations are considered bilinear after degradation and the geometrical nonlinear strain is neglected. To check the importance of the multimodulus behavior in structures where the sign of the stress state varied, three different tests are performed (1) With only tensile modulus. (2) With only compressive modulus. (3) With the multimodulus (model developed in Section 3.3). Fig. 8 shows the maximal displacement at the boundary of the ring for the three tests. In the first increment only the initial compression load exists. Then the rotations are incremented to 31.6, 400.0, 700, 1000, 1200 and 1400 rad/sec and produce then a tension state in the structure. Table 1 Material properties of the composite
layer
Longitudinal tensile modulus E: Transverse tensile modulus E; Longitudinal compressive modulus E’, Transverse compressive modulus ES Shear modulus G12 Tension Poisson’s ratio Compression Poisson’s ratio Density
8.500 560 6.000 700 210
kg/mm’ kg/mm’ kg/mm’ kg/mm’ kg/mm’
0.34 0.3 1.56 gm/cc
Longitudinal tensile strength Transverse tensile strength Longitudinal compressive strength Transverse compressive strength Shear strength
141 2.8 2.8 14.1 4.5
kg/mm’ kg/mm’ kg/mm2 kg/mm’ kg/mm’
S. Idelsohn et al., Degradation analysis of composite materials
___o_._
CCMPRESSIVE
-o-
TENSILE
+
MULTlMOWUJS
147
MODUWS
MODULUS
Fig. 8. Spherical
ring; maximal
displacement.
At a rotational velocity of 1000 rad/sec, 2 layers reach the degradation in the multimodulus material instead of 1 layer with the traction test and 5 layers in the compression one. At the final rotational velocity of 1400 rad/sec 12 layers reach the degradation instead of 18 layers using the compression modulus. The maximal displacement obtained in this case is 34% above the value corresponding to the multimodulus theory. In fact, a small difference between the tensile and compressive moduli introduces for a given load a difference in the number of degraded layers. When new increments are performed this difference produces a new break-up and so on. This behavior explains the large differences encountered between a multimodulus analysis and a single modulus one in the presence of matrix degradation. 7.3. Geometrically nonlinear cantilever beam To test the element behavior in geometrically nonlinear problems a very simple cantilever beam is discretized by a plane stress element made up of four layers each oriented at O”, 90”, 45” and -45”. Table 2 Geometrically
nonlinear
cantilever
beam
Number of geometrically nonlinear iterations Increment 1 2 3
Tsai-Hill 6 5 4+4
Sandhu 6+1 5+1 5+1+2+1
Number of material nonlinear iterations Tsai-Hill 1 1 2
Sandhu 2 2 4
Matrix degradation non non yes
148
S. Idelsohn
.75-
et al., Degradation
analysis
SIRESS-SlRAIN RELATIONSFORWWU
Fig. 9. Cantilever
of composite materials
CRIlERK)N
beam.
All the mechanical properties before and after degradation are the same as those used in the axisymmetric cylinder described in Section 7.1. Two different failure criteria are tested here: Tsai-Hill and Sandhu. In the Sandhu criterion, nonlinear stress-strain relations are used (Fig. 9). Nevertheless these nonlinearities are mild and the two failure criteria can be compared. The external load is a uniform pressure p distributed on the top of the beam. Three increments are applied to obtain respectively: 2.5, 5.0 and 7.5 hb. Table 2 summarizes the performances obtained using the two failure criteria. It can be seen that in both cases, 2 layers reach the matrix degradation at the 3rd increment. Six geometrically nonlinear iterations are necessary for the first increment and 5 for the second. The Sandhu criterion needs a second nonlinear material iteration for these 2 increments. The third increment shows clearly the uncoupled behavior between the geometrical and material iterations: when the Tsai-Hill failure criterion is used, 4 geometrically nonlinear iterations are required to obtain equilibrium, In the subsequent mechanical iteration two layers degrade and to achieve a new equilib~um, 4 geometrical iterations are again needed. A similar behavior is observed by using Sandhu’s failure criterion.
8. Conclusions Some difficulties inherent to the development of a finite element for the analysis of fiber reinforced composite materials have been presented. The proposed element has been conceived for use in a research laboratory in order to compare numerical solutions with experimental results. For these reasons a large number of possibilities and criteria are present in the model. The desire of the authors is to present a mathematical tool to aid a better unde~tanding of the complicated composite material behavior. Unfortunately, very few experimental results are available in the literature in order to confirm the analysis. Nevertheless, the examples show that the solution of a problem is
S. Idelsohn et al., Degradation analysis of composite materials
149
sensitive to the failure criteria used and that the differences between the tension and compression moduli must not be neglected when change in the stress signs occurs. Future research and numerical tests should be devised in order to adjust the coupling terms cy, /3 in the strain-stress relation of material with different moduli in tension and compression, the degraded elastic and shear moduli E, and Gd and the differences between tension and compression stiffness values. The effect of shear stresses in the multimodulus material should be explored to ascertain their influence in the coupling terms. The analysis of decreasing load and numerically decreasing stresses after the matrix degradation are actually unsolved problems.
References [l] G. Sander and C. Nyssen, Modelisation des materiaux composites dans les methodes d’elements finis, L.T.A.S. SF-90, presented at XIVth Internat. Aeronautical Congress, Paris, 1979. [2] R.M. Jones, Mechanics of Composite Materials (McGraw-Hill, New York, 1975). [3] E. Hinton, The flexural analysis of laminated composite using a parabolic isoparametric plate bending element, Internat. J. Numer. Meths. Engrg. 11 (1979) 174-179. [4] C. Nyssen, Modelisation par elements finis du comportement non lineaire de structures aerospatiales, Doctoral Thesis, Aerospace Laboratory, University of Liege, 1979. [5] C. Nyssen and P. Beckers, Finite element linear and nonlinear analysis of composite materials, Proc. Conf. Nonlinear Finite Element Analysis and ADINA, MIT, Cambridge, MA, 1979. [6] R.S. Sandhu, Ultimate strength analysis of symmetric laminates, Tech. Rept. AFFDL-TR-73-137, 1974. [7] R.M. Jones, Stress-strain relations for materials with different moduli in tension and compression, AIAA J. 15(l) (1977) 16-23. [8] R.M. Jones and D.A.R. Nelson, Material models for nonlinear deformation of graphite, AIAA J. 14(6) (1976) 709-717. [9] G. Laschet, S. Idelsohn and M. Hogge, Element multicouche pour l’analyse des r&sines armees, L.T.A.S. SF-91, Aerospace Laboratory, University of Liege, 1979. [lo] S.W. Tsai, Strength characteristics of composite materials, NASA-CR-224, NASA, 1965. [ll] S.W. Tsai and E. Wu, A general theory of strength of anisotropic materials, J. Comput. Materials 5 (1971) 58-80. [12] J.L. Walsch, J.H. Ahlberg and E.N. Wilson, Best approximation properties of the spline fit, J. Math. Mech. 2(2) (1962) 225-233.