Singularities near three-dimensional corners in composite laminates

International Journal of Fracture 115: 361–375, 2002. © 2002 Kluwer Academic Publishers. Printed in the Netherlands.

Singularities near three-dimensional corners in composite laminates A. DIMITROV, H. ANDRÄ and E. SCHNACK Institute of Solid Mechanics, Karlsruhe University 76128 Karlsruhe, Germany Received 19 February 2001; accepted in revised form 8 April 2002 Abstract. The high interlaminar stresses, which appear in laminated composites due to the boundary layer effect near the free edge, play an important role in the analysis and design of advanced structures. Moreover, they are also the dominant effect causing delamination. Even if the singular behavior of such structures is investigated in many works, most of them deal either with 2D, or with pseudo-3D problems, i.e. problems of two variables in a threedimensional space. However, some numerical and experimental findings indicate that laminated plates exhibit a tendency to delaminate at corners, an effect impossible to be determined by a two-dimensional analysis. The aim of the present paper is to investigate stress singularities in a laminated composite wedge under consideration of real three-dimensional corner effects. A weak formulation, as well as a finite element approximation technique introduced in the past for isotropic problems is extended here to cover anisotropic material properties. This formulation leads to a quadratic eigenvalue problem, which is solved iteratively using the Arnoldi method. The first singular terms in the asymptotical expansion of the linear-elastic solution near the vertex of the wedge are obtained as eigenpairs of this eigenvalue problem. The order and mode of singularity are reported for all wedge angles and different fiber orientations for angle and cross-ply laminates. All calculations are based on a typical for some high modulus graphite-epoxy systems orthotropic material model. Key words: Arnoldi method, boundary layer, composite laminates, corner singularity, free edge, laminated composite wedge, quadratic eigenvalue problem

1. Introduction The use of fiber-reinforced composites in various application areas permanently increases during the last years, especially in the aircraft and automotive industries. Features as high stiffness-to-weight and strength-to-weight ratios, especially at temperatures over 2000 ◦ C as in the case of carbon fiber – carbon matrix composites make them attractive for applications in many advanced structures like space vehicles, airplanes, reactors, etc. Nevertheless, the applicability of such materials to primary structures is difficult due to an essential shortage of design data – in the first place an unpleasant disposition to delamination should be named. The delamination as a critical failure mechanism for a wide range of fiber-reinforced composites has received without any doubt appreciable attention. Since the original paper by Pipes and Pagano (1970) various numerical (Rybicki, 1971; Wang and Crossman, 1977; Wang and Dickson, 1978; Kassapoglou, 1986) and experimental (Pipes and Daniel, 1971; Oplinger et al., 1974; Czarnek et al., 1983) investigations have been performed to clarify the origin of this phenomenon. Some high inter-laminar stress concentrations near the free edge (the so called ‘free edge’ or ‘boundary layer’ effects) are now well established as the dominant effect causing delamination.

362 A. Dimitrov et al. However, even if early works (Bogy, 1968; Hein and Erdogan, 1971; Wang and Choi, 1982) have shown the singular nature of stresses at the free edge (in the framework of linear elasticity), sometimes the importance of such studies is underestimated and singularities are considered as a ‘mathematical artifact’, without any relation to the ‘real’ world. This leads unfortunately to some misapprehensions – in a number of cases one tries to calculate numerically the stresses at the edge using ‘suitable elements and grids’ (Altus et al., 1980; Spilker and Chou, 1980; Herakovich et al., 1985) and to validate the obtained results by comparison with experimental findings (Herakovich et al., 1985) or by convergence considerations (Spilker and Chou, 1980). Unfortunately, there is no numerical approximation scheme (except some augmented methods) which would be able to reach the infinite values of the theoretical solution at the edge and thus arbitrary results between a constant and infinity are possible, in dependence on the mesh size. This fact finds a excellent illustration in the controversial results obtained by various investigators solving the free-edge problem, see (Whitcomb et al., 1982) for a survey. It is evident that the type of the singular solution plays an important role – the application of stress intensity factors (SIFs) or energy-based failure criteria is not possible if it is unknown. Moreover, some failure concepts based on SIFs can be used not only for crack problems, but for arbitrary problems containing singularities, as bodies with polyhedral corners or free edges (Groth, 1985, 1988), similar to strategies in the linear-elastic fracture mechanics (LEFM). The singular solution in various problems became therefore the focus of attention in many works, see Leguillon and Sanchez-Palencia (1987), Gu and Belytschko (1994), Pageau and Biggers (1996), Yosibash (1997), and the references therein. However, most of these investigations deal with two-dimensional or pseudo three-dimensional problems. In a typical work the solution is investigated in the vicinity of an edge, but sufficiently away from a vertex, for instance near to point B in Figure 1. However, what happens in the neighborhood of a vertex, i.e. near points A,C,D,E? The list of publications becomes significantly shorter in this case, because some really three-dimensional effects should be taken into account. Rigorous (by means of an asymptotical expansion) studies of this kind are relatively rare to our knowledge, we can mention (Somaratna and Ting, 1986; Ghahremani, 1991; Leguillon and Sanchez-Palencia, 1999). However, some numerical analyses as well as experimental findings (Griffin and Roberts, 1983; Griffin, 1988) indicate that laminated plates exhibit a tendency to delaminate at corners, an effect which cannot be determined by two-dimensional (plane strain or plane stress) analysis. This underlines once more the importance of taking into account all relevant three-dimensional effects. The aim of the present paper is to investigate stress singularities in laminated composites under consideration of real three-dimensional corner effects. A weak formulation, as well as a finite element approximation technique introduced in a previous work (Dimitrov et al., 2000) for isotropic problems is extended here to cover anisotropic material properties. This formulation leads to a quadratic eigenvalue problem, which is solved iteratively by the Arnoldi method. The first singular terms in the asymptotical expansion of the linear-elastic solution near the vertex of the wedge are obtained as eigenpairs of this eigenvalue problem. A description of the used FE-formulation and the iterative solution technique is briefly outlined in Section 2. In Section 3 some results for the laminated composite wedge are reported. The order and mode of singularity are calculated for all wedge angles and different fiber orientations for angle and cross-ply laminates. For all calculations a typical for some graphite-epoxy systems orthotropic material model is used. A conclusion is provided in Section 4.

Singularities near three-dimensional corners in composite laminates 363

Figure 1. Typical singular points in the interface of a laminated composite.

2. Formulation of the problem From a theoretical point of view it is well known that the linear elastic solution may contain gradient singularities, if the domain of consideration includes re-entrant corners like cracks or sudden changes in the material properties as in the case of composites, see for instance Kondratiev (1967), Leguillon and Sanchez-Palencia (1987), Dauge (1988), Grisvard (1992) and the references therein for a survey. For the treatment of such situations it can be helpful to expand the solution in the neighborhood of the singular point in an asymptotical series
Ki r λi fi (θ, ϕ), (1) u= i

where r, θ, ϕ are the polar coordinates, λi are the singularity exponents, fi are the so called angular functions and Ki are the corner stress intensity factors (CSIFs). The asymptotical series (1) can be explicitly constructed, particularly if one considers special geometries and material properties. In this case the exponents λi are obtained as solutions of some transcendent equations. However, for general three-dimensional problems such approaches do not work and some numerical methods are needed to obtain λi and fi . In the following we briefly describe a technique introduced in our previous works (Dimitrov et al., 2000; Dimitrov and Schnack, 2001) for isotropic problems and extended here to anisotropic materials. It is based on a weak formulation and a finite element approximation technique and results in a quadratic eigenvalue problem, solved iteratively. 2.1. W EAK FORMULATION In order to keep our considerations sufficiently general, we introduce an abstract mixed boundary value problem of linear elasticity in a bounded domain ⊂ R3 which coincides with a cone K in the -vicinity UO of the origin O, so that O := K ∩ UO = ∩ UO , see Figure 2. We denote with S the cross-section of the cone on the unit sphere. By this general model every geometry containing a conical singular point can be described (note that the shape of S is not specified), including the problem of a bimaterial wedge. Our aim is to find all solutions, often called the eigenstates, which satisfy the differential equation of elasticity and the boundary conditions in a sufficiently small neighborhood O of O, i.e. on 0 ∪ 1 , but not the boundary conditions remote from O, i.e. on T . 0 , 1 , T denote the Dirichlet, Neumann and transmissional part of the boundary ∂ O , respectively. We obtain the eigenstates by solving the Lamé system on O with only the local boundary

364 A. Dimitrov et al.

Figure 2. Boundary value problem with a singular point O (a). For obtaining the order and mode of the singularity the -neighborhood of O is considered with only local boundary conditions on 0 ∪ 1 (b).

conditions, Lu := D T C D u = 0 on O , u = 0 on 0 , T u := t(u) = 0 on 1 ,

(2)

where D(∂x , ∂y , ∂z ) is the symmetrical gradient operator in matrix notation, C denotes the matrix of elastic moduli and t are the boundary tractions. In other words we consider a local problem with zero body forces and homogeneous boundary conditions. Because (2) is not a full boundary value problem (no conditions are specified on a part of the boundary T , i.e. we consider the solution itself as boundary condition on T ), the requirements of the uniqueness-theorem are not fulfilled and the solution is therefore nonunique. We obtain a set of solutions from which we should exclude these with unbounded elastic energy: We only look for solutions in the Sobolev space [H 1 ( O )]3 in order to have a square-integrable first derivative1 . The unique solution of the problem in the vicinity of O can be obtained afterwards as a linear combination of all eigenstates, so that all boundary conditions including the transmissional one on T to be fulfilled. Since we are interested in a numerical solution, the boundary value problem (2) should be formulated in a weak form. Let us denote the space of admissible (with respect to the boundary conditions) displacement fields by [H01 ( O )]3 . Then, introducing some test function v ∈ [H01 ( O )]3 , we formulate: PROBLEM 1. Find u ∈ [H 1 ( O )]3 , so that B(u, v) = 0, ∀v ∈ [H01 ( O )]3 , where the bilinear form B(u, v) is defined in the linear elasticity by  σ T (u) ε(v) d . B(u, v) :=

(3)

(4)

O

The right-hand side in (3) vanishes because of the homogeneous boundary conditions. The stress and strain vectors, which contain all independent components of the stress and strain 1 We denote by [H m ( )]3 the usual Sobolev space of order m in three dimensions over a domain .

Singularities near three-dimensional corners in composite laminates 365 tensors respectively, are defined by εT := [εx , εy , εz , 2εxy , 2εxz , 2εyz ], σ T := [σx , σy , σz , τxy , τxz , τyz ].

(5) (6)

2.2. F INITE E LEMENT A PPROXIMATION We now consider a finite-element approximation of the weak problem (3). Introducing two different finite-dimensional subspaces Uh  = Vh as trial and test spaces we obtain PROBLEM 2. Find an approximate solution uh ∈ Uh ⊂ [H 1 ( O )]3 , so that B(uh , vh ) = 0,

∀vh ∈ Vh ⊂ [H01 ( O )]3 .

(7)

Such a formulation is called Galerkin–Petrov method. The spaces Uh , Vh are obtained by the partitioning of S into triangular finite elements i . According to (1) the trial and test functions in a typical space sector (r, θ, ϕ) ∈ [0, ] × i can be expressed by uhi (r, θ, ϕ) = r λ N(θ, ϕ) T−1 d di , vhi (r, θ, ϕ) =

(8)

(r) N(θ, ϕ) T−1 d bi ,

ˆ N, ˆ N ˆ and the nodal displacement vector di = with the block-diagonal matrix N = N, i i i i i i T [u1 , v1 , w1 , . . . , uK , vK , wK ] , which contains the displacements in the K nodes of the i-th ˆ = [N1 , N2 , · · · , NK ] collects K shape functions, which can be detriangle. The vector N fined in terms of the usual triangular coordinates (Dimitrov and Schnack, 2001). The Boolean matrix Td reorders the components of di , so that N can be defined in the convenient blockdiagonal form. Equation (8) comprises the main idea of the method: By a semi-continuous approach, in which the discretization is performed only with respect to θ, ϕ, a separation of variables can be enforced. The final problem is consequently independent of r. For the formulation of the problem one needs a constitutive relation between stresses and strains, which is introduced via the generalized Hooke’s law (with respect to the principal material coordinate system) ˆ εˆ (uhi ). σˆ (uhi ) = C

(9)

ˆ depends for a general anisotropic domain on 21 independent The symmetric 6 × 6 matrix C constants. In the important case of an orthotropic material this matrix is given by   C11 C12 C13   C12 C22 C23     C13 C23 C33 , ˆ :=  (10) C   C44     C55 C66

366 A. Dimitrov et al.

Figure 3. Relation between the material (1, 2, 3) and the global (x, y, z) coordinate systems.

where the independent elastic coefficients are reduced to 9 (Reddy, 1996) (Note that the components of the stress and strain vectors have a different order in (Reddy, 1996)). The coefficients Cij can be expressed through the engineering constants Ei , νij and Gij by 1 − ν23 ν32 ,  1 − ν13 ν31 , = E2 

ν21 + ν31 ν23 ,  ν32 + ν12 ν31 , = E2 

ν31 + ν21 ν32 ,  1 − ν12 ν21 , = E3 

C11 = E1

C12 = E1

C13 = E1

C22

C23

C33

C44 = G12 ,

C55 = G13 ,

(11)

C66 = G23 ,

 = 1 − ν12 ν21 − ν23 ν32 − ν31 ν13 − 2ν21 ν32 ν13 , where: – E1 , E2 , E3 are Young’s moduli in the principal material directions 1,2,3, – νij , i, j = 1, 2, 3 are Poisson’s ratios, defined as ratio of the transverse strains −j /i if a stress is applied only in the direction i and all other stress components are zero, – Gij , i, j = 1, 2, 3 are the shear moduli in the i − j planes. The above constitutive relations were written in terms of stresses and strains in the principal material coordinate system. If the global coordinate system doesn’t coincide with the material one, a transformation rule should be established. Let (x, y, z) denotes the global coordinate system and let (1, 2, 3) be the principal material coordinate system, in which the constitutive equations are defined. We consider the following transformation from the global to the material coordinate system: First a rotation around the z-axis on angle ζ , and then a rotation around the y-axis on angle η, see Figure 3. The material matrix C with respect to the global coordinate system can be calculated by ˆ TT (ζ, η), C = T(ζ, η) C

(12)

Singularities near three-dimensional corners in composite laminates 367 where the transformation matrix T(ζ, η) is defined by 

cζ2 cη2

sζ2

cζ2 sη2

−2 sζ cζ cη

2 cζ2 sη cη

−2 sη sζ cζ



   c2 s 2 cζ2 sη2 sζ2 2 sζ cζ cη 2 sζ2 sη cη 2 sη sζ cζ  η ζ     2 2   sη 0 cη 0 −2 sη cη 0   T(ζ, η) =  , 2 2 2 2  cη sζ cζ −sζ cζ sη sζ cζ cη (2cζ − 1) 2 sη cη sζ cζ sη (2cζ − 1)      2   −cζ sη cη 0 cζ sη cη sη sζ cζ (2cη − 1) −cη sζ   2 −sζ sη cη 0 sζ sη cη −sη cζ sζ (2cη − 1) cη cζ

(13)

sζ = sin(ζ ), cζ = cos(ζ ), sη = sin(η), cη = cos(η). So, we obtain with respect to the global coordinate system a new constitutive equation σ (uhi ) = C ε(uhi ).

(14)

Introducing (8) into the expression for the discrete bilinear form (7) and considering the constitutive relation (14) we obtain after some additional transformations and the substitution λ = λ¯ − 1/2 a quadratic eigenvalue problem

(15) P + λ¯ Q + λ¯ 2 R d = 0, where the definition of the matrices P, Q and R can be found in Dimitrov et al. (2001), Dimitrov and Schnack (2002). This problem does not depend on r. From the finite element approximation we therefore obtain an algebraic quadratic eigenvalue problem (15), from which the singular exponents λi can be calculated as eigenvalues and the angular functions fi as eigenvectors. The corresponding linear eigenvalue problem is not symmetric (Dimitrov et al., 2000), thus complex roots are possible. Since there is no a-priori information about real or complex roots, the numerical scheme always searches for complex values. In the most general case the eigenvalues may also be defective, i.e. the algebraical and the geometrical multiplicities may not coincide and some logarithmical terms can appear in the solution, which will be slightly ‘more singular’ than the numerical approximation. A suitable solution technique for (15), based on a linearization procedure and an iterative Arnoldi solver has been proposed in Dimitrov et al. (2001), Dimitrov and Schnack (2002). It is most appropriate for large structured matrices because it requires only 2n · O(k) + O(k 2 ) storage (k  n, n denotes the dimension of P, Q, R) and no explicit knowledge of the corresponding standard eigenvalue problem. So the present method requires only one direct factorization of the matrix P for relatively small systems or an incomplete factorization of P for large systems, as well as few matrix-vector products with Q, R in both cases to find all eigenvalues (λ) ∈ (−0.5, 1.0) as well as the corresponding eigenvectors simultaneously. The interval of interest (−0.5, 1.0) is obvious, due to the restriction that we are only interested in solutions with finite energy ((λ) > −0.5) that have a gradient singularity ((λ) < 1.0). 3. Numerical examples In this section we consider the three-dimensional bimaterial wedge problem, see Figure 4. We introduce a global coordinate system (x, y, z) with origin in the wedge vertex, such that the

368 A. Dimitrov et al.

Figure 4. The bimaterial wedge problem. The angle ζ defines the orientation of the fibers with respect to the x-axis.

material interface lies in the x − y plane. The wedge takes the space between −α/2 and +α/2 measured from the positive x-axis. The numerical scheme discussed above allows the study of non-homogeneous problems with anisotropic material properties. For the calculations in the following we introduce for each layer of the wedge an orthotropic material defined via E1 = 13.8 × 104 MPa, E2 = E3 = 1.45 × 104 MPa, G12 = G13 = G23 = 0.586 × 104 MPa, ν12 = ν13 = ν23 = 0.21,

(16)

which is typical for some high-modulus graphite-epoxy systems. We define the above material constants with respect to a positive principal material coordinate system (1, 2, 3), where 1 denotes the fiber direction, 2 a direction normal to 1 in the interface plane and 3 a direction normal to 1 and 2. In order to study the effect of different fiber orientations in the layers we introduce the angle ζ from the positive x-direction to the fiber direction. This means, that the material coordinate system is rotated about the z-axis on angle ζ . The relations for the transformation of the material matrix from the material to the global coordinate system are given with Equation (12) in the previous Section. Our aim is to investigate the asymptotical solution in the neighborhood of the origin O, i.e. to find the exponents λi and the corresponding angular functions fi . Both of them depend in general on the elastic constants and the geometry. Three different cases will be analyzed in the following: − Angle-play composites ζ = [±30◦ ], [±45◦ ], [±60◦ ] for all angles α ∈ (0, 2π ). − Cross-ply composites ζ = [0◦ /90◦ ], [−30◦ / + 60◦ ] for all angles α ∈ (0, 2π ). − A comparison between the corner case α = π/2 and the edge problem α = π for all fiber orientations [±ζ ], ζ ∈ (0, π/2). 3.1. A NGLE - PLY COMPOSITES Consider a bimaterial wedge with opening angle α, which is composed of two orthotropic layers with material properties defined in (16), see Figure 4. The fiber orientation ζ in the first layer differs from the orientation in the second layer. We investigate first the case +ζ and −ζ in layers 1 and 2 respectively, i.e. we consider the so-called ‘angle-ply composite’ [±ζ ].

Singularities near three-dimensional corners in composite laminates 369 The domain which has to be discretized is (θ, ϕ) ∈ (0, π ) × (−α/2, +α/2). The two subdomains with different material properties are localized at θ ∈ (0, π/2) and θ ∈ (π/2, π ). Meshes with linear elements, which are a-priori refined at (θ, ϕ) ∈ {π/2} × {±α/2} have been used in order to obtain a better approximation due to the additional edge singularity. The following boundary conditions are applied: homogeneous Neumann b.c. at ϕ = ±α/2 and a rigid-body coupling condition at theartificial boundary θ = 0 and θ = π , i.e. u1 = u2 = · · · = um for the displacements and m i fi = 0 for the nodal forces for the m nodes at each of the poles. In order to estimate the discretization error some convergence considerations have been performed. The singular exponents λi have been calculated in a number of cases (exactly for α = {π/2, π, 3π/2, 2π }) for a series of locally refined grids (h-method). For instance, in Figure 5 top the smallest exponent for the case α = π/2 is shown over the degrees of freedom (DOF) of different meshes. If we assume, that 1/DOF is proportional to h2 (h denotes the characteristic mesh size), which is true for sufficiently small h, the straight line in Figure 5 top means a constant convergence rate of order two. This corresponds very well to the theoretical estimation (Apel et al., 2001) and is the best possible for the linear elements used here. Similar results have been obtained for all cases under consideration. The convergence of constant known order means that Richardson-extrapolation for h → 0 will improve the accuracy significantly. Assuming this extrapolated, very accurate solution to be the true one and comparing it with the actual solution for a given mesh an estimation of the discretization error can be obtained, see for instance Figure 5 bottom for the case α = π/2 and α = 3π/2. In order to bound the numerical effort it is assumed that the solution for all other wedge angles α behaves similar. For the following calculations only these meshes have been used, for which the relative error was smaller than 0.5%, which is sufficient for our purposes. The eigenvalues λi for different wedge angles α and the fiber orientation [±ζ ] as a parameter are shown in Figure 6 top. All eigenvalues are found to be real within the numerical accuracy. In the range α ∈ (0, π ) only one eigenvalue λ < 1 exists. Moreover, this value is close to 1 and shows a slow variation with respect to α and the fiber orientation ζ , so that the singularity is relatively weak. With this in mind points A, B, D and E in Figure 1 are not critical for the structure. However, in the range α ∈ (π, 2π ) additional eigenvalues are found. Two of them appear immediately after α = π and are significantly smaller than the eigenvalue near 1. Moreover, the smallest one reaches very quickly its minimum near 0.5, which is the well known value for crack singularity in homogeneous materials. This means that a rectangular hole for instance, or a sharp concave corner (point C in Figure 1) will have a singular behavior just as a crack in a homogeneous material, which is without any doubt a potential danger for such a structure. The fourth eigenvalue appears near α = 1.5π . It is significantly larger than the other two eigenvalues and is therefore not critical for the behavior. As can be seen from Figure 6 top, the fiber orientation plays an important role also – the greater the difference between the fibers, the stronger the singularity. This behavior is valid for two of the three additional exponents in the interval α ∈ (π, 2π ). However, the middle eigenvalue does not depend on the fiber orientation, almost the same values for the exponents have been obtained for each of the considered fiber angles ζ . The three deformation modes, corresponding to the three dominant eigenvalues for α = 1.5π, ζ = ±45◦ are presented in Figure 7 in a x − y view (bottom) and a x − z view (top). Shown are the distorted shapes of the unit sphere, i.e. each point on the unit sphere

370 A. Dimitrov et al.

Figure 5. Convergence of calculated eigenvalues λ of angle-ply laminates [±45◦ ]. Extrapolated solution for α = π/2 (top) and relative errors with respect to the extrapolated value (bottom).

Singularities near three-dimensional corners in composite laminates 371

Figure 6. Eigenvalues λ for different wedge angles α for angle-ply laminates [±ζ ] (top) and cross-ply laminates [−ζ1 / + ζ2 ], ζ1 + ζ2 = π/2 (bottom).

372 A. Dimitrov et al.

Figure 7. Modes of deformation (approximated eigenfunctions) for α = 1.5π, ζ = [±45◦ ]. Mode I*/II*/III* corresponds to λ1 = 0.548, λ2 = 0.673 and λ3 = 0.887, respectively.

has been shifted with the corresponding displacement multiplied by a scaling factor (in order to show more clearly the deformation). In order to avoid a mismatch with the well known crack opening modes I/II/III we denote the present modes with I*/II*/III*. The first mode (or mode I*) seems to be a combination of mode I and III-like opening, as can be seen from Figure 7. A mode I-like opening leads apparently to an out-of-plane deformation (with respect to the x − y plane) due to the unsymmetric layer configuration. Whereas in the case of mode II* a pure mode III-like (out-of-plane) opening can be observed. This explains the fact, that mode II* is independent of the fiber orientation — because of the in-plane variation of the fibers the out-of-plane deformation is not affected. Mode III* seems to be more or less similar to mode II, however it is not symmetric with respect to the x − y plane, because of the different fiber orientation within the layers. 3.2. C ROSS - PLY COMPOSITES Similar calculations as described in the previous Section 3.1 have been performed also for the case of cross-ply composites [0◦ /90◦ ] and [−30◦ / + 60◦ ] with material properties defined in (16). The same domain and boundary conditions are used, so that the convergence considerations in Section 3.1 are still valid. The calculated eigenvalues λ for different wedge angles α are shown in Figure 6 bottom, the results for [±45◦ ] are also given for comparison. Again, all eigenvalues are found to be

Singularities near three-dimensional corners in composite laminates 373

Figure 8. Comparison between the eigenvalues λ for the edge problem (α = π) and the corner problem (α = π/2) for different fiber orientations [±ζ ]. The reference solution (small rectangles) is by Yosibash (1997).

real within the numerical accuracy in the whole range of investigation. As in the case of angleply composites in the range α ∈ (0, π ) only one eigenvalue λ < 1 exists, which varies slowly with respect to α. The singularity is relatively weak. In the range α ∈ (π, 2π ) three additional eigenvalues appear. Moreover, the middle eigenvalue coincides with the solution for angle-ply composites and seems to be really independent of the fiber orientation. For angles α > π but close to π the smallest eigenvalue here is stronger than the one in the angle-ply case, however this behavior changes for larger angles. This means that not only the difference in the fiber orientation between the layers plays an important role, but also the orientation of the fibers within the body. 3.3. C OMPARISON BETWEEN THE EDGE AND THE CORNER SINGULARITY IN ANGLE - PLY COMPOSITES In order to compare the behavior of the solution in the important from a practical point of view case of a corner α = 0.5π (point A in Figure 1) and an edge α = π (point B in Figure 1) some calculations are realized for all fiber orientations ζ ∈ (0, π/2). The results are shown in Figure 8. Because of the refined scale for the solution visualization, meshes with higher density have been used in these calculations. The expected relative error is around 0.15%, see Figure 5. The case α = π is a typical two-dimensional free-edge problem and some very accurate numerical (Yosibash, 1997) and analytical results (Wang and Choi, 1982) exist in the literature. The comparison with Yosibash (1997) (small rectangles in Figure 8) shows a very good agreement, which confirms our convergence considerations.

374 A. Dimitrov et al. Figure 8 shows, that the corner singularity is stronger than the edge one in the whole range of ζ , even if both exponents are close to 1 and the singularity is therefore relatively weak. This important conclusion is confirmed also by other numerical and experimental investigations (Griffin and Roberts, 1983; Griffin, 1988), where it is clearly shown, that laminated plates exhibit a tendency to delaminate at corners. This effect has a real three-dimensional nature and cannot be determined by an two-dimensional (plane strain or plane stress) analysis. 4. Conclusion We have demonstrated a numerical procedure for the asymptotical analysis of the linear-elastic solution in the neighborhood of some three-dimensional singular points. This procedure allows among others the consideration of some multi-material joints with anisotropic properties, so the problem of a composite wedge with two orthotropic layers have been investigated. The results show a strong dependence of the singular exponents on the wedge angle: for wedge angles smaller than π (convex wedges) the singularity is relatively weak, whereas for angles greater than π (concave wedges) the dominant singularity is significantly stronger and reaches quickly its minimum near 0.5. This means, that holes with sharp edges or concave corners are much more dangerous for the composite structures as convex corners or edges. The fiber orientation plays also an important role – in general, the greater the difference in the fiber angles, the stronger the singularity. After comparison of the results for a right-angled convex corner with these for a free edge it can be concluded that the corner exhibits a stronger singular behavior than the free edge for all angle-ply fiber orientations and is therefore more dangerous. This fact underlines once again the necessity to take into account all relevant three-dimensional effects in the analysis and design of composite structures. References Altus, E., Rotem, A. and Shmueli, M. (1980). Free edge effect in angle ply laminates – a new three dimensional finite difference solution. Journal of Composite Materials 14, 21–30. Apel, T., Sändig, A.-M. and Solov’ev, S. (2001). Computation of 3D vertex singularities for linear elasticity: Error estimates for a finite element method on graded meshes. Preprint SFB393/01-33, Technical University Chemnitz, Germany. Bogy, D. (1968). Edge-bounded dissimilar orthogonal elastic wedges under normal and shear loading. Journal of Applied Mechanics 35, 460–466. Czarnek, R., Post, D. and Herakovich, C. (1983). Edge effect in composites by moire interferometry. Experimental Techniques 7, 18–21. Dauge, M. (1988). Elliptic boundary value problems on corner domains. smoothness and asymptotics of solutions, Lecture Notes in Mathematics, Vol. 1341, Springer-Verlag, Berlin. Dimitrov, A. and Schnack, E. (2002). Asymptotical expansion in non-Lipschitzian domains – a numerical approach. Numerical Linear Algebra with Applications (to appear). Dimitrov, A., Andrä, H. and Schnack, E. (2001). Efficient computation of order and mode of corner singularities in 3D-elasticity. International Journal for Numerical Methods in Engineering 52, 805–827. Ghahremani, F. (1991). A numerical variational method for extracting 3d singularities. International Journal of Solids and Structures 27, 1371–1386. Griffin, O. (1988). Three-dimensional thermal stresses in angle-ply composite laminates. Journal of Composite Materials 22, 53–69. Griffin, O. and Roberts, J. (1983). Numerical/experimental correlation of three dimensional thermal stress distributions in graphite/epoxy laminates. Journal of Composite Materials 17, 539–548.

Singularities near three-dimensional corners in composite laminates 375 Grisvard, P. (1992). Singularities in Boundary Value Problems, Recherches en Mathematiques Appliquées. Masson, Springer-Verlag, Berlin. Groth, H. (1985). A method to predict fracture in an adhesively bonded joint. International Journal of Adhesion Adhesives 5, 19–22. Groth, H. (1988). Stress singularities and fracture at interface corners in bonded joint. Int. J. Adhesion Adhesives 8, 107–113. Gu, L. and Belytschko, T. (1994). A numerical study of stress singularities in a two-material wedge. International Journal of Solids and Structures 31, 865–889. Hein, V. and Erdogan, F. (1971). Stress singularities in a two-material wedge. International Journal of Fracture Mechanics 7, 317–330. Herakovich, C. et al. (1985). Free edge strain concentrations in real composite laminates: experimental-theoretical correlation. Journal of Applied Mechanics 52, 787–793. Kassapoglou, C. (1986). An efficient method for the calculation of interlaminar stresses in composite materials. Journal of Applied Mechanics 53, 744–750. Kondratiev, V. (1967). Boundary problems for elliptic equations in domains with conical or angular points. Transactions Moscow Mathematics Society 16, 227–313. Leguillon, D. and Sanchez-Palencia, E. (1987). Computation of Singular Solutions in Elliptic Problems and Elasticity, Masson, Wiley, New York. Leguillon, D. and Sanchez-Palencia, E. (1999). On 3D cracks intersecting a free surface in laminated composites. International Journal of Fracture 99, 25–40. Oplinger, D., Parker, B. and Chiang, F. (1974). Edge effect studies in fiber-reinforced laminates. Experimental Mechanics 14, 347–354. Pageau, S. and Biggers, S. (1996). A finite element approach to three-dimensional singular stress states in anisotropic multi-material wedges and junctions. International Journal of Solids and Structures 33, 33–47. Pipes, R. and Daniel, I. (1971). Moire analysis of the interlaminar shear edge effect in laminated composites. Journal of Composite Materials 5, 255. Pipes, R. and Pagano, N. (1970). Interlaminar stresses in composite laminates under uniform axial extension. Journal of Composite Materials 4, 538–548. Reddy, J. (1996). Mechanics of Laminated Composite Plates: Theory and Analysis, CRC Press, Boca Raton. Rybicki, E. (1971). Approximate three-dimensional solutions for symmetric laminates under inplane loading. Journal of Composite Materials 5, 354–360. Somaratna, N. and Ting, T. (1986). Three-dimensional stress singularities in anisotropic materials and composites. International Journal of Engineering Sciences 24, 1115–1134. Spilker, R. and Chou, S. (1980). Edge effects in symmetric composite laminates: importance of satisfying the traction-free-edge condition. Journal of Composite Materials 14, 2–20. Wang, A. and Crossman, F. (1977). Some new results on edge effect in symmetric composite laminates. Journal of Composite Materials 11, 92–106. Wang, J. and Dickson, J. (1978). Interlaminar Stresses in symmetric composite laminates. Journal of Composite Materials 12, 390–402. Wang, S. and Choi, I. (1982). Boundary-layer effects in composite laminates: Part 1 and 2. Journal of Applied Mechanics 49, 541–560. Whitcomb, J., Raju, I. and Goree, J. (1982). Reliability of the finite element method for calculating free edge stresses in composite laminates. Computers and Structures 15, 23–37. Yosibash, Z. (1997). Computing edge singularities in elastic anisotropic three-dimensional domains. International Journal of Fracture 86, 221–245.