Potential of a simple variational analysis in predicting shear modulus of laminates with cracks in 90-layers

Abstract

The potential of accurate modelling of the shear modulus reduction of laminates with cracked 90-layers using models based on the minimization of the complementary energy with improved stress description in the constraint layers is evaluated.

This group of models refine Hashin’s model by introducing shape functions with unknown parameters to represent the out-of-plane shear stress distribution across the constraint layer thickness. The Hashin’s model becomes a particular case of the presented when the shape parameter approaches to zero. The most accurate shape parameters are found in the result of minimization. Three models are compared: the present variational model, Hashin’s model and the shear lag model introduced by Soutis which assumes linear out-of-plane shear stress distribution over an unknown part of the layer. It is shown in this paper that the size of this part may be determined by minimization of the complementary energy. The present model is the most accurate amongst the three, whereas Soutis’ model is more accurate than the Hashin’s model for laminates with constraint layer, thicker than the cracked layer. The comparison with finite element method results shows reasonable agreement. Agreement can be improved developing models with better description of the stress state in the cracked layer.

Keywords

Shear modulus intralaminar cracking complementary energy

Introduction

In laminated composite materials, many different types of micro-damage entities may evolve without causing final failure of the structural element. However, they are responsible for stiffness reduction and initiating of other, more critical, damage modes. One of the first and most common damage mode in laminates, causing degradation of thermo-elastic properties, is intralaminar cracking in layers with off-axis orientation with respect to the main loading direction. Large amount of experimental data has been gathered regarding the reduction of the elastic modulus and Poisson’s ratio of damaged laminates. Numerous models based on analytical approaches and/or numerical routines have been developed to study the relationship between the damaged state and the modulus of the laminate.

The situation is entirely different regarding the in-plane shear modulus of the damaged laminate. Only a few experimental investigations are described in the literature and often the data are not reliable. For example, Tsai et al.¹ developed a methodology for shear modulus determination using an experimental set-up where a laminate pre-damaged in tension is subjected to in-plane tangential displacement in the middle part. The shear modulus was calculated using a Timoshenko beam approximation. The problem, in general, is that it is difficult to introduce intralaminar cracks in a controlled way during the shear test. Cracks have to be introduced, for example, during a tensile test with transverse tensile stress in layers. Unfortunately, the requirements on specimen geometry in tension and in shear are not compatible.

Another problem of using experimental data for model validation is related to the non-linear stress–strain response to shear of unidirectional (UD) composites. For a layer in the laminate, this response is non-linear not only because of the evolving damage state, but also because of non-linear viscoelastic and viscoplastic effects.² For the above reasons, analytical model predictions are usually validated compared with the results of numerical solutions.

Hashin³ used the principle of minimum complementary energy. The admissible stress field, which satisfies equilibrium as well as the boundary and interface conditions in tractions, was constructed assuming that the in-plane shear stress distribution across each ply thickness is uniform and the single unknown function describing the in-plane shear stress distribution depends on the distance from the crack only. Using the stress distributions obtained by minimization, the lower bounds for the shear modulus of the laminate were obtained.

Two entirely different approaches were used^1,4–7: (a) Tan et al.^4,5 obtained expressions for axial modulus, Poisson’s ratio and shear modulus of the cross-ply laminate with 90-cracks, by integrating the equilibrium and constitutive expressions over the ply thickness and obtaining a second order differential equation for stress distributions; (b) Tsai et al.^1,6 and Abdelrahman et al.⁷ reduced the 3-D elasticity problem to 2-D problem in terms of displacements¹ or stresses⁷ averaged over the ply thickness. In both cases, the obtained set of differential equations with constant coefficients was solved analytically. Due to the assumed linear through the thickness dependence of the out-of-plane shear stresses, the results of these models coincide with Hashin’s result.

Henaff-Gardin et al.⁸ analysed a double-cracked cross-ply laminate in a similar manner as in Ref. 1 just using a simpler shear lag model with parabolic opening displacement and uniform sliding displacement distributions. For the same problem, Kashtalyan et al.^9,10 adjusted the effective properties of the constraint layer for damage when analysing the local stresses in another layer. This leads to an iterative procedure when cracks are present in both 0- and 90-layers of the cross-ply laminate. It was shown that the interaction of cracks in two layers leads to considerable additional reduction of the laminate’s shear modulus.

Fan et al.¹¹ presented constitutive equations for a layer with cracks containing the so-called ‘in-situ damage effective functions (IDEF)’ which depend on the crack density of the lamina and on the neighbouring layer constraints. In order to determine IDEF, they introduced ‘an equivalent constraint model’, which assumes that the constraint of the lay-ups above and below the analysed lamina is the same as from two homogenized orthotropic sublaminates and the actual laminate was replaced by a cross-ply. The constitutive relationships for the damaged layers were used in the framework of the CLT to obtain the stiffness matrix of the damaged laminate.

This approach was further refined by the Soutis research group¹² where the local stress problem was solved using an improved shear lag model which assumed that the intralaminar shear stress in the 0-layer differs from zero only in a zone close to the interface. This idea which is supported by finite element calculations was not consistently implemented: in Ref. 12, the authors claimed that the perturbation zone covers one prepreg layer, whereas in more recent publications^9,10,13 they are extending the perturbation zone over the whole homogenized sublaminate. Obviously, the problem is that the size of this zone becomes a fitting parameter.

Their hypothesis that the out-of-plane shear stress in the constraint layer can be described by a decreasing linear function of the distance from the interface (the maximum is at the interface and the value is zero at distance z₁) will be inspected in this paper. To find z₁, a simple calculation routine based on the principle of minimum complementary energy will be described and used.

The objective of the current work is to address the potential of improving the conservative predictions for [S/90_n]_s and [90_n/S]_s laminates using variational models with more refined stress description in the constraint layer. A more general variational model is introduced to predict shear modulus reduction of laminates with intralaminar cracks. This model includes Hashin’s model as particular case and can be used to find the shear zone size in the model based on Soutis assumptions. In contrast to Hashin’s model, which assumes linear distribution, the novel model allows an arbitrary shape of the out-of-plane shear stress distribution in the constraint layer. Selecting a specific form of the shape function with unknown parameters, the principle of the minimum complementary energy can be used to find these parameters. The bilinear distribution¹² is one particular case considered and the shape parameter to find is z₁. The most accurate model, between the ones that were compared, can be identified by comparing the complementary energy values at the minimum; the lowest value corresponds to the most accurate stress distribution. It is shown that the improved stress description in the constraint layer is not sufficient for accurate predictions of the laminate shear modulus. The approach for further improvement of the model is outlined.

Repeating element and assumptions

Model formulation

We consider a repeating element of the [S,90_n]_s or [90_n S]_s laminate with a periodic system of cracks in the 90-layers as shown in Figure 1. S is the notation used for the homogenized adjacent lay-up. Its elastic properties are as for a balanced ‘sub-laminate’ with the same lay-up. For simplicity and in order to maintain symmetry, with respect to the mid-plane, we assume that the lay-up of the sub-laminate on the top and on the bottom of the 90-layer is the same. In the simplest case, the sublaminate is a UD 0-layer. The crack spacing is noted as $2 l_{0}$ , the 90-layer thickness and the sub-laminate’s thickness is shown in Figure 1. The laminate’s overall thickness is noted as $2 h = 2 (b + d)$ . The in-plane shear modulus of the undamaged laminate, $G_{xy 0}^{LAM}$ can be calculated using classical laminate theory (CLT).

G_{xy 0}^{LAM} = \frac{1}{h} (G_{xy}^{S} b + G_{12} d)

(1)

In equation (1),

G_{xy}^{S}

is the in-plane shear modulus of the sub-laminate which may be calculated using CLT and

G_{12}

is the shear modulus of the UD composite. The laminate is subjected to average in-plane shear stress

τ_{0}

. Then the undamaged laminate shear strain,

γ_{xy 0}^{LAM}

and stresses in layers,

τ_{0}^{90}, τ_{0}^{S}

are:

γ_{xy 0}^{LAM} = τ_{0} / G_{xy 0}^{LAM}, τ_{0}^{90} = G_{12} γ_{xy 0}^{LAM}, τ_{0}^{S} = G_{xy}^{S} γ_{xy 0}^{LAM}

(2)

Figure 1.

Intralaminar cracks in symmetric and balanced laminates. (a) [S/90n]s laminate with cracks in central 90-layer and (b) [90n/S]_s laminate with cracks in surface 90-layer.

Assumed stress profiles

Independent on the location of the cracked 90-layer (on the surface or in the center of the laminate), we use the same assumptions and the same method of solution. It was shown in Refs. 3 and 14 that for the case when the stress distributions does not depend on the y-coordinate (infinite material in y-direction) the two stress components $σ_{xy}$ and $σ_{yz}$ are separated in equilibrium equations from the rest of stress components. If in addition the layers in the laminate are orthotropic in the used system of coordinates, the stress distribution problem in macroscopic in-plane shear loading does not include other stress components. If, in contrary, some of the layers are monoclinic (for example off-axis layers) the complementary energy expression will include all stress components (due to coupling terms in stress-strain relationships). Hence, the homogenization over $\pm θ$ layers suggested in this paper to obtain homogenized sub-laminate is not entirely correct. The same interaction appears if the sub-laminate is not balanced.

We assume that the in-plane shear stress in the cracked 90-layer, $σ_{xy}^{90}$ does not depend on the thickness coordinate z, whereas in the sub-laminate supporting this layer the z-distribution is given by an unknown function $f' (z)$ . Here and in following, $()'$ denotes derivative of the function of one variable. In the sub-laminate, the shear stress is assumed in a form of product of two functions; one depending on x and the second on z. After introducing two unknown stress perturbation functions $Φ (x)$ and $Φ_{1} (x)$ , to describe the stress distribution in x-direction in 90-layer and in the sub-laminate, respectively, we use the following form of the in-plane shear stress distribution:

σ_{xy}^{90} = τ_{0}^{90} \cdot (1 - Φ (x))

(3)

σ_{xy}^{S} = τ_{0}^{S} \cdot (1 - d \cdot Φ_{1} (x) \cdot f' (z))

(4)

These forms are a generalization of the expressions used by Hashin,³ in which these stress components were assumed to be z-independent.

The shear stress in the 90-layer has to satisfy the traction-free condition on crack surfaces:

σ_{xy}^{90} (x = \pm l_{0}) = 0

(5)

The solution for

σ_{yz}

has to satisfy the stress-free boundary conditions on

z = 0; h

. In addition to that, the equilibrium equations in each point and in average, as well as the interface conditions in tractions have to be satisfied. The principle of minimum complementary energy will be used in the following section to insure displacement continuity in average. The derivation of all stress expressions is now demonstrated in detail for the case of central transverse 90-layer only.

Expressions for stress components

Cracked central 90-layer

From the shear force balance

\int_{0}^{h} σ_{xy} d z = τ_{0} (d + b)

(6)

we obtain

Φ_{1} (x) = - Φ (x) \frac{G_{12}}{G_{xy}^{S}} \frac{1}{f (h) - f (d)}

(7)

leading to

σ_{_{xy}}^{S} = τ_{0}^{S} + τ_{0}^{90} \cdot d \cdot \frac{f' (z)}{f (h) - f (d)} \cdot Φ (x)

(8)

Substituting

σ_{xy}^{S}

in the equilibrium equation

\frac{\partial σ_{yx}}{\partial x} + \frac{\partial σ_{yz}}{\partial z} = 0

(9)

and integrating it with respect to x we obtain after applying stress free condition on the surface

z = h

σ_{yz}^{S} = τ_{0}^{90} \cdot d \cdot [\frac{f (h) - f (z)}{f (h) - f (d)}] \cdot Φ' (x)

(10)

The same equilibrium equation + boundary condition at

z = 0

in 90-layer yield

σ_{yz}^{90} = τ_{0}^{90} \cdot Φ' (x) \cdot z

(11)

It can be checked that continuity of

σ_{yz}

at the layer interface is satisfied. We introduce a new unknown function

ϕ (z)

as:

ϕ = \frac{f (h) - f (z)}{f (h) - f (d)}

(12)

Then, the following expressions can be written for stress components in the case of cracked central layer (lower index c is used for ‘central’)

σ_{xy}^{90} = τ_{0}^{90} \cdot (1 - Φ_{C} (x))

(13)

σ_{xy}^{S} = τ_{0}^{S} - τ_{0}^{90} \cdot d \cdot Φ_{C} (x) \cdot ϕ'_{C} (z)

(14)

σ_{yz}^{90} = τ_{0}^{90} \cdot Φ'_{C} (x) \cdot z

(15)

σ_{yz}^{S} = τ_{0}^{90} \cdot d \cdot ϕ_{C} (z) \cdot Φ'_{C} (x)

(16)

Finally, the unknown functions

ϕ_{C} (z)

and

Φ_{C} (x)

have to satisfy the boundary conditions:

ϕ_{C} (h) = 0, ϕ_{C} (d) = 1, Φ_{C} (\pm l_{0}) = 1

(17)

Cracked surface 90-layer

The derivation of expressions is very similar to the case of a central cracked layer and all stress components are expressed through unknown functions $Φ_{S} (x)$ and $ϕ_{S} (z)$

σ_{xy}^{90} = τ_{0}^{90} (1 - Φ_{S} (x))

(18)

σ_{xy}^{S} = τ_{0}^{S} - τ_{0}^{90} \cdot d \cdot Φ_{S} (x) \cdot ϕ'_{S} (z)

(19)

σ_{yz}^{90} = - τ_{0}^{90} \cdot Φ'_{S} (x) \cdot (h - z)

(20)

σ_{yz}^{S} = τ_{0}^{90} \cdot d \cdot Φ'_{S} (x) \cdot ϕ_{S} (z)

(21)

In this case

ϕ_{S} (b) = - 1, ϕ_{S} (0) = 0, Φ_{S} (\pm l_{0}) = 1

(22)

Minimization of the complementary energy

Thermal stresses do not affect the calculated elastic properties and therefore thermal terms can be neglected in the expression for complementary energy. The complementary energy expression becomes equal to the strain energy. In spite of this equality, there is a fundamental difference in using minimum principle for complementary energy and minimum principle for strain energy. In the former case, the used admissible stress field for minimization has to satisfy automatically the equilibrium conditions and all boundary and interface conditions in tractions, whereas in the latter case conditions of displacement continuity are automatically satisfied by the assumed functions. The strain energy of a repeating element between two cracks in 90-layer is

U S = \frac{1}{2} \cdot \int_{- l_{0}}^{l_{0}} \int_{h_{S}} [\frac{(σ_{xy}^{S}) 2}{G_{xy}^{S}} + \frac{(σ_{yz}^{S}) 2}{G_{yz}^{S}}] d z d x

(23)

U 90 = \frac{1}{2} \cdot \int_{- l_{0}}^{l_{0}} \int_{h_{90}} [\frac{(σ_{xy}^{90}) 2}{G_{12}} + \frac{(σ_{yz}^{90}) 2}{G_{12}}] d z d x

(24)

Symbols

\int_{h_{S}}

and

\int_{h_{90}}

are used for integration over the corresponding layer thickness. After substituting equations (13)–(16) (or equations 18–21) in equations (23) and (24) and performing z-integration, we obtain the following expression for complementary energy of the laminate in case of a cracked 90-layer

\begin{matrix} U LAM = U S + U 90 = U_{0}^{LAM} + \frac{(τ_{0}^{90}) 2 d}{2 G_{12}} \\ \times \int_{- l_{0}}^{+ l_{0}} [A_{1} \cdot Φ 2 (x) + A_{2} \cdot d 2 \cdot [Φ' (x)] 2] d x \end{matrix}

(25)

Since equation (25) is valid also for laminate with cracks in surface 90-layer, indexes C or S in Φ are omitted. For cracked central 90-layers:

A_{1} = 1 + \frac{G_{12}}{G_{xy}^{S}} K_{1 C} A_{2} = \frac{1}{3} + \frac{G_{12}}{G_{yz}^{S}} K_{2 C}

(26)

K_{1 C} = d \int_{d}^{h} [ϕ'_{C} (z)] 2 d z K_{2 C} = \frac{1}{d} \int_{d}^{h} [ϕ_{C} (z)] 2 d z

(27)

For cracked surface 90-layers:

A_{1} = 1 + \frac{G_{12}}{G_{xy}^{S}} K_{1 S} A_{2} = \frac{1}{3} + \frac{G_{12}}{G_{yz}^{S}} K_{2 S}

(28)

K_{1 S} = d \int_{0}^{b} [ϕ'_{S} (z)] 2 d z K_{2 S} = \frac{1}{d} \int_{0}^{b} [ϕ_{S} (z)] 2 d z

(29)

In the following, indices C and S are omitted only if the expression has the same form for both locations of the 90-layer.

The function Φ that minimizes equation (25) has to satisfy the Euler–Lagrange equation (30):

\frac{d}{d x} (\frac{\partial L}{\partial Φ'}) - \frac{\partial L}{\partial Φ} = 0

(30)

which leads to the following equation for

Φ (x)

determination:

d 2 \cdot A_{2} \cdot Φ ″ (x) - A_{1} \cdot Φ = 0

(31)

The solution of equation (31) which satisfies the boundary condition equation (17) or (22) is

Φ = \frac{\cosh (μ \frac{x}{d})}{\cosh (μ \frac{l_{0}}{d})}

(32)

where

μ = \sqrt{\frac{A_{1}}{A_{2}}} = \sqrt{\frac{1 + \frac{G_{12}}{G_{xy}^{S}} K_{1}}{\frac{1}{3} + \frac{G_{12}}{G_{yz}^{S}} K_{2}}}

(33)

The equation (32) gives the minimum of the complementary energy for the used selection of

ϕ (z)

which means that it can be further minimized by selection of

ϕ (z)

. For this purpose, the complementary energy equation (25) can be rearranged as follows:

\begin{matrix} U LAM = U_{0}^{LAM} + \frac{(τ_{0}^{90}) 2 d}{2 G_{12}} {\int_{- l_{0}}^{+ l_{0}} A_{1} \cdot Φ 2 (x) d x \\ + A_{2} \cdot d 2 [Φ (x) \cdot Φ' (x) |_{- l_{0}}^{+ l_{0}} - \int_{- l_{0}}^{+ l_{0}} Φ (x) Φ ″ (x) d x] \end{matrix}}

(34)

By replacing

Φ ″

in the second integral via Φ, as follows from equation (31), we obtain

U LAM = U_{0}^{LAM} + \frac{(τ_{0}^{90}) 2 d 3}{2 G_{12}} A_{2} 2 Φ' (l_{0})

(35)

Substituting in equation (35) A₂ given by equations (26) or (28) and expressing

Φ' (l_{0})

from equation (32) we obtain

U LAM = U_{0}^{LAM} + \frac{(τ_{0}^{90}) 2 d 2 μ}{G_{12}} (\frac{1}{3} + \frac{G_{12}}{G_{yz}^{S}} K_{2}) \tanh (μ \frac{l_{0}}{d})

(36)

By selecting different functions

ϕ (z)

, we will have different values of K₂ and μ, leading to different values of

U_{LAM}

. The lowest value will correspond to the most accurate stress distribution in average.

Shear modulus expression

The shear modulus of the damaged laminate can be determined using $u LAM$ – the strain energy density of the repeating element:

γ_{xy}^{LAM} = \frac{\partial u LAM}{\partial τ_{0}}

(37)

The strain energy density can be obtained by dividing equation (36) with the volume of the repeating element, while the shear stress

τ_{0}^{(90)}

in equation (36) can be expressed through

τ_{0}

using equation (2). The undamaged laminate strain energy density is

τ_{0}^{2} / 2 G_{xy 0}^{LAM}

. After performing differentiation according to equation (37), we obtain

\begin{matrix} γ_{xy}^{LAM} = \frac{τ_{0}}{G_{xy 0}^{LAM}} + \frac{τ_{0} G_{12}}{(G_{xy 0}^{LAM}) 2} \\ \times \frac{d 2}{(d + b) l_{0}} μ (\frac{1}{3} + \frac{G_{12}}{G_{yz}^{S}} K_{2}) \tanh (μ \frac{l_{0}}{d}) \end{matrix}

(38)

Since the Hook’s law for the damaged laminate has to have the form

γ_{xy}^{LAM} = \frac{τ_{0}}{G_{xy}^{LAM}}

(39)

the following shear modulus expression for the damaged laminate is obtained:

\frac{G_{xy}^{LAM}}{G_{xy 0}^{LAM}} = \frac{1}{1 + \frac{G_{12}}{G_{xy 0}^{LAM}} \frac{d}{d + b} \frac{d}{l_{0}} μ (\frac{1}{3} + \frac{G_{12}}{G_{yz}^{S}} K_{2}) \tanh (μ \frac{l_{0}}{d})}

(40)

For cross-ply laminates,

G_{xy}^{S} = G_{12}

G_{yz}^{S} = G_{23}

and

G_{xy 0}^{LAM} = G_{12}

. Then the expression for the shear modulus can be re-written as

\frac{G_{xy}^{LAM}}{G_{xy 0}^{LAM}} = \frac{1}{1 + \frac{d}{d + b} \frac{d}{l_{0}} μ (\frac{1}{3} + \frac{G_{12}}{G_{23}} K_{2}) \tanh (μ \frac{l_{0}}{d})}, μ = \sqrt{\frac{1 + K_{1}}{\frac{1}{3} + \frac{G_{12}}{G_{23}} K_{2}}}

(41)

The shear modulus reduction is usually¹⁴ presented as a function of the crack density, ρ (cr/mm) or the normalized crack density,

ρ_{n}

ρ = \frac{1}{2 l_{0}}, ρ_{n} = \frac{t_{90}}{2 l_{0}}

(42)

where

t_{90}

is the thickness of the 90-layer. For surface layers

t_{90} = d

, whereas for central layers

t_{90} = 2 d

. This definition of the normalized crack density means that for the same relative distance between cracks

l_{0} / d

the

ρ_{n}

in the cracked central layer is two times larger than for the cracked surface layer.

Selection of the shape function $ϕ (z)$

Since it is expected that the shear stress perturbation in the sub-laminate will decay with increasing distance from the crack tip (layer interface), the absolute values of the $ϕ (z)$ function and its first derivative have to be decreasing monotonously with the increase of the distance from the interface.

Exponential shape function

A natural choice with these features is exponential function. The shape function will contain one shape parameter α, which will be found as the one which produces the lowest value of the complementary energy given by equation (36). The expected shape is shown schematically in Figure 2 as ‘LTU’.

Figure 2.

Profile of the shape function in the sub-laminate according to Hashin,¹ the present model (LTU) and Soutis¹²: (a) for cracked central 90-layer and (b) for cracked surface 90-layer.

Central 90-layer with cracks

The shape function that satisfies boundary conditions in equation (17) is selected in the form:

ϕ_{C} (z) = \frac{1}{e \frac{α \cdot b}{d} - 1} \cdot (e - \frac{α \cdot (z - h)}{d} - 1) ϕ'_{C} = - \frac{α}{d} \cdot \frac{e - \frac{α \cdot (z - h)}{d}}{e \frac{α b}{d} - 1}

(43)

This function increases exponentially approaching the transverse crack tip. Using equation (27), the parameters K₁ and K₂ are calculated

K_{1 C} = \frac{α}{2} \frac{e \frac{α b}{d} + 1}{e \frac{α b}{d} - 1} K_{2 C} = \frac{\frac{b}{d} - \frac{1}{2 α} (1 - e \frac{2 α b}{d}) + \frac{2}{α} (1 - e \frac{α b}{d})}{(e \frac{α b}{d} - 1) 2}

(44)

Surface 90-layer with cracks

The shape function that satisfies the boundary conditions in equation (22) is selected in the form

ϕ_{S} (z) = - \frac{e \frac{α \cdot z}{d} - 1}{e \frac{α \cdot b}{d} - 1} ϕ'_{S} = - \frac{α}{d} \frac{e \frac{α \cdot z}{d}}{e \frac{α \cdot b}{d} - 1}

(45)

By using equation (29), the parameters K₁ and K₂ are calculated. For this particular selection of

ϕ_{S}

K_{1 S} = K_{1 C} K_{2 S} = K_{2 C}

(46)

This means that the same value of n will give the most accurate solution for both ([S, 90_n]_s and [90_n S]_s) lay-ups.

Hashin's assumptions

It can be proven that the Hashin’s solution³ is a partial case of the selected functions (43) and (45). In Hashin’s model, the function $ϕ_{C} (z)$ representing the $σ_{yz}$ distribution is linear

ϕ_{Hashin} (z) = \frac{h - z}{b}

(47)

Writing the Taylor expansion of equation (43) for small α values we obtain

{lim}_{α \to 0} ϕ_{C} (z) = \frac{\neg 1 - \frac{α \cdot (z - h)}{d} - \neg 1}{\neg 1 + \frac{α \cdot b}{d} - \neg 1} = \frac{h - z}{b} = ϕ_{Hashin} (z)

(48)

In Hashin’s case:

K_{1 C, Hashin} = \frac{d}{b} K_{2 C, Hashin} = \frac{b}{3 d}

(49)

Bilinear shape functions by Soutis

In Ref. 12, the stress distribution in the 0-layers was assumed to be bilinear as schematically presented in Figure 2. The distance from the interface, where the stress perturbation becomes negligible, is denoted by z₁. The shape functions for cracking in central and in surface 90-layers are

\begin{matrix} ϕ_{c} = {\begin{matrix} 0, & d + z_{1} < z < h \\ \frac{z_{1} + d - z}{z_{1}}, & d \leq z \leq d + z_{1} \end{matrix} \\ ϕ_{s} = {\begin{matrix} 0, & 0 < z < z_{1} \\ \frac{b - z_{1} - z}{z_{1}}, & z_{1} \leq z \leq b \end{matrix} \end{matrix}

(50)

Using equation (50) in equations (27) and (29) we obtain:

K_{1 C, Soutis} = K_{1 S, Soutis} = \frac{d}{z_{1}} K_{2 C, Soutis} = K_{2 S, Soutis} = \frac{z_{1}}{3 d}

(51)

The value of the complementary energy depends on the unknown z₁.

Results and analysis

Minimum of the complementary energy

During the minimization of equation (36), the parameter α was increased with a small step and at each value the perturbation complementary energy (the second term in equation 36) was calculated. Since the minimum point does not depend on the applied stress and on the geometrical scale, this term was normalized dividing by $(τ_{0}^{90}) 2 d 2$ . An example of these calculations for [0/90]_s laminate with $ρ_{n} = 0.3$ is shown in Figure 3(a). There is only one minimum in the curve and it is usually in the region $α \in (0, 1.5)$ . It was found that the α value, which corresponds to the minimum, becomes smaller with increasing crack density as shown in Figure 3(b).

Figure 3.

Cross-ply laminate with $d / b = 1$ , G₁₂ = 4.58 GPa, G₂₃ = 3.4 GPa, $ρ_{n} = 0.3$ . (a) Perturbation energy dependence on the shape parameter α showing minimum in small α-value region and (b) α dependence on crack density.

For very high crack density, the shear modulus from equation (40) has to approach to predictions of the ply-discount model, which assumes that transverse and shear modulus of the damaged layer is zero. One can show that it is possible only if the shape parameter $α \to 0 .$ Thus, our solution in the very high crack density region coincides with the Hashin’s solution.

Comparison with Hashin’s model

The complementary energy, which was calculated by using the exponential shape function, was always lower than using Hashin’s model. As a consequence (a) the present model is more accurate for all geometrical and elastic parameters and (b) Hashin’s model predicts larger shear modulus reduction than the more accurate model.

However, as shown in Figure 4 for [S/90_n]_s laminates, the numerical differences in the predicted shear modulus are rather small. The well-known and expected trends shown in Figure 4 are as follows: higher shear modulus reduction when the cracked layer is relatively thick and the predictions are not sensitive to the layer shear modulus change of 20% or less (compare Figure 4(b) with Figure 4(c)); in the laminate, where the 0-layer is replaced by a sub-laminate with high in-plane shear modulus, the shear modulus reduction is much smaller. The comparison of the models for the cases presented in Figure 4 shows that (a) the improved model gives more noticeable difference for laminates with relatively thin cracked 90-layer; (b) the difference between models is slightly larger if the anisotropy in composite shear moduli increases, see Figure 4(b) and (c) (carbon fiber composites) and (c) for laminates with stiff sub-laminates, the difference is rather large in an intermediate crack density region.

Figure 4.

Reduction of the shear modulus of a cross-ply laminate due to cracking according to Hashin’s and the present model. The upper index LAM is omitted. (a) $d / b = 1, d = 0.3 mm$ , $G_{12} = 4.58 GPa$ , $G_{23} = 3.4 GPa$ ; (b) $d / b = 0.5$ , $G_{12} = 4.58 GPa$ , $G_{23} = 3.4 GPa$ ; (c) $d / b = 0.5$ , $G_{12} = 5.0 GPa$ , $G_{23} = 2.5 GPa$ and (d) $d / b = 0.5$ , $G_{12} = 5.0 GPa$ , $G_{23} = 2.5 GPa$ , $G_{xy}^{S} = 20 GPa$ , $G_{yz}^{S} = 2.5 GPa$ .

Parametric analysis of the [0_m/90_n]_s laminate using exponential shape functions

First, the thickness of the 0-layer was kept constant and equal to 0.3 mm and the thickness of the 90-layer was given different values. Specifically, values from 0.01–3 mm were studied with the step of 0.3 mm, i.e. the same as the assumed ply thickness. The values of the in-plane and out-of-plane shear modulus were kept constant as well; G₁₂ = 3.4 GPa and G₂₃ = 4.58 GPa. The results are presented in Figure 5. The arrow in Figure 5 indicates the direction of the 90-layer thickness increase; larger thickness leads to higher $G_{xy}^{LAM}$ reduction. It is shown that the thickness of the transverse layers affects both the saturation value of the $G_{xy}^{LAM}$ at high crack density and the slope of the curve.

Figure 5.

Effect of the 90-layer thickness on the shear modulus of the damaged [0/90n]s laminate.

In order to quantify the effect of the thickness of the 0°-layer on the laminate’s shear modulus reduction, the thickness of the transverse layer was kept constant (0.3 mm) and the thickness of the 0-layer was given the values 0.01–3 mm with an increment of 0.3 mm. The results are presented in Figure 6. As expected from the ply-discount model, the thickness of the longitudinal ply has a strong influence on the final shear properties of the laminate: for 3-mm thick 0-layer, the modulus reduction is only 10% at saturation, whereas for 0.3 mm (i.e. only one 0-ply) the reduction is close to 53%.

Figure 6.

Effect of the 0-layer thickness on the shear modulus of [0n/90]s cross-ply laminate.

Apart from the geometrical parameters and the in-plane shear modulus, the current model contains the out-of-plane shear modulus G₂₃ of the constraint layer as input. This property is usually not measured; it is estimated from the transverse isotropy if the constraint layer is a UD layer. The sensitivity of the results to the value of the out-of-plane shear modulus is shown in Figures 7 and 8. The geometrical parameters were kept constant (the thickness 0.3 mm for both 0- and the 90-layer) and G₁₂ = 4.58 GPa. G₂₃ was given values from 0.01–15 GPa. Figure 7 proves that for high crack density, the value of G₂₃ is not of great importance (the curve approach to the ply-discount value calculated using CLT). Figure 8 is the detail of Figure 7 at more realistic, low crack density. With higher G₂₃, the laminate’s shear modulus reduction versus crack density is much slower. This may be related to the decreased stress transfer length (size of the perturbation zone).

Figure 7.

Effect of $G_{23}$ on the shear modulus reduction of [0/90]s laminate.

Figure 8.

Effect of $G_{23}$ on the shear modulus reduction of [0/90]s laminate at low crack density.

There was no need to perform parametric analysis for similar laminates with cracks in the surface layer; as shown in the section Exponential shape function $K_{1 S} = K_{1 C}$ , $K_{2 S} = K_{2 C}$ and the value of α that gives the lowest value of complementary energy is the same. Hence, calculations using equations (40) or (41) for the same crack density (the same $l_{0} / d$ ) would give identical shear modulus reduction. Presenting the shear modulus reduction as a function of the normalized crack density, $ρ_{n}$ the predicted reduction is different due to reasons described in the section Shear modulus expression.

The conclusion that in both cases shown in Figure 1, the laminate’s shear modulus reduction with crack density is the same, was unexpected because the shape of the crack faces in the deformed state is different. The same result is obtained also for Hashin’s model (as a special case of our model) and also for Soutis model, see equation (51). Since this result could be due to the approximate nature of the discussed analytical solutions, results from finite element method (FEM), presented in¹⁵ were used to resolve this question. More exactly, in Ref. 15 fitting functions for the normalized average crack face sliding displacement were presented in the form:

u_{1 an} = A + B (\frac{G_{12}}{G_{xy}^{S}}) n

(52)

In Ref. 15 it was found from data fitting that

n = 0.82

in both cases, whereas

A = 0.3

B = 0.066 + 0.054 (d / b)

for central cracks and

A = 0.6

B = 0.134 + 0.105 (d / b)

for surface cracks. So, the

u_{1 an}

for surface cracks is almost exactly two times larger. Since in Ref. 15 the normalization of the

u_{1 an}

is with respect to the layer thickness (2d for central cracks and d for surface cracks) the real ‘non-normalized’ average sliding displacement is almost identical in both cases. Therefore, the effect of both cracks on the shear modulus reduction is also almost identical. This conclusion confirms our findings and should be incurraging for researchers developing approximate analytical models.

Comparison with FEM

The shear modulus reduction for GF/EP laminate with layer elastic properties and geometrical parameters given in Table 1 was compared with results in Ref. 15 obtained using FEM. In calculations ANSYS 8.1 with 21,000 3_D 8-node elements representing the upper part of the repeating element were used (more details regarding boundary conditions, etc. are given in Ref. 15). The accuracy of the FEM model was confirmed by complete agreement with FEM results presented in Ref. 16. The analysis involved three different stacking sequences; [0_n/90₂]_s, where n = 1, 2 and 4. The results obtained using the present model show correct trends but the values are lower than the FEM predictions, see Figure 9. That is expected since the current model is based on the principle of minimum complementary energy, which always gives a lower bound to exact results. The agreement with FEM is acceptable but not perfect. The main reason for differences is that the used exponential shape function with shape parameter obtained by minimization improves the stress state description mostly in the constraint layer (sub-laminate). In result, see Figure 4, the predictions are improved more comparing with Hashin’s model when the 0-layer is relatively thick. The agreement with FEM results improves for laminates with larger $b / d$ . For example, for the considered GF/EP with $b / d = 3$ the values at $ρ_{n} = 0.5$ are 0.88 and 0.91 according to the developed model and FEM. Reducing the in-plane shear modulus of the 90-layer to 3 GPa while keeping for 0-layer the value in Table 1, the analytical value at $ρ_{n} = 1$ is 0.66 whereas the FEM value is 0.68. Both these examples with a good agreement are for cases when the role of the 90-layer is lower. Apparently, the stress description in the 90-layer also has to be improved to come closer to the exact results in this case represented by the FEM solution.

Figure 9.

Comparison between FEM results and analytical modelling for [0_n/90₂]_s laminate.

Table 1.

Parameters used for the comparison in Figure 9.

Material	G12 (GPa)	G23 (GPa)	d (mm)	b (mm)
GF/EP	5.0	5.36	0.1	0.1

Accuracy of Soutis et al. model¹²

As described in the introduction, the research group by Soutis et al.¹² suggested to improve the shear lag assumption-based stress distribution models by replacing the common assumption of linear out-of-plane shear stress distribution across the constraint layer thickness by another assumption which was based on analysis of FEM results. It was assumed that $σ_{yz}^{S}$ is a linear function of z with the maximum at the layer interface. In contrast with the common assumption, stating that the stress is non-zero over the whole layer thickness, they suggested that the perturbation stress $σ_{yz}^{S}$ becomes equal to zero at some distance z₁ from the interface. The suggested distributions are shown in Figure 2. The value of the distance z₁ is unknown, but in Soutis et al.¹² it was assumed to be equal to one prepreg layer thickness. Much more reasonable, from the elasticity point of view, would be to assume that the perturbation zone z₁ is proportional to the size of the crack (the thickness of the 90-layer d).

Since the shape of the $σ_{yz}^{S}$ z-distribution is given by $ϕ_{S}$ and $ϕ_{C}$ , the implementation of the bilinear shape functions is very simple. The parameters for complementary energy calculation according to this model are given by equations (50) and (51). Since parameters K₁ and K₂ for surface cracks and central cracks are the same, there is no difference in the value of the complementary energy and in the shear modulus reduction for these two cases.

The best value of z₁ can be found as the one that produces the lowest value of the complementary energy. The complementary energy values can also be used to resolve the question about the most accurate among the three models (Hashin’s,³ Soutis¹² and the present). The comparison in this section is performed for central cracks. The complementary energy values according to these models are given in Tables 2 and 3 for several values of the layer thickness ratio. The z₁ values found by minimization are also presented there.

Table 2.

Normalized complementary energy of $[0_{n} / 90_{4}]_{s}$ at $ρ_{n} = 0.1$ , $G_{12} = 3.4 GPa$ and $G_{23} = 4.58 GPa$ .

n	Present	Hashin	Soutis	z₁/b
1	0.41228	0.4134	0.4134	1.00
2	0.34100	0.3444	0.3444	1.00
3	0.31718	0.3236	0.3236	1.00
4	0.30678	0.3170	0.3170	1.00
8	0.29723	0.3278	0.3161	0.57
12	0.29657	0.3522	0.3161	0.38
16	0.29652	0.3783	0.3161	0.29

Table 3.

Normalized complementary energy of $[0_{n} / 90_{4}]_{s}$ , $ρ_{n} = 0.8$ , $G_{12} = 3.4 GPa$ and $G_{23} = 4.58 GPa$ .

n	Present	Hashin	Soutis	z₁/b
1	0.41218	0.41333	0.41333	1.00
2	0.34007	0.34327	0.34327	1.00
3	0.31486	0.32043	0.32043	1.00
4	0.30293	0.31092	0.31092	1.00
8	0.28911	0.30591	0.30518	0.82
12	0.28713	0.31121	0.30518	0.55
16	0.28683	0.31704	0.30518	0.41

The complementary energy value according to the present model is the lowest for all configurations and crack densities which means that this model is the most accurate. Using Soutis model in cases when the 0-layer is thinner than the 90-layer $(b < d)$ results show that $z_{1} / b = 1,$ which means that the best solution is obtained when the $σ_{yz}^{S}$ action zone covers the whole layer. In this case, Hashin’s and Soutis models give identical results. If $b > d$ we obtain that $z_{1} / b < 1$ which means that the $σ_{yz}^{S}$ action zone covers only a part of the 0-layer as was earlier suggested in Ref. 12. The predictions of the Soutis model in this case are more accurate than the results under the assumption that the linear distribution is over the whole 0-layer. It is interesting to note that z₁ depends on the crack density: with increasing crack density, the stress $σ_{yz}^{S}$ action zone tends to expand in the thickness direction to cover the whole 0-layer.

It has to be emphasized that the values of the complementary energy in Table 2 (for low crack density, $ρ_{n} = 0.1$ ) and in Table 3 (for $ρ_{n} = 0.8$ ) can be used only for ranking the models with respect to the accuracy. They should not be used to determine the magnitude of the improvement. It can be done using equation (40) and also comparing the stress distributions as shown in Figure 10 for the case of [0₃/90₂]_s with elastic properties from Table 1 and with normalized crack density $ρ_{n} = 0.5$ . Figure 10(a) shows that the shape of the out-of-plane distribution given by $ϕ_{C} (z)$ is rather different. Nevertheless, the differences in the normalized in-plane shear stress distribution in Figure 10(b) are much smaller and they are almost negligible between Soutis and the present model. Since the value of the average shear stress decrease between cracks determine the amount of shear modulus reduction, the Hashin’s model predicts the largest reduction, the Soutis model predicts smaller reduction and the present model even slightly smaller. FEM profiles are not presented here because they depend on the distance from the symmetry axes and the average value of stress which is of relevance in this paper is not available.

Figure 10.

(a) Shape of the out-of-plane distribution $ϕ_{c} (z)$ and (b) in-plane shear stress distributions.

In the analysis, we assumed that cracks in surface layers are located symmetrically with respect to the middle-plane. Due to statistical failure properties distribution along the 90-layer, more probable is random location, especially for low crack density. For higher crack density, it may be more staggered (the crack in the bottom layer is located in the middle between cracks in the top layer). The variational stress model for the latter case was introduced by Nairn.¹⁷ Since the variation of failure properties along the 90-layer transverse direction is much larger than the stress perturbation in the bottom layer of the laminate due to crack in the top layer, the assumption of staggered cracks and the assumed symmetry of the damage with respect to the mid-plane used in this paper may be two extreme cases.

Conclusions

The developed analytical model for the shear modulus of laminate with intralaminar cracks in 90-layer is based on minimization of the complementary energy. It is shown that this model is more general and accurate than the Hashin’s variational model, which was developed using an assumption that the out-of-plane shear stress has linear distribution across the thickness of the damaged layer and also across the thickness of the constraint layer. The refinement presented here consists of introducing an exponential shape function to characterize this distribution, with the shape parameters determined during minimization. Hashin’s model becomes a particular case of the present model.

As expected from an improved model based on minimization of complementary energy, the predicted laminate shear modulus is slightly higher than according to Hashin’s model. The difference is larger for laminates with relatively thin 90-layers. Laminates with damaged 90-layer in the central part as well as on the surface are analysed. The predicted shear modulus reduction is the same in both cases, when exponential shape function is used. The same conclusion follows from Hashin’s model and it is confirmed by FEM results.

Finally, it is shown that the shear lag model introduced by Soutis¹² with an unknown perturbation zone can be considered as a special case of the developed variational model and the size of the zone may be found. In this case, the shape function in the 0-layer is bilinear and the out-of-plane shear stress is active only in a part of the 0-layer. It is demonstrated that for relatively thick 0-layer, this model is more accurate than Hashin’s model but less accurate than the presented model. For laminates with relatively thin 0-layer, the results of the Hashin’s and Soutis models coincide.

A comparison with FEM results shows that the obtained shear modulus values are conservative and rather insensitive to the refinement of the stress description in the constraint layer. Therefore, improvement of the model requires introducing non-linear shape functions also in the damaged 90-layer.

Footnotes

Conflict of Interest

None declared.

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

References

Tsai

Daniel

(1992) The behaviour of cracked cross-ply composite laminates under shear loading. Int J Solids Struct, 29(24): 3251–3267.

Giannadakis

Mannberg

Joffe

(2011) The sources of inelastic behavior of glassfibre/vinylester non-crimp fabric [45/-45]s laminates. J Reinf Plast Compos, 30(12): 1015–1028.

Hashin Z. Analysis of cracked laminates: a variational approach. Mechanics of Materials 1985; 4: 121–136.

Tan

Nuismer

(1989) A theory for progressive matrix cracking in composite laminates. J Compos Mater, 23: 1029–1047.

Nuismer

Tan

(1988) Constitutive relations of a cracked composite lamina. J Compos Mater, 22: 306–321.

Tsai

Daniel

Lee

. Dvorak

Lagoudas

(1990) Progressive matrix cracking of crossply composite laminates under biaxial loading. Microcracking-induced damage in composites, New York: ASME, pp. 9–18.

Abdelrahman

Nayfeh

(1999) Stress transfer and stiffness reduction in orthogonally cracked laminates. Mech Mater, 31: 303–316.

Henaff-Gardin

Lafarie-Frenot

Gamby

(1996) Doubly periodic matrix cracking in composite laminates Part 1: general in-plane loading. Comp Struct, 36: 113–130.

Kashtalyan

Soutis

(2000) Stiffness degradation in cross-ply laminates damaged by transverse cracking and splitting. Compos Part A, 31: 335–351.

10.

Kashtalyan

Soutis

(2002) Mechanisms of internal damage and their effect on the behavior and properties of cross-ply composite laminates. Int Appl Mech, 38(6): 641–657.

11.

Fan

Zhang

(1993) In situ damage evolution and micro/macro transition for laminated composites. Compos Sci Technol, 47: 107–118.

12.

Zhang

Fan

Soutis

(1992) Analysis of multiple matrix cracking in [±θm/90n]s composite laminates. Part 1. In-plane stiffness properties. Composites, 23(5): 291–304.

13.

Kashtalyan

Soutis

(2005) Analysis of composite laminates with intra- and interlaminar damage. Prog Aerosp Sci, 41: 152–173.

14.

Vinogradov

Hashin

(2010) Variational analysis of angle-ply laminates. Compos Sci Technol, 70: 638–646.

15.

Lundmark

Varna

(2006) Crack face sliding effect on stiffness of laminates with ply cracks. Compos Sci Technol, 66: 1444–1454.

16.

Tao

Sun

(1996) Effect of matrix cracking on stiffness of composite laminates. Mech Comp Mater Struct, 3: 225–239.

17.

Nairn

(1992) The formation and effect of outer-ply microcracks in cross-ply laminates: a variational approach. Eng Fract Mech, 41(2): 203–221.

Potential of a simple variational analysis in predicting shear modulus of laminates with cracks in 90-layers

Abstract

Keywords

Introduction

Repeating element and assumptions

Model formulation

Assumed stress profiles

Expressions for stress components

Cracked central 90-layer

Cracked surface 90-layer

Minimization of the complementary energy

Shear modulus expression

Selection of the shape function ϕ ( z )

Exponential shape function

Central 90-layer with cracks

Surface 90-layer with cracks

Hashin's assumptions

Bilinear shape functions by Soutis

Results and analysis

Minimum of the complementary energy

Comparison with Hashin’s model

Parametric analysis of the [0m/90n]s laminate using exponential shape functions

Comparison with FEM

Accuracy of Soutis et al. model 12

Conclusions

Footnotes

Conflict of Interest

Funding

References

Selection of the shape function $ϕ (z)$

Parametric analysis of the [0_m/90_n]_s laminate using exponential shape functions

Accuracy of Soutis et al. model¹²