## Abstract

The boundary integral equations for the stretching–bending coupling analysis of general laminated plates are derived in this paper with the aid of the reciprocal theorem of Betti and Raleigh. No restriction is placed on the location of the point load or moment, and so it can be an internal or external point or a smooth or non-smooth boundary point. The fundamental solutions derived from Green’s function for an infinite laminated plate subjected to concentrated forces/moments are presented in the complex matrix form. Through the use of some identities converting complex form to real form, the free term coefficients, which are important for the boundary integral equations, are obtained explicitly and expressed in the real form. Alternative formulae for calculating the free term coefficients are also derived using five rigid body movements. By using the explicit real expressions obtained recently for the Stroh-like formalism of coupled stretching–bending analysis, the free term coefficients are further reduced to the cases of isotropic plates and are then compared with known expressions published in the literature. Since stretching and bending decouple for isotropic plates, comparison is made separately for the in-plane problem and plate bending problem, and it is shown that the boundary integral equations derived in this paper agree with previous results for these two special cases.

## 1. Introduction

The main feature of a boundary integral equation is that only the boundaries of the region being investigated are involved, which therefore leads to much smaller systems of algebraic equations for final solution than other methods such as the finite-element method. Based on the boundary integral equations for different problems, several different boundary elements used in solid mechanics or fluid mechanics have been developed (Brebbia *et al.* 1984). Among these boundary elements, the development of plate bending problems (Stern 1979; Hartmann & Zotemantel 1986; Shi & Bezine 1988) started relatively later than that of two-dimensional or three-dimensional elastostatics (Rizzo 1967; Brebbia *et al.* 1984). For example, even though vast advances have been made in boundary element methods for two- and three-dimensional fracture problems (Li *et al.* 1998; Frangi *et al*. 2002; Rungamornrat & Mear 2008*a*,*b*), relatively few studies have been done on boundary element formulation for cracks in plates subjected to transverse loading or bending moments. As to the stretching–bending coupling analysis of composite laminates, as far as the author is aware, no associated boundary integral equation has ever been derived in the literature. In practical applications, to take advantage of the designable characteristics of composite laminates, there is always the possibility to design an unsymmetrically laminated plate. In that case, the stretching–bending coupling may occur no matter what kind of loading is applied on the laminated plates. In this sense, similar to the advancement of elastostatic and fracture problems in two and three dimensions, providing engineers with an efficient approach other than the popular finite-element method to analyse general composite laminates motivates the present study.

In order to establish a boundary element code for the analysis of general composite laminated plates, the derivation of its related boundary integral equations and associated fundamental solutions become important and necessary, especially when none of them has been presented in the literature. In this paper, with the aid of the reciprocal theorem of Betti and Raleigh (Sokolnikoff 1956), the general boundary integral equations for the coupled stretching–bending analysis of composite laminates are derived. The basic variables considered in the formulation are in-plane displacements or tractions in *x*- and *y*-directions, deflection or effective transverse shear force, normal slope or bending moment and corner force.

To make the boundary integral equations work for the programming of boundary element codes, one of the important tasks is finding the appropriate fundamental solutions, which are associated with Green’s functions, i.e. solutions for problems subjected to concentrated forces and moments at arbitrary points. Many analytical closed-form Green’s functions have been provided in the literature for several different problems such as two- or three-dimensional infinite spaces, half spaces for isotropic materials (Brebbia *et al.* 1984), isotropic plate bending problems (Stern 1979) and two-dimensional infinite spaces, half-spaces, bimaterials, cracks, holes or inclusions for general anisotropic materials (Hwu & Yen 1991, 1993; Ting 1996). For composite laminates with coupled stretching–bending deformation, detailed discussion and solutions have been provided by using the complex variable formulation, such as that by Becker (1995). However, owing to mathematical complexity, most of the solutions presented in the literature left a system of linear algebraic equations to be solved numerically. This is inconvenient when we employ the fundamental solutions to the boundary element formulation for solving practical engineering problems. Recently, by employing Stroh-like formalism for the coupled stretching–bending analysis of composite laminates (Cheng & Reddy 2002, 2005; Hwu 2003; Lu 2004), several Green’s functions have been obtained, such as those for infinite laminates (Cheng & Reddy 2003; Hwu 2004), laminates with holes/cracks (Hwu 2005) and laminates with inclusions (Hwu & Tarn 2006). To establish the boundary integral equations for general laminates with coupled stretching–bending deformation, in this paper the fundamental solutions are derived on the basis of Green’s functions for infinite laminates obtained in Hwu (2004). The success of getting the fundamental solutions for the problem studied in this paper depends heavily on the development of Stroh-like formalism, which is an elegant and powerful complex variable formulation for static problems or steady motions in anisotropic elasticity. One may refer to Fu & Brookes (2006) and Fu (2007) for a more recent interpretation of Stroh-like formalism and its application to the study of edge wave propagation in asymmetrically laminated plates.

With the derived fundamental solutions, it is found that the boundary integral equations should be further modified to treat the singular integrals caused by the point forces/moments and to suit for all possible boundary situations such as corners (i.e. non-smooth boundary points). The associated free term coefficients are then obtained in an explicit real form or alternatively calculated from five different rigid body movements. With the boundary integral equations and the explicit closed-form fundamental solutions derived in this paper, the boundary element for the coupled stretching–bending analysis of composite laminates can then be established and coded in the next step. Since only boundary meshes are needed in the boundary element formulation, it is expected that the boundary element implemented based on the present boundary integral equation is more computational efficient than the standard finite-element method. Moreover, if Green’s functions for specific problems such as holes/cracks (Hwu 2004) or inclusions (Hwu & Tarn 2006) are utilized, the present boundary integral equation can be further improved for the stretching–bending laminated plates with holes/cracks/inclusions.

## 2. Boundary integral equations—internal points

The reciprocal theorem of Betti and Raleigh in terms of stresses and strains can be expressed as
2.1where *σ*_{ij},*ε*_{ij} and , *i*,*j*=1,2,3 are the stresses and strains induced by two different loading systems on the same elastic body whose region is denoted by Ω. If the elastic body is a thin laminated plate with plane region *A*, according to the Kirchhoff’s assumptions, the integration of equation (2.1) with respect to the thickness may lead to
2.2where (*N*_{x},*N*_{y},*N*_{xy}), (*M*_{x},*M*_{y},*M*_{xy}) and (*Q*_{x},*Q*_{y}) are the stress resultants, bending moments and transverse shear forces (figure 1); (*ε*_{x0},*ε*_{y0},*γ*_{xy0}), (*κ*_{x},*κ*_{y},*κ*_{xy}) and (*γ*_{xz},*γ*_{yz}) are the midplane strains, curvatures and transverse shear strains, which are related to the midplane displacements *u*_{0},*v*_{0} and *w* by
2.3

With relations (2.3), the surface integral of equation (2.2) can be reduced to a line integral by taking integration by parts, for example,
2.4in which *Γ* is the boundary of area *A* whose normal direction is denoted by (*n*_{1},*n*_{2}). After taking integration by parts as that shown in equation (2.4), equation (2.2) can be separated into two parts. One is a surface integral and the other is a line integral. To reduce the surface integral further, we now consider the equilibrium equations of the plates, which are
2.5awhere
2.5band *f*_{x},*f*_{y} and *f* are body forces in *x*-, *y*- and *z*-directions. Thus, *q*_{x},*q*_{y},*q* and *m*_{x},*m*_{y} represent the total distributed loads and moments applied on the upper and lower surfaces of the plates including the forces/moments induced by the body forces (figure 1). By employing the equilibrium equation (2.5) and the transformation relations (A1), equation (2.2) can now be reduced to
2.6where the subscripts *n* and *s* denote the values in normal and tangential directions, respectively (figure 2).

For thin plates, transverse shear deformations are usually ignored, i.e. 2.7

With this assumption and the last two equations of (2.3), the rotation angles *β*_{x} and *β*_{y} will not be independent but related to the deflection *w* by
2.8which also leads to
2.9Substituting equation (2.9)_{2} into equation (2.6) and taking integration by parts as that shown in equation (2.4), we get
2.10where *V*_{n} is the effective shear force defined as
2.11and *Γ*_{0} and *Γ*_{f} represent, respectively, the starting and final points of the boundary *Γ*. If the *M*_{ns}*w** or value is continuous, the last terms of both sides of equation (2.10) will vanish. Otherwise, the addition of these two terms becomes necessary, which may occur at corners where the *M*_{ns} value is discontinuous. When the boundary has more than one corner, the last terms of both sides of equation (2.10) may be represented as
2.12where the subscript *k* stands for the value in the *k*th corner, the superscripts + and − denote, respectively, the value ahead of and behind the corner and *N*_{c} is the number of corners (figure 1).

Through the above derivation, it should be noted that the integral equation (2.6) is applicable for relatively thick plates whose transverse shear deformation is taken into account, whereas equation (2.10) is applicable for thin plates in which the effects of transverse shear deformation are ignored. Since the fundamental solutions to be presented in the next section are obtained for thin laminated plates, in the following derivation, we will start from equation (2.10) instead of equation (2.6).

Consider four independent unit point loads or moments applied at point ** ξ** inside the body, e.g. the load applied in each of the three orthogonal directions

**e**

_{i},

*i*=1,2,3, and the moment applied on the surface with normal

**n**, i.e. 2.13where

*δ*(

**−**

*ξ***x**) represents the Dirac delta function,

**is the singular load point and**

*ξ***x**∈

*A*is the field point. Substituting each point load of equation (2.13) independently into equation (2.10) and using equation (2.12), we get 2.14in which new notations are used for the convenience of later presentation. They are 2.15and , and ,

*i*=1,2,3,

*j*=1,2,3,4, represent, respectively,

*u*

_{j},

*t*

_{j}and

*t*

_{c}at point

**x**corresponding to a unit point force acting in the

**e**

_{i}direction applied at point

**, whereas , and ,**

*ξ**j*=1,2,3,4, represent

*u*

_{j},

*t*

_{j}and

*t*

_{c}at point

**x**corresponding to a unit point moment acting on the surface with normal

**n**applied at point

**. The normal**

*ξ***n**associated with the boundary

*Γ*(

**x**) is the normal of the surface boundary at point

**x**, whereas the normal

**n**associated with the interior domain

*A*(

**x**) is the normal compatible with the direction of the distributed moment

*q*

_{4}(

**x**)(=

*m*

_{n}(

**x**)) at point

**x**. For example, if the plate is subjected to a uniform moment expressed by

*m*

_{x}and

*m*

_{y}, an equivalent expression with

*m*

_{s}=0,

*m*

_{n}≠0 should be found to determine the direction

**n**of the applied moment.

Since are the stress resultants corresponding to the unit point load applied at ** ξ** directing in

**e**

_{i}, when

**x**approaches to

**they will become singular. Therefore,**

*ξ**if*

*ξ**goes to the boundary Γ*, the integral given in equation (2.14) should be modified. Detailed derivation of the modified boundary integral equations for smooth or non-smooth boundaries will be given in §4. It should also be noted that the surface integral appearing in the second term of the right-hand side of equation (2.14) will not influence the discretization of the boundary element since its integrand contains only known functions. One is the fundamental solution provided in the next section, and the other is the given surface load

*q*

_{j}. The unknown basic variables of the present paper are

*u*

_{j},

*t*

_{j},

*t*

_{c},

*j*=1,2,3,4, and they do not appear in the surface integral of equation (2.14).

## 3. Fundamental solutions

Consider an infinite laminate subjected to a concentrated force and moment at point as shown in figure 1. The elasticity solution of this problem is generally called *Green’s function* and has been obtained in the explicit closed form (Hwu 2004) using a Stroh-like formulation (Hwu 2003). To facilitate the following discussions, we note that Green’s function is given by
3.1awhere Re stands for the real part of a complex number, **A** and **B** are material eigenvector matrices, **u**_{d} and **ϕ**_{d} are the generalized displacement and stress function vectors defined by
3.1band **f**(*z*) is the complex function vector that depends on the location of the concentrated force/moment and is given by
3.2ain which
3.2bThe angular bracket stands for the diagonal matrix whose components vary according to the subscript *α*, *α*=1,2,3,4, i.e. 〈*f*_{α}〉=diag[*f*_{1},*f*_{2},*f*_{3},*f*_{4}]. *μ*_{α},*α*=1,2,3,4, of equation (3.2b) are the material eigenvalues that have been assumed to be distinct and have positive imaginary part; *x*_{1} and *x*_{2} are the alternative expressions of *x* and *y*. In equation (3.1b), *u*_{1} and *u*_{2} are the midplane displacements in *x*- and *y*-directions; *β*_{1} and *β*_{2} are the alternative expressions of *β*_{x} and *β*_{y}, which are related to deflection *w* by equation (2.8); *ϕ*_{1},*ϕ*_{2} and *ψ*_{1},*ψ*_{2} are the stress functions related to the stress resultants *N*_{ij}, shear forces *Q*_{i}, effective shear forces *V*_{i} and bending moments *M*_{ij} by
3.3In the tangent–normal (*s*–*n*) coordinate system, the counterparts of the above relations are (Hwu 2004)
3.4awhere
3.4band *θ* is the angle directed counterclockwise from the positive *x*_{1}-axis to the tangential direction **s** (figure 2). Note that when *θ*=0, equation (3.4) reduces to equation (3.3).

From equations (2.14) and (2.15), we know that the basic functions for the boundary integral equations are *u*_{0},*v*_{0},*w*,*β*_{n} and *T*_{x},*T*_{y},*V*_{n},*M*_{n}, which are different from the generalized displacement and stress function vectors **u**_{d} and **ϕ**_{d} shown in equation (3.1b). Hence, to get the fundamental solutions , and , we need some relations to transform the basic functions between Stroh-like formalism and boundary integral equation. With the aid of relations (3.3) and (3.4), some useful transformation relations are obtained in equations (A2) and (A3). The fundamental solutions can then be obtained by using equations (A2) and (A3) and Green’s functions (3.1) and (3.2), which are
3.5awhere
3.5band
3.5c

In the above, the subscripts *p*, *b*, *t* and *c* denote, respectively, the values related to the in-plane, bending, twisting and corner responses. The superscripts *T*, + and − signify matrix transposition and evaluation immediately ahead of and behind the corner, respectively. The normal and tangential vectors, **n** and **s**, will depend on the location of **x** or ** ξ**. A prime and a superimposed tilde signify differentiation and integration, respectively. For example,
3.6Note that in deriving the fundamental solutions (3.5), the relations obtained in equation (A3) play an important role especially when the normal and tangential directions of the boundary are allowed to vary point by point, i.e.

*s*′

_{1}and

*s*′

_{2}shown in equation (A5) are not zero for general curvilinear boundaries. Without relations (A3), the fundamental solutions obtained by the conventional approach may become awkward and lengthy. Similar to the calculation of the stress components by the conventional approach, one usually applies the transformation law of second-order tensors to obtain the bending moments (

*M*

_{n},

*M*

_{s},

*M*

_{ns}) and in-plane forces (

*N*

_{n},

*N*

_{s},

*N*

_{ns}) and that of first-order tensors to obtain the transverse shear force (

*Q*

_{n},

*Q*

_{s}) and then uses the definition of the effective transverse shear force to obtain (

*V*

_{n},

*V*

_{s}).

## 4. Boundary integral equations–boundary points

The boundary integral equation shown in (2.14) is valid only when the unit point force or moment is applied at a point ** ξ** inside the body. In order to derive a general boundary integral equation valid for any location such as internal, external, boundary or corner points, we now place the source point

**on the boundary and for generality suppose that**

*ξ***is a corner point with an internal angle Δ**

*ξ**φ*(=

*φ*

^{+}−

*φ*

^{−}) (figure 3); should

**be on a smooth boundary, we merely set Δ**

*ξ**φ*=

*π*.

When the source point ** ξ** is placed on the boundary, the fundamental solutions will become singular when

**x**=

**. However, even if when**

*ξ***x**→

**, its associated integral in equation (2.14) may not be infinite in the sense of Cauchy principal value and can be written as 4.1In equation (4.1), the integral path**

*ξ**Γ*(=

*Γ*

_{ρ}+

*Γ*

_{*}) of the left-hand side is replaced by of the right-hand side, where the paths , and are shown in figure 3. Through the new path , the source point

**is located outside the body, and hence the first term of the left-hand side of equation (2.14) vanishes. Moreover, the corner point**

*ξ***also disappears and is replaced by two new corners**

*ξ*

*ξ*^{+}and

*ξ*^{−}. Thus, the corner number shown in the boundary integral equation (2.14) should be modified to , which excludes the corner point

**if it happens to be the source point, whereas the contribution of the new corners**

*ξ*

*ξ*^{+}and

*ξ*^{−}will be assigned to the coefficient of free term

*u*

_{i}(

**).**

*ξ*Through the replacement of the integral path shown in equation (4.1), the boundary integral equation (2.14) can now be modified as
4.2awhere
4.2bNote that the first integral on the left-hand side of equation (4.2a), which comes from the second integral on the right-hand side of equation (4.1), is taken in the sense of Cauchy principal value and represented by . The remaining integrals in equation (4.2a) do not contain any singularities and can be interpreted in the normal sense of integration. Furthermore, from the fundamental solutions provided in equation (3.5), we see that is proportional to 1/*ρ*^{2} when *ρ*→0, where *ρ* is the distance from the integration point **x** to the source point ** ξ**, and hence Δ

_{4}(

**) will not converge. To remedy this situation, like the plate bending problems discussed in Stern (1979), the fourth equations of (4.2a,**

*ξ**b*) should be modified by the replacement of

*u*

_{3}(

**x**) with

*u*

_{3}(

**x**)−

*u*

_{3}(

**), i.e. 4.3awhere 4.3bIn order to evaluate Δ**

*ξ*_{i}(

**),**

*ξ**i*=1, 2, 3, 4 from equations (4.2b) and (4.3b), we assume that is a circular path with radius

*ρ*and polar angle

*φ*, which starts from

*φ*

^{−}and ends at

*φ*

^{+}(figure 3). Along this path, the fundamental solutions (3.5) can be written as (appendix B) 4.4where the explicit expressions of and are given in equation (B3). Substituting equation (4.4) into equations (4.2b) and (4.3b) and using equation (B2) with

*ρ*→0, we get 4.5awhere 4.5band 4.5c

4.5dBy employing the results of (4.5) in (4.2) and (4.3), the modified boundary integral equations can now be written as
4.6in which the first two terms of the left-hand side of equation (4.6) can be expressed in the matrix form as
4.7Note that the extra term *g*_{43} related to Δ_{4}(** ξ**), which comes from the terms associated with

*u*

_{3}(

**) in (4.3a), is 4.8In the next section, we will prove that**

*ξ**c*

_{i5}(

**)≠0,**

*ξ**i*=1,2,3,4 only when the source point

**is a corner point. In other words, the corner points possess an extra degree of freedom**

*ξ**β*

_{s}(

**) in contrast to all the other points, which is not convenient for boundary element formulation. Since in boundary element formulation one corner is usually represented by two nodes, it is suggested that**

*ξ**β*

_{s}(

**) and**

*ξ**β*

_{n}(

**) of the corner points be transformed to the normal slopes of the boundaries ahead of and behind the corner point**

*ξ***,**

*ξ**β*

_{n+}(

**) and**

*ξ**β*

_{n−}(

**). Their transformation can be done by using the following relations 4.9which leads to 4.10Alternatively, one may also consider the transformation of**

*ξ**β*

_{s}(

**) and**

*ξ**β*

_{n}(

**) to the global**

*ξ**x*–

*y*coordinate,

*β*

_{x}(

**) and**

*ξ**β*

_{y}(

**), which can be done by the following transformation law, 4.11**

*ξ*## 5. Free term coefficients *c*_{ij}(*ξ*)

*ξ*

The free term coefficients *c*_{ij}(** ξ**) of the modified boundary integral equations have been provided in equation (4.7) when the nodal degrees of freedom are set to be (

*u*

_{1},

*u*

_{2},

*u*

_{3},

*β*

_{n},

*β*

_{s}). Substituting equation (4.10) or (4.11) into equation (4.5a), the coefficients

*c*

_{ij}(

**) for the alternative degrees of freedom can be obtained as follows: 5.1awhere 5.1band 5.1c**

*ξ*### (a) Explicit real form solutions

Although the closed-form solutions for the coefficients *c*_{ij}(** ξ**) have been obtained in equation (4.7), their calculation depends on the integrals of complex functions given in equation (B3). For the in-plane or pure bending problem associated with an isotropic plate, it is well known that

*c*

_{ij}=

*δ*

_{ij}/2 for a smooth boundary,

*c*

_{ij}=

*δ*

_{ij}for an internal point, in which

*δ*

_{ij}is the Kronecker delta, and

*c*

_{ij}=0 for a point outside the body (Brebbia

*et al.*1984). Therefore, with the knowledge of the results for the special cases—isotropic plates, it is of interest to derive explicit real form solutions of the coefficients

*c*

_{ij}(

**) for the general coupled stretching–bending problems.**

*ξ*In order to get the explicit real form solutions, some identities converting the complex form to real form become important. Because the Stroh-like formalism developed by Hwu (2003) was purposely organized into the form of Stroh formalism (Ting 1996), all the identities established for two-dimensional problems can be applied to the present coupling case without any further rigorous proof. With this understanding, for the convenience of readers’ reference, several identities that are useful for the present derivation are collected in appendix C. Moreover, some explicit expressions for coupled stretching–bending problems, which are useful for the reduction of our present results to the special case—isotropic plates, can be found in Hwu (in press).

Using the identities shown in appendix C, the real form solutions of of (B 3*a*–*c*) for *general composite laminates*, which are related to the fundamental solutions , can now be obtained as
5.2a
5.2b
5.2cwhere
5.2dIn (5.2b,*c*), *S*_{ij}, *X*_{ij}, *X*_{ij}(*φ*), , *X*_{ij}(*φ*), *Y*_{ij}(*φ*), , *i*,*j*=1,2,3,4, are the components of the real matrices defined in equations (C3)–(C6).

Note that the complex form solutions shown in equation (B3) are valid for the laminates whose material eigenvector matrices **A** and **B** exist, i.e. the material eigenvalues *μ*_{α} are distinct or a complete set of independent eigenvectors exists when *μ*_{α} are repeated. For degenerate laminates that have repeated material eigenvalues but do not have a complete set of eigenvectors, a small perturbation of the material constants is usually made for the numerical calculation to yield a complete set of independent material eigenvectors. In contrast, when we apply the real form solutions obtained in equation (5.2), the calculation of the material eigenvalues and eigenvectors has been avoided. Hence, the explicit real form solutions (5.2) are valid for any kind of laminated plates including degenerate laminates such as isotropic plates discussed next.

### (b) Isotropic plates

By using the explicit expressions shown in Hwu (in press) for isotropic plates, the results of (5.2a–*c*) can be further reduced to
5.3a
5.3bSubstituting equation (5.3) into (4.5c,*d*), we get
5.4a
5.4b
5.4cin which *k*_{11}, *k*_{12}, *k*_{21} and *k*_{22} agree with the coefficients for plane stress problems (Hartmann 1983). Note that Hartmann’s solution was presented for the plane strain case, which should be converted to the plane stress case with *ν* replaced by *ν*/(1+*ν*). By equation (4.5b), we have
5.5which *broadly* agree with the solutions presented by Stern (1979) for isotropic plates. Here, the word ‘broadly’ is used because Green’s functions used in our problem are *not* exactly the same as those used by Stern. When reduced to the isotropic plates with bending rigidity *D* in terms of the polar coordinates *r* and *φ*, it is found that the deflection *w* associated with the transverse load is (Hwu in press)
5.6which is different from given by Stern (1979). Moreover, the fundamental solution of the deflection *w* associated with the bending moment , or say , is (Hwu in press)
5.7which is different from given by Stern (1979). Because the closed contour integrals of the stress functions associated with the last two terms of equation (5.6) and the second term of equation (5.7) vanish, with or without these terms will not influence the satisfaction of force equilibrium condition and all the other conditions for solving Green’s functions (Hwu 2004), the fundamental solutions used in this paper as well as those used by Stern (1979) should all be applicable to their corresponding problems.

In addition to noticing the difference in the expressions for the Green’s functions, when making comparison with Stern’s solutions, the coordinates should also be adjusted. This will lead *φ*^{−}−*α* to (*π*/2)−* γ* and

*φ*

^{+}−

*α*to

*kπ*+(

*π*/2)−

*, where*

*γ**k*and

*are the symbols used in Stern (1979).*

*γ*From figure 3, we see that at point *ξ*^{+}, the polar angle *φ*^{+} and the direction angle *α*^{+} are related by *φ*^{+}=*α*^{+}(in some other cases they may be related by *φ*^{+}=*α*^{+}−2*π*), whereas at point *ξ*^{−}, *φ*^{−}=*α*^{−}−*π*. If the point moment is applied in the direction *α*=(*α*^{+}+*α*^{−})/2, it can be proved that . Together with the results obtained in the last equation of (5.4b), we see that all the coefficients related to the extra degree of freedom *β*_{s}(** ξ**) are zero for isotropic plates.

For a smooth boundary, we have *φ*^{+}=*φ*^{−}+*π*. With this relation, the coefficients obtained in equations (5.4) and (5.5) can be further reduced, which leads to the well-known result that *c*_{ij}=*δ*_{ij}/2. Similarly, *c*_{ij}=*δ*_{ij} for an internal point and *c*_{ij}=0 for an external point can also be established by letting *φ*^{+}=*φ*^{−} for a closed contour.

Knowing the coefficients *c*_{ij}(** ξ**) for the nodal degrees of freedom (

*u*

_{1},

*u*

_{2},

*u*

_{3},

*β*

_{n},

*β*

_{s}), the coefficients

*c*

_{ij}(

**) for the alternative degrees of freedom (**

*ξ**u*

_{1},

*u*

_{2},

*u*

_{3},

*β*

_{n−},

*β*

_{n+}) and (

*u*

_{1},

*u*

_{2},

*u*

_{3},

*β*

_{x},

*β*

_{y}) can then be obtained by substituting equations (5.4) and (5.5) into equations (5.1b) and (5.1c). We have 5.8a 5.8bWith

*α*being the direction angle of the source point

**whose applied moment is in the opposite tangent direction (figure 2), it can be proved that the results of equation (5.8b) are identical to those presented in Hartmann & Zotemantel (1986). Note that in Hartmann & Zotemantel (1986), a typing error occurs in the last term of their equation for in which should be corrected as .**

*ξ*If the point moment is applied in the direction *α*=(*α*^{+}+*α*^{−})/2, it can be proved that
5.9which also *conceptually* agrees with the solutions presented by Stern (1979).

### (c) Calculated from rigid body motion

Usually, in boundary element formulation even though no explicit closed-form solutions are provided for the coefficients *c*_{ij}(** ξ**), they still can be computed by applying the boundary integral equations to represent rigid body movements. With this understanding, we now consider the following five different rigid body movements:
5.10Since

*β*

_{n}and

*β*

_{s}defined in equation (2.9) are, respectively, the negative slopes in the normal and tangent directions, they are dependent on the direction of the boundary. Through vector transformation, we have (figure 4) 5.11Substituting the conditions set in the last two rigid body movements of equation (5.10) into the first relation of (5.11) and integrating with respect to

*x*

_{1}and

*x*

_{2}, we can get

*u*

_{3}(

**x**) for the whole plate. Also, substituting equation (5.10)

_{4,5}into equation (5.11)

_{2}, we can get

*u*

_{4}(

**x**), i.e.

*β*

_{n*}(

**x**) for all

**x**along the boundary

*Γ*. Their final expressions are 5.12awhere 5.12bKnowing that

*t*

_{j}=

*q*

_{j}=

*t*

_{c}=0,

*j*=1,2,3,4, for rigid body movements, the substitution of equations (5.10) and (5.12) into the modified boundary integral equations (4.6) will then give us 5.13After getting the coefficients

*c*

_{ij}(

**),**

*ξ**i*=1,2,3,4,

*j*=1,2,3,4,5, for the modified boundary integral equations (4.6), the alternative coefficients for different degrees of freedom shown in equation (5.1a) can then be calculated through the relations shown in (5.1b,

*c*).

From the results shown above, we see that the free term coefficients *c*_{ij}(** ξ**) can be calculated from equation (4.5) or (5.13). The former is an integral around the source point

**whose integrands**

*ξ**t*

_{ij}(

*φ*) have been obtained explicitly in equation (B3) with the complex form and in equation (5.2) with the real form, whereas the latter is an integral around the whole boundary of the body and is obtained from the application of five different rigid body movements. From the computational viewpoint, equation (4.5) is more efficient and accurate than equation (5.13). The relations shown in equation (5.13) provide a useful check on the boundary element coding especially when the fundamental solutions (3.5) obtained in this paper are expressed in terms of multivalued logarithmic functions of complex variables.

## 6. Conclusions

Boundary integral equations suitable for coupled stretching–bending analysis of general composite laminates are derived in this paper. The unknown variables associated with the integral equations are in-plane displacements or tractions in *x*- and *y*-directions, deflection or effective transverse shear force, normal slope or bending moment and corner force. The fundamental solutions derived from the Green’s function for infinite laminates subjected to concentrated forces/moments are presented in the complex matrix form. Through the use of some identities of Stroh-like formalism converting complex form to real form, the free term coefficients, which are important for the integral equations, are obtained explicitly and expressed in the real form. Moreover, an alternative approach using five rigid body movements is also presented in this paper, which can then be used as a numerical check for the free term coefficients. Owing to the involvement of the corner forces and normal slopes, the relations (5.13) obtained from five rigid body movements are slightly more complicated than those for the classical two- or three-dimensional elastostatic problems.

From the modified boundary integral equations shown in (4.6), we see that the corner points possess an extra degree of freedom *β*_{s} in contrast to the other points. To eliminate this extra degree of freedom, a transformation based on the normal slopes of the boundaries ahead of and behind the corner points is suggested in equation (4.10). Since the explicit expressions of the free term coefficients are obtained for a general laminate here for the first time, the comparison can only be made with the special case—isotropic plates, in which in-plane problem and plate bending problem decouple. Under this specialization no matter what kind of degrees of freedom is considered, the explicit expressions of the free term coefficients are shown to agree with those presented in the literature.

## Acknowledgements

The author would like to thank the National Science Council, Taiwan, Republic of China for support through grant NSC 96-2221-E-006-174.

## Appendix A: some transformation relations

*(i) Relations between quantities in x–y and s–n coordinates*
A1*(ii) Relations between the basic functions of Stroh-like formalism and boundary integral equations*

Displacement and slope:A2Stress resultants, bending moments and corner force:

The following relations are obtained from equations (2.15), (3.3) and (3.4).A3In the above, *s*_{1} and *s*_{2} are the first and second components of tangent vector **s**, and *n*_{1} and *n*_{2} are the first and second components of normal vector **n**. They are related to the direction angle *θ* by that shown in equation (3.4b)_{2,3} and have the following relations:A4As to *s*′_{1} and *s*′_{2}, they can be calculated byA5where *R* is the radius of curvature of the considered point. If the point is located on a straight boundary, and *s*′_{1}=*s*′_{2}=0.

## Appendix B: derivation of and

Along the circular path (figure 3),B1andB2where *α* is the value of direction angle *θ* of the source point ** ξ**. If the source point is a corner,

*α*=(

*α*

^{+}+

*α*

^{−})/2, where

*α*

^{+}and

*α*

^{−}are, respectively, the direction angles of the boundary ahead of and behind the corner point

**. In figure 3, the polar angle**

*ξ**φ*of point

**x**on route is related to its direction angle

*θ*by

*φ*+(

*π*/2)=

*θ*. Substituting equation (B1) into equation (3.5), we get the results of (4.4), in whichB3a

B3b

B3candB3d

## Appendix C: identities converting complex form to real form

The following are some identities expressing various combinations of **A**, **B** and *μ*_{α} in terms of real matrices. Since most of them can be found in the book (Ting 1996) on anisotropic elasticity, we list them below without giving their detailed derivations.C1whereC2andC3andC4In the above, **N**_{i} and **N**_{i}(*φ*), *i*=1,2,3, are, respectively, the fundamental elasticity matrices and the generalized fundamental elasticity matrices, which are related by **N**_{i}(0)=**N**_{i}; **S**, **H** and **L** are Barnett–Lothe tensors, which are the matrices defined byC5**S, H** and **L** defined by equation (C5) have been proved to be real and are related to the fundamental elasticity matrices byC6

## Footnotes

- Received August 16, 2009.
- Accepted November 3, 2009.

- © 2009 The Royal Society