JCg^ United States Department of Commerce

National Institute of Standards and Technology

NAT'L INST. OF STAND 4 TECH R.I.C.

A111D3 5 M M 3 Mb

PUBLICATIONS

NIST Special Publication 802

Effect ofa Free Surface on Stress

Distribution in a Bimaterial Composite

V. K. Tewary and R. D. Kriz

QC

100

.U57

#802

1991

C.2
NATIONAL INSTITUTE OF STANDARDS &

TECHNOLOGY

Research Information Center

Qethersburg, MD 20899
091

mi

NIST Special Publication 802

Effect ofa Free Surface on Stress

Distribution in a Bimaterial Composite

V. K. Tewary* and R. D. Kriz**

Materials Reliability Division

Materials Science and Engineering Laboratory

National Institute of Standards

and Technology

Boulder, Colorado 80303-3328

* On attachment from Ohio State University, Columbus, OH

**Present address: Dept. of Engineering Science and Mechanics,

Virginia Polytechnic Institute and State University, Blacksburg, VA 24061

Issued March 1991

U.S. Department of Commerce

Robert A. Mosbacher, Secretary

National Institute of Standards and Technology

John W. Lyons, Director
National Institute of Standards U.S. Government Printing Office

and Technology Washington: 1991

Special Publication 802

Natl. Inst. Stand. Technol.

Spec. Publ. 802

185 pages (March 1991)

CODEN: NSPUE2

For sale by the Superintendent

of Documents

U.S. Government Printing Office

Washington, DC 20402
CONTENTS

Page

Preface v

Abstract xi

PART I

Elastic Green's Function for a Composite Solid Containing

a Free Surface Normal to the Interface

1. Introduction 3

2. Calculation of the SI Green's Function 3

3. Series Representation of the H-Function 19

4. Expressions for the Green's Functions 37

5. Application to the Antiplane-Strain Problem 42

6. Application to a Cubic Crystal Containing a Z-5

Grain Boundary—Plane-Strain Problem 53

7. Summary of Main Results 57

Appendix I-A. Elastic Green's Function for an Infinite

Bimaterial Composite 58

Appendix I-B. Solution of Generalized Inhomogeneous Vector

Hilbert Problem 62

Appendix I-C. Evaluation of Integrals 72

Figures (Part I) 83-93

Part II

Generalized Plane-Strain Analysis of a Bimaterial Composite

Containing a Free Surface Normal to the Interface

1. Introduction 97

2. Displacement and the Stress Field in an Infinite

Composite Solid Subjected to Generalized Plane Strain 98

3 . Displacement and the Stress Field in a

Composite Solid Containing a Free Surface 107
4. Application to E-5 Grain Boundary in a Cubic Crystal 126

5. Application to Fiber-Reinforced Laminated Composites 130

6. Summary of Results and Discussion 137

Appendix II-A. Transform of a Constant and Some Numerical

Considerations 141

Figures (Part II) 145-174

Acknowledgments 175

References (For Parts I and II) 175

iv
PREFACE

The elastic behavior of composite materials is a subject of

strong topical interest because of the technological applications

of these materials. Since no real solid is infinite, it is

important to study the effect of a free surface on the elastic

properties of composite materials. A topic of special interest is

the stress distribution due to an applied load in such solids.

Many papers have already been published on this topic. The

behavior of stress near the intersection of the free surface and

the interface in a bimaterial composite is singular and quite

complex (see, for example, [1-14] and other papers quoted in these

references) . This is usually referred to as the free-edge effect

in composite materials.

The analysis of the free-edge effect is particularly

important because damage in a composite material initiates at the

free-edges, often as ply cracks and delaminations . Usually, these

delaminations are the source of compressive failure in reverse

fatigue loading [15]. The free-edge effect is a characteristic

property of composite materials and, in general, not exhibited by

homogeneous solids. Such a distinction between homogeneous and

composite solids should not be surprising because a similar

distinction is also found in the behavior of the stress near an

interfacial crack tip. Unlike that in a homogeneous solid, the

stress near a crack tip, shows strong oscillations [16,17].

There is some confusion in the literature regarding the

nature of singularities in the stress field near the free surface.

In general it has been found by using different numerical

techniques that the stress behaves like r near r=0 where r is the

radius vector and 8 is a negative number—usually a fraction.

v
Zwiers et al. [10] have shown that in case of an out-of-plane

load (the generalized plane-strain problem) , logarithmic

singularities are also present in the stress field. Further, it

has been pointed out [10] that additional singularities like

[log(r)] n(n-1, 2 , 3 , . . . ) may also arise in the stress field.

Moreover, the behavior of the nonsingular terms has not been

properly investigated so far except numerically. There is clearly

a need for a reliable complete solution which gives the stress

distribution in the composite solid containing a free surface.

The stress distribution as well as many other elastic

properties of a solid can be easily calculated by using the

elastic Green's function of the solid. The Green's function gives

the solution of the elastic equilibrium equations for a unit force

subject to all the prescribed boundary conditions. Once the

Green's function is known, the solution of the elastic equilibrium

equations can then be obtained for any arbitrary force

distribution by a straightforward integration. The Green's

function is therefore a very powerful mathematical tool for

calculating the elastic response of a solid.

The objective of the present work is to calculate the elastic

Green's function for an anisotropic composite solid containing a

free surface and apply it to calculate the stress distribution in

the solid to an applied load. We shall consider a semi-infinite

bimaterial composite solid containing a planar interface and a

single free surface normal to the interface. The solid is assumed

to extend to infinity in other directions. This problem will be

referred to mathematically as the surface-interface problem and

the corresponding Green's function as the surface-interface

(abbreviated as the SI) Green's function.

We shall calculate the SI Green's functions for the

plane-strain and the generalized plane-strain problems. The

out-of-plane components of the strain tensor are zero in the

vi
plane-strain case and are nonzero but constant in the generalized

plane-strain case. The out-of-plane components of the displacement

and the stress tensor will be, in general, nonzero, in both cases.

The SI Green's function as calculated here is applicable to

different types of bimaterial composites such as, for example,

those containing different materials separated by a phase boundary

or different orientations of the same crystal separated by a grain

boundary. It can also be applied to laminated fiber composites in

which different layers have different fiber orientations.

Normally, the Green's function method is used to calculate

the stress distribution in a solid due to an specified force. We

have, however, extended the Green's function method to the case

when instead of a force, an out-of-plane displacement field is

specified. This enables us to study the technologically relevant

generalized plane-strain problems.

In two earlier papers [18,19], henceforth referred to as I

and II, we have calculated the elastic plane-strain Green's

function for an infinite bimaterial composite solid. This Green's

function gives the solution of the elastic equilibrium equations

and satisfies the continuity conditions at the interface. In this

paper we shall use this Green's function to calculate the SI

Green's function which, in addition, will also satisfy the free

surface condition at the free surface.

The calculation of the SI Green's function would require the

solution of an inhomogeneous generalized vector Hilbert problem.

The properties and the solution of the corresponding scalar

problem have been discussed in the excellent treatise by

Muskhelishvili [20] and the vector problem by Vekua [21] and also

briefly in [20] . However, in the present case, we find it more

convenient to solve the aforementioned Hilbert problem by an

alternative method using a complex transform. This gives the

vii
result in a closed integral form which can be calculated quite

easily. This method should be generally useful in several problems

in stress analysis in composite materials containing cracks,

interfaces, and free surfaces.

The earlier methods [1-14] for calculation of the stress

distribution near a free surface in a laminated composite are

either purely numerical or use a series representation. The series

representation is valid only near the origin. Zwiers et al. [10]

have used the Stroh's method [22] for obtaining a series

representation of the stress in terms of the roots of a 12 x 12

complex matrix. All these calculations were done for a solid

subjected to certain specific loadings.

The purely numerical methods are not reliable near the origin

where the stress is singular and quite complex. As pointed out in

[10], a series representation obtained by using an eigenfunction

expansion method [7] or Lekhnitskii method suffers from lack of

completeness of eigenfunction expansion as well as numerical

convergence. Moreover, some of these methods do not account for

out-of-plane components of the displacement and the stress tensor

which, in general, will exist even in plane-strain problems.

The Green's function method can be applied to solids

subjected to any loading or force distribution. In addition to the

stress distribution, it can be used to calculate the interaction

between various defects such as dislocations, surfaces, and

cracks. The interaction between defects determines several

important physical properties of materials such as, for example,

the work hardening of materials.

The Green's function method accounts for the out-of-plane

components of the displacement field and the stress tensor in the

plane-strain problems. It gives the result in a closed integral

• • • Vlll
representation which can be evaluated numerically quite easily.

The result is therefore valid for the entire solid and not only

near the origin. The integral for the Green's function can also be

evaluated analytically, which is useful in identifying the exact

nature of all the singularities which may be present in the

stress. This method is also numerically more convenient than other

methods because it involves manipulation of only a 6 x 6 complex

matrix and not 12 x 12

This paper is in two parts. In Part I, we have calculated the

SI Green's function and applied it to plane-strain problems in

cubic solids. Specifically, we have calculated the stress

distribution due to a line load in a cubic solid containing a Z 5

grain boundary and a free surface normal to the grain boundary

interface. Both the plane-strain and antiplane-strain problems

have been considered. In Part II, we have described an extension

of the Green's function method to the generalized-plane-strain

problem which enables us to calculate the elastic response of a

solid to a prescribed out-of -plane strain. We have applied it to

calculate the stress distribution in a composite solid subjected

to an out-of-plane load.

Our final results are in the form of an analytical closed

integral representation which can be calculated numerically as

well as analytically. The analytical evaluation of the integral by

using the contour integration method gives the result in the form

of a series. This series contains singular as well as non-singular

terms. In general, it is difficult to ascertain the convergence of

the series representation. The numerical evaluation of the

integral does not suffer from such problems and gives the result

to arbitrary accuracy.

The analytical series representation can be used to identify

precisely the nature and the weight of each singularity in the

stress distribution. This is useful in developing a numerical

ix
finite-element technique for the stress analysis in a real solid.

A purely numerical finite-element technique, when applied to

composite materials containing free surfaces, suffers from poor

convergence near the origin because of the singularities in the

stress. A knowledge of the analytical behavior of the stress near

the origin can be used to generate the so-called enriched elements

near the origin, which substantially improves the convergence of

the numerical technique [14].

x
ABSTRACT

This paper is in two parts. In part I, the elastic

plane-strain Green's function is calculated for an anisotropic

bimaterial composite solid containing a free surface normal to the

interface. An exact integral representation is obtained for the

Green's function which is evaluated numerically. The integral is

also evaluated analytically, which gives a series representation

for the Green's function. For illustration, the formalism is

applied . to the antiplane-strain problem and the plane-strain

problem in a cubic solid containing a Z-5 grain boundary. In part

II, the Green's function is extended to the case of generalized

plane strain and applied to calculate the stress field in a cubic

solid containing a Z 5 boundary and to fiber-reinforced composites

with various fiber orientations subjected to an out-of-plane load.

It is found that, as predicted by earlier authors, the stress

is singular at the intersection of the free surface and the

interface in both plane-strain and generalized plane-strain case.

These singularities in the stress field associated with the

presence of the free surface are identified and discussed. The

singularities in the stress field are found to be of the type

—5

r , ln(r) as well as containing higher powers of ln(r), where 5

is between 0 and 1 and r is the radial distance from the

intersection of the free surface and the interface. In addition,

it is found that, in general, the stress field will also contain

an oscillatory factor of the type exp [ igln (r) ], where g is a

constant which depends upon the material parameters of the two

solids. However, for the case of the generalized plane strain in

fiber-reinforced composites, the factor g is almost 0 so that the

oscillatory behavior is absent.

Part I also gives a simple and convenient method for solving

an inhomogeneous generalized vector Hilbert problem which contains

xi
a nonsingular kernel in addition to the usual singular Hilbert

kernel. The solution of this Hilbert problem is required for

obtaining the Green's function in the present case, as well as

many other cases dealing with the stress analysis of composite

materials containing surfaces, interfaces, and cracks.

The solution of the Hilbert problem is obtained by using a

complex transform of the type y

L<^ 0,5

/ where y and q are variables

on the real axis. It is shown that this function is an

eigenfunction of the Hilbert operator and is orthogonal for y,

varying in the range between 0 and ± » and q between -oo to + ».

xii
PART I

Elastic Green's Function for a Composite Solid Containing

a Free Surface Normal to the Interface

by

V. K. Tewary

1. Introduction

In this part, we describe the calculation of the SI Green's

function. We give detailed results for the displacement and the

stress Green's function and a discussion of their singularities.

For the purpose of illustration, we apply the method to a simple

antiplane-strain problem for a cubic solid containing a £-5 grain

boundary and also to the corresponding plane strain.

In calculating the SI Green's function, we have used the

Green's function for an infinite bimaterial composite as the

starting Green's function. A summary of these functions, as taken

from two earlier papers [18,19], has been given in Appendix I-A.

For calculating the SI Green's function, we need to solve an

inhomogeneous generalized vector Hilbert problem. Our method for

solving this Hilbert problem has been described in Appendix I-B.

The evaluation of certain integrals as required for various

formulae in the paper has been described in Appendix I-C.

2. Calculation of the SI Green's Function

We consider a bimaterial composite with a planar interface

and a free surface normal to the interface as shown in figure 1-1.

This figure also shows the coordinate system used in our

calculations. We take the origin of our coordinate system at the

intersection of the free surface and the interface. The X-axis is

taken along the interface and the Y-axis is taken along the

surface. The Z-axis is taken normal to the plane of paper. We

assume that there is no variation of the field quantities

(stresses, displacement etc.), which enables us to consider a

two-dimensional problem as shown in figure 1-1. Although the field

quantities will remain constant in the Z-direction, they can and,

in general, will have a nonvanishing Z-component.

3
The notation in this paper is generally the same as in I and

II. The solid in the upper half plane (UHP) and its parameters

will be labelled by the superscript A and in the lower half plane

(LHP) by the superscript B. The X-, Y-, and Z- components of a

quantity will be identified by subscripts 1, 2, and 3,

respectively. Following the usual notation in the theory of

elasticity, we shall denote the radius vector of a point by the

vector x with its X-, Y- , and Z- components denoted by x^ x

2

, and

x

3

, respectively. The x

2

coordinate may also be represented by y.

The Roman indices i,j, k, etc. will take the values 1, 2 or 3

and will denote the Cartesian components. Unless otherwise stated,

we shall assume the summation convention over repeated Roman

indices but not over repeated Greek indices a, /3, A etc., which

will also take the values 1, 2 or 3. We shall denote by

Greek i and not by the more usual symbol i which, as mentioned

earlier, has been used for labelling Cartesian components.

We consider a case when a line force 0 is applied to a solid

in the XY- plane at the position vector x' • Then the displacement

field u(x) at point x is given by the following equation for

elastic equilibrium:

a 2

u.

c

ikj l 3_ = 0 5(x _ x <) 6 (x - x

2

) , (1)

Sx^ Sx^

where c denotes the fourth rank elastic constant tensor and 5 (x)

is Dirac's delta function. We assume that 0 is applied in the UHP

(x

2

' > 0) . Equation (1) is written separately for UHP and LHP with

c and u labelled appropriately by A or B for UHP and LHP,

respectively. In the equation for LHP, the RHS will be 0.

We prescribe the following boundary conditions (for i=l-3)

4
u

A

i [x1/ 0] = u

B

.[x1/ 0] (0 ^ x

x

^ co), (2)

x

A

i2 [ Xl/ 0] = T

B

i2 [x1 ,0] (0 s xx s »), (3)

T

A

i:L

[0,x

2

] =0 (0 * x

2 * co) , (4)

and

T

B

il [0 / x2

] =0 (0 * x

2 * -oo) (5)

where t is the stress tensor. The stress components z . . and z. 0

can be obtained from the displacement field by using

and

aukauk

T

i2 = c i21k ~ + c i22k ~ ' (7)

ax

l

dx

2

Equations (2) and (3) are the continuity conditions which

imply that the two solids are perfectly welded at the interface.

Equations (4) and (5) are the free surface conditions at x

1 = 0.

Our object is to solve eq (1) subject to the boundary

conditions given by eqs (2) - (5) . The Green's function will be

obtained from this solution by taking the force £ equal to unity.

The corresponding Green's function for the case when the unit line

force is in LHP can be obtained by replacing x

2

' by and

appropriate relabelling of the solids in the UHP and the LHP.

The Green's function, which is a solution of eq (1) and which

satisfies the boundary conditions given by eqs (2) and (3) , has

5
been obtained in I (see also Appendix A of II) and quoted in

Appendix I -A of this paper. We shall use this result to obtain the

solution of eq (1) , which would also satisfy the boundary

conditions given by eqs (4) and (5)

For this purpose we shall use the standard technique for

Green's function calculations (see, for example, I and II) and

apply a hypothetical distribution of line forces just outside the

region in which the equation has to be solved. We then determine

the force distribution by requiring that the solution satisfies

the prescribed boundary conditions. In addition, we shall impose

the condition that the total space integral of the forces is 0.

A

We denote the force distributions by F (t) (0 * t * ») and

F (t) (-co < t ^ 0) in the UHP and LHP. As shown in figure 1-1,

these forces are applied just outside the free surface at a

A continuous set of points in the UHP and LHP viz. F (t) at x=-e

and x

2

=t (0 < t * oo ) and F (t) at x

1

=-e and x

2 = t (-00 < t * 0) , where e is a small positive number which is taken to be zero in

the limit.

From the definition of the Green's function, we can write the

displacement field in UHP and LHP; i.e., a solution of eq (3), in

the following form:

o

+

—00

(8)

and

(x;t)

A

(t)dt

0

6
(9) —00

where G

Q

denotes the Green's function for an infinite bimaterial

composite which was derived in I. In eqs (8) and (9), we have not

explicitly written -e , the x.^ component of the position vector of

the second variable in the arguments of the Green's functions

under the integral signs, for reasons of notational brevity.

A

It may be emphasized that the so-called force functions F (t)

and F (t) are arbitrary functions which have to be determined so

that the boundary conditions given by eqs (2) — (5) are satisfied.

These functions have the dimensions of force but are not physical

forces applied to the solid.

Using the result given in eqs (A.1)-(A.4) of Appendix I -A, we

obtain the following expressions for the displacement field from

eqs (8) and (9)

a

a/3

0

00

JT

1

0

Re zA

+ c-p^ t) F

A

(t) dt

7
„0

+ i

7T Re

—oo ZA, . ^III , A B . , 1 («) ln(z

a

e ~p/3 t}

a/3

F

B

(t) dt

and

(for 0 s x

2

s ») (10)

u (x) =

1 n V" -B. . _II , ,-B -A ,

,

- ^

Re

2^ Z

(a)

2/3 ln(z

a~ p0

^

a/3

00 1

TT Re Y I

B

(a) Q

-a/3

II i ,-B, -A .

_ ln(z + e-p0 t)

/3

v a ^/3 ' F (t) dt

.0

1

TT Re ^

J

B

(a) ln(z B

+ e-p

B

t)

—00 a

,B F (t) dt

.0

+ i

TT Re V" -B, . .IV , ,-B^ B

> 1 ( a ) ln ( z ~+ c_Po t)

/3

—00

La/3

F (t) dt

(for -oo < x

2 * 0) , (11)

where Re denotes the real part, overhead bar denotes complex

conjugate and other symbols have been defined in Appendix I-A.

Equations (10) and (11) will satisfy eq (1) everywhere in the

region of solution, since the line of application of the forces is

just outside this region, which is ensured by the presence of c.

In the end we shall take the limit c = 0.

8
Further, since the Green's function used in writing eqs (10)

and (11) satisfies the continuity conditions at the interface as

given by eqs (2) and (3) , the solution as given above will also

A B satisfy these conditions for all values of F (t) and F (t) . Our

main task, therefore, is to determine these two force functions so

that the free-surface boundary conditions as given by eqs (4) and

(5) are also satisfied. The resulting solution will then satisfy

all the prescribed boundary conditions.

In addition, we impose the following condition on the force

functions to ensure that the displacements are single valued [23]

everywhere except possibly at the cut:

F

A

(t) dt + F

B

(t) dt = 0 . (12)

0 -co

As described in I and II, a simple method for writing the

stress component x.- from eqs (8) and (9) is to differentiate the

log functions in these equations with respect to their arguments

A B A B and replace the matrix y ' by the matrix Q '

, where Q is defined

below (for superscripts A or B)

fi.^a) = L

1 ^ (a) rkj (a) (13)

and

L

1

.. (a) = c. 1iv + p„ c. 10V . (14)

The stress components t^

2

can similarly be obtained by

replacing y by a where a has been defined in Appendix I-A. In what

follows we shall denote the three components t ^ by the vector T^.

Similarly, T

2

will denote the components t^

2

-

9
Thus, we obtain the following expressions for

(8) and (9) . from eqs

a

a/3

+ L I ,V, h-(!^L] + 5V)„M^L]

«e L pp

1 p£ J

and

- L I

a/3 l P/3

-A

if (15)

3Tj(x) = -

a/3

+ L

r -B B

a0 L p~

+ L a Pa

-B , z + c

a

10
- L I

a/3

K

BB (a,/3) H

-B

Bf z a )

v P

B

/3

+ KBB (a,/3) H

B

B

P/3

(16)

where

3

A

< Z > = WZ

"M

F

A

(t)

t - z dt , 0

J

(17)

5

B

( Z

) = 2WT

-° F

B

(t)

t - z dt , —00

(18)

K^ia.P) = QA

(a) / p

A

, (19)

KAB ( a/ ^) = Q

A

(oO Qp

11

/ p* # (20)

KBA (a^) = fi B

(a) Q* 1

/ p

A

, (21)

KBB (a^) = n

B

(a) Q*V

/ p

B

, /3 (22)

,A,B, . oA,B, . / A,

J ' (a) = Q («)/Pa

(23)

and z, in eqs (17) and (18) , is any complex variable. We have

neglected e in the arguments of the H-functions in eqs (15) and

(16) except in the third term on the RHS, because these are the

only terms which will be sensitive to e in the limit c=0.

A B We now equate u.^x) and z. , (x) to zero at x = 0. At x.. = 0,

A B the variables z and z assume the following values:

11
A.B A.B

z = p y , (24)

where we have denoted the x^- component of x by the variable y

Using eq (24) in eqs (15) and (16) , we obtain the equations

(i) For 0 ± y 2= w

TAHA(y-ie) + F

A

HA(y+ie)

a(3 >- K^Ca,^) H^a^y) + KAa (a,/3) H

A

(5/vft y) a/3-1 ' -AA a/3-

"I

af3

K

AB (a,P) H

B

(b

og

y) + KAB (a,p) H

B (b^y)

= -^T Re Z ^<a > ^ Pa/Pa]"' i

a

a0

(25)

and

(ii) For -oo < y < 0

I

a£ l

KBA («,0> HA (d^y) +KRa(a/ /3) HA (d^y)

a/3-7 ' -BA a/3-

+ T

B HB (y-Le) + TB HB

(y+Le)

12
-z

af3

K

BB («,0) H

B

(e

aey) + K

BB ( a ,(3) H

B

(e

agy) 'a/3- BB

- 1

a0

(26)

where

r

A ' B = £

JA ' B(a) , a

A,

(27)

P(a,/3) = JA

(«) Qj£ » (28)

R(a,/3) = J

B

(a) Q* 1

, (29)

A , -A

aa0

= P<x / P/3

' (30)

b

«0

= Pa / P

^

(31)

B = Pa / P13 (32)

-B B

e _ = p / p_

a/3 a ' ^/3 (33)

In eqs (25) and (26) , we have used the fact that, in view of

eq (A. 20), the sign of the imaginary part of the terms of the type

e/p^ (for superscripts A or B) will be negative. In the limit e =

0, the functions H(y ± e /Pa

) will be sensitive only to the sign of

the imaginary part and not its magnitude, which approaches 0 any

way. Accordingly, in the argument of the H-functions, we have

written e/p' as -lc and e/ p as +lc. *a ' *a

It may be remarked that the quantities J, K, T, P, etc. (with

various subscripts and superscripts) are 3x3 square matrices,

13
whereas H, F, and <£ are 3x1 column matrices or vectors.

Equations (25) and (26) together represent a 6 x 6 matrix equation

or a set of 6 simultaneous equations in 6 unknowns. The unknowns

A B

are the three components each of F (t) and F (t) . These equations

constitute a generalized inhomogeneous vector Hilbert problem. Its

solution has been described in Appendix I-B.

Following the method described in Appendix I-B, we write the

two force functions in terms of their complex transforms as given

by eqs (B.25) and (B.26). In Appendix I-B, the functions F(t) and

H(z) are assumed to be vectors having two components. In the

present case, these are vectors having 6 components. However, the

formulae given in Appendix I-B are ABA formally valid in the present

B

case as well, provided we regard F (t) , F (t) , H (z) , and H (z)

themselves as 3 component vectors.

Thus, as in eqs (B.27) and (B.28), we write the following

representation for the H-functions:

TA

r ( z ) V^qjz^"0,5 dq (34)

—00

and

00 H

B

(z) = TTB / x tq-0.5 , ,. c . V (q)z M dq , (35)

where (q) and V (q) are the unknown functions of the real

variable variable q, which have to be determined. For each q,

these are 3x1 column matrices like the H-functions.

Once the H functions are known, the force functions can be

obtained in a simple manner in terms of the discontinuity in the H

function over the real axis by using the Plemelj relation (see,

for example, [20]) as given by eq (B.42).

14
We now substitute for the H functions from eqs (34) and (35)

in eqs (25) and (26) and use the method described in Appendix I-B

for deriving eqs (B.31) and (B.32). This would transform each

element of the matrices in eqs (25) and (26) into a function of q

alone. As remarked earlier, we can represent eqs (25) and (26) by

a single 6x6 matrix equation. Accordingly, we regard V(q) as a

single 6x1 column matrix or a six-dimensional vector. Thus, eqs

(25) and (26) , after the transformation can be written in the

following matrix representation:

M(q) Y(q) = N(q)/E(q)

,

(36)

where M(q) and N(q) are respectively the 6x6 and 6x1 matrices

formed by the transformed elements of the LHS and RHS of eqs (25)

and (26)

.

In order to represent the matrices M(q) , N(q) and V(q) , it

would be convenient to introduce a block matrix notation in terms

of A and B which correspond to the solids in UHP and LHP

respectively. According to the block scheme, we first write the

6x1 matrix V(q) in terms of two 3x1 blocks V^q) and V

B

(q) and

then write M(q) and N(q) in terms of the same block

representation. Their matrix elements will thus be denoted in

terms of the block matrices M

AB , , N^, etc. For notational

brevity, we shall not show their functional dependence on q except

where explicitly needed.

The block structure of eq (36) can thus be written as

M

-AA

M ^ -BA

M

-AB

M

-BB

1

> Y

B

(q),

E(q)

-A

-B

(37)

where

15
= -r

A

exp(-27rq) + fA

+ Zccfi

(38)

a/3>-

a/3 (39)

M

-BA = Z

a/3 L

(40)

M

BB

= Lr

Bexp(Trq)-LfB

exp(-7rq)

-z

a/3

(41)

NA exp (-2Trq)

2n

a

JA(a)

' Ay h"

Pa/Pa

u -A, % f -A /-Al

+ J"(a) /

exp(-2Trq) V"

271

a/3

P (a,/3) r-A/ a) u VP

«J

+ P(a,/3)

r A/-ai^

<t , (42)

N

B _ l exp (-Tig)

27T E(q) z

a/3

R (a,/3) (-A ,-B)

V

1 J + R(a,/3) i / (43)

and

E(q) = 1 + exp(-2Trq),

l± = tq - 0.5 .

(44)

(45)

16
Using the same representation as given in eq (37) , the block

structure of the matrix H(z) is

H(z) =

f A

H (Z)

I H

B

(z)J

(46)

Equation (37) can be solved for V(q) by matrix inversion:

-1

V(q) = M N / E(q) . (47)

For later use, we shall also write eq (47) in block representation

corresponding to eq (36) as

I Y

B

(q)J

E(q)D(q)

c 1 -BA

Sab

Sbb

> ' 1 Sa

J

(48)

where C are the block elements of the matrix of cofactors of M in

the same block representation and D(q) is the determinant of M.

The matrix of cofactors is proportional to the inverse of M

which, in the present block representation, can be formally

expressed as

1 StofcrtW*) = 1 ?Wi>saB <q> " D(q) 5 ab - ( 49 > a a

where 6,_, is Kronecker's delta and A, a. and B are block indices AB '

' which label the blocks in the corresponding matrix representation.

As in the case of M and N matrices, we shall show the functional

dependence of the elements of C on q only where explicitly needed.

27
In general the inverse of M in eq (48) has to be obtained

numerically, and there is no need to calculate the cofactor

matrix. The two H functions can then be calculated from eqs (34)

and (35) by integrating over q. The integrals over q in these

equations can be calculated numerically in a straightforward

manner. It is also possible to evaluate the q integral

analytically by using contour integration. For this purpose eq

(49) is more convenient. This would give a series representation

of the H-functions, which will be derived in the next section.

Equation (47) gives the particular solution of Hilbert's

equation. To obtain the general solution of Hilbert's equation, as

described in Appendix I-B, we need to add the following

homogeneous solution

_ • H

Q

(Z) = Y h(Q

r

) z LQr"°* 5, (50)

r

where Q is a root of the determinantal equation

D(Q

r

) = 0, (51)

and h(Q

r

) is an eigenvector of M(q) corresponding to zero

eigenvalue for q = Qr

; it is a solution of the matrix equation:

M(Q

r )h(Qr

) = 0 . (52)

In analogy with vectors H(z) and V(q) the two block

components of h(Q

r

) will be labelled by the superscripts A and B.

According to the theory of homogeneous equations, one of the

components of h will be arbitrary. This provides a set of

arbitrary constants (one for each value of Q

r

) , which have to be

chosen so that eq (12) is satisfied. These constants can also be

chosen to satisfy any other boundary condition on the stress field

18
such as those which may be required at the outer surfaces in a

finite solid.

Thus, finally, we can write the H- functions which give the

complete solution of eqs (25) and (26) as follows:

H(z)

HA(z)

H

B

(z)

J

(53)

where the subscript p denotes the particular solution which is

given by eqs (34) , (35) and (47) . The homogeneous solution,

labelled with the subscript 0, is given by eq (50)

.

Henceforth, for notational brevity, we shall omit the

subscript p from the particular solution except where we need to

distinguish it from the homogeneous solution. Accordingly, H(z)

will refer to the particular solution for the H-function.

The expression for the H-functions as given above is the main

result of this paper. The stress distribution in the solid is

given by eqs (15) and (16) in terms of the two H-functions. The

displacement field is given by eqs (10) and (11) in terms of the

logarithmic integrals which can be obtained from the H- functions

by a simple integration over z. Equations (10) and (11) give the

required Green's function when £ is taken to be a unit force.

In the next section, we shall derive a series representation

for the H-functions and discuss their singularities.

3. Series Representation of the H-Function

A

In the previous section we obtained the two H functions H

and H , as integrals over q. In this section we shall obtain

19
A B series representations for H and H in the Cartesian space and

discuss their singularities. The singularities in the H-functions

give the singularities in the stress. We shall also obtain the

A B

space integral of the force functions F and F in order to ensure

that eq (12) is satisfied.

A B First we evaluate H (z) and H (z) from eqs (34) and (35) for

A

a complex variable z. Using eq (48), we write H (z) in the form

HA(z) =

-co'

" m

T-.A, . iq-0.5 , £ <q> 2 dq

(54)

E(q) D(q)

where

?

A(q) = [^AA (q) ^A (q) + ^AB (q) ^B (q)

]

(55)

g with a similar expression for H (z)

The q-integral in eq (54) can be evaluated by using the

method of contour integration as described in Appendix I-C for eq

(C.16). From the result given by eq (C.24) we obtain the following

A B expressions for H (z) and H (z) (in the limit e = 0)

HA(z) =Y^Y}*A{n,e) 2 (n+1) "C ^(«,n)

oc n

(r,e) z lQ -0.5+e

r

20
a£ L n

,-(11+1) -e

£ (af, #n)

r—^ iQ -0.5+e

+ ) ?AA (r ' c) Z

r

£2

(«/3/LQr

)

a/3 L n

,-(n+l) -e

£3

(a/3,n)

^ lQ -0.5+e

+ 2^ ^AB (r ' C) Z

r

g3 (aP/t-Qr

) r

s

z + Vz > (56)

and

H

B

(z) = ]T ^BA(n ' e) z " (n+1) C ^(^n)

L n

-, lQ -0.5+e

a/3 L n

,-(n+l)-e

£2

(a/3,n)

21
^—i lQ -0.5+c

+ 2^ ^BA (r ' C) 2

r

g2

(a/3,(.Q

r

)

a/3 L

-(n+l)-c

£,(*P»n)

LQ -0.5+C

2^

*BB (r,e) z r

g3 (a/3,LQr

) s

z + RH

(z) f (57)

where Re denotes the real part and other symbols are as defined

below:

-AA (n ' C) =

- LexP("2TrLC ) 9aa (qn+LC) /D(qTi 4 +2tc) , n

n

-AB (n/C) = L(_1) exP("TrLe )2 aR (^+l- e )

/

D (q„+2Le ) / AB v ^n l n

$ BA (n,e) = -Lexp(-27TLC) CRa (q^+Lc) /D(qri +2(,c) , BA ^n l n

n

$ BB (n,c) = l (-1) exp(-TTLc) gBB (qn+Lc ) /D (qn

+2 lc ) ,

(58)

(59)

(60)

(61)

2TTLexp[-2Tr(Q

r

-LC) ] 9AA (Qr-LC)

E(Q

r

-2LC) D' (Qr

)

(62)

*AB (r,c) = -

27T eXp[-7T(Q

r

-LC) ] CAB (Qr-LC)

E(Q

r

- 2 L c ) D' (Qr

)

(63)

22
27TL exp[-27T(Q

r

-LC) ] 9BA (Qr-Le)

E(Q

r

-2LC) D' (Qr

)

(64)

*BB (r,e) =

2n exp[-Tr(Q

r

-Le) ] 2gg (Qr

~ LC )

E(Q

r

-2te) D' (Qr

)

^(oc f n) = 2 Re JA(a) [p

A

/ p

A

]

n+e

I^OMQ,.) = J*(«) [p* / p*]

-LQ

r

~0 . 5-e

-lQ -0.5-e

r + SN«> [p* / i£]

£ 2

(<x/9,n) = 2 Re P(a,/3) [p

A

/ p

A

J

n+e

(65)

(66)

(67)

(68)

£ 2 (af3,LQr

) = P(a,/3) [p

A

/p

A

]

-lQ -0.5-e

r

+!(«,(»

[

P

* / p*]

-lQ -0.5-e

r

(69)

£ 3

(«0,n) = 2 Re R(a,|3) [p

A

/ p*

]

n+e

(70)

Z 3 (otfl,LQr

) = R(a,/3) [p

A

/p*]

-lQ -0.5-e

r

23
(ot,0) [p

A

/ p*]

-lQ -0.5-c

r

+ R (71)

The vectors £ 2

, £ 3 , £ lf and £

2

are in units of £/27T. The

factor s

z is +1 if the integration contour was chosen to be in the

UHP (that is, if the expansion parameter is greater than 1) and -1

otherwise. In deriving the above equations, we have taken the

following values of the branches of the complex exponential

As described in Appendix I-C [see eqs (C.22) and (C.23)], the

allowed values of n and Qrin each sum depend upon the magnitude

of the expansion variable. For example in the first term on the

RHS of eq (61) , n will take 0 or any positive integral value if

A A

mod (z p / pa

) is larger than unity and negative integral values

otherwise. Thus n and also of course Qr depend upon the indices of

the outer sum. The factor s in eqs (56) and (57) , which defines

z the sign of each expansion term, also depends upon the indices of

the outer sum. Because of this dependence we cannot write the

above equations in a compact form for a general value of z.

However, in the two limiting cases when mod (z) -> 0 or mod

(z) -> oo, n as well as Qr

become independent of a or jS. In such

cases the above equations are considerably simplified. Since these

cases are of substantial interest, we give below the matrix

representations for the two H functions in the AB block

representation in the two limiting cases:

exp(TTL/2) = i and exp(-7ri,/2) = -l . (72)

H(z) =±£ $(n,e) £(n,e)z -(n+1) -e

n

lQ -0.5+e

± BH ( Z )

(73)

r

24
where $ and * are square matrices, and H, R^, £ and £ are column

vectors with elements labelled by A and B in the same block

representation as defined by eq (37). As in eq (C.24), the plus

sign is to be taken if mod(z) > 1 and negative if mod(z) < 1. The

matrix elements of $ and ¥ have already been defined in eqs (58) — (65). The elements of £ and £ are defined as

^A (n) = I Ma ' n) +

I £2 («0' n )L (74)

£ B

(n) = £ £3(«/3,n) , (75)

a/3

^A (LQr

} = 1 £1 («' L Qr > + 1 22

(a,3 ' LQ

r

)

' (76)

a a/3

^B (LQ

r

) = I M0^'^^ ' (77)

aft

The dependence of the various parameters in eqs (74) -(77) on

e is not shown for reasons of notational brevity. As described in

Appendix I-C, e can be set to zero immediately for those values of

n and r for which there are no higher order poles. When the higher

order poles are present, their contribution has to be calculated

separately by evaluating the limit at e = 0. In eq (73), this

contribution has been denoted by R„(z) . Equations (73) , (56) and

( 57 ) , therefore, provide a convenient representation for the H

functions which is formally valid whether higher order poles are

present or not. In section 5 we shall consider the

antiplane-strain case in which all qn

's are also roots of eq (51)

and are therefore second-order poles.

As explained in Appendix I-C and given by eq (C.22) and

(C.23), the position and the contribution of poles depend upon

whether the contour of integration has been chosen in the UHP or

LHP. This, in turn, depends upon the magnitude of z or the

25
expansion variable. In sections 3B and 3C, respectively, we have

calculated the contribution of the poles in the two limiting cases

when mod(z) -> » and mod(z) -» 0.

We see from eq (73) , and also from eqs (56) and (57) , that

the H-functions contain basically two series representations, one

in qn

and the other in Q , and a remainder term which is

nonvanishing when there are higher order poles in the integrand.

The terms containing qnare actually those which give the

particular solution of eqs (25) and (26) . The terms containing Qr

are the solution of the homogeneous equation because, when

substituted in eqs (25) and (26) , they lead to the homogeneous

matrix equation (36) (unless qn

= Qr

for some n and r)

Thus, although eqs (56) and (57) are particular solutions of

the Hilbert's equation in the q representation, they may still

contain a homogeneous part which separates out when integrated

over q. We shall, however, refer to these terms also as part of

the particular solution to distinguish it from the homogeneous

part which has to be separately introduced in order to satisfy eq

(12) . This part will be discussed a little later in section 3D.

3A. Behavior of H^(z) for large values of mod(z)

We shall now discuss the behavior of H(z) given by eq (73)

for large mod(z) . In this limit, as given by eq (C.22), only the

poles in the UHP will contribute to the integral. The

corresponding allowed values of n and Qrare then given by eq

(C.22). For a moment if we ignore any higher order poles, we find

from eqs (73) and (C.22), that H(z) will consist of only negative

powers of z. This shows that the stress will vanish as mod(z) -> «.

The first series in eq (73) contains negative integral

powers, whereas the second series, in general, will contain

26
negative fractional powers of z. Let us write Qr

in terms of real

and imaginary parts as

Q = g + t.k , *r ^r r' (78)

with the constraint

k

r = Im(Qr

) > 0. (79)

The second series in eq (74) can then be written in the form

The subscripts p2 on the vectors V and H identify them as

belonging to the second series of the particular solution. In eqs

(80) and (81) , we have put c = 0 consistent with our hypothesis

that the particular value of Q is not a higher order pole.

We notice from eq (81) that each term falls off with

increasing mod(z) . The exponential factor for each term will be

nonunity whenever gr

, the real part of Q corresponding to that

term is nonzero.

This discussion is valid for those terms which arise from

only simple poles. However, as shown below, at least one pole,

that at n=0, is a second-order pole. For this pole, therefore, we

have to evaluate the proper limit at c = 0.

It can be verified that D(q) is zero at q = q where q has

(80)

[for mod(z) -» oo]

,

where

Y 2 (Qr

,c) = *(Q

r ) £(LQr

). (81)

27
been defined in eq (C.25). For this purpose, we note from eqs

(38)-(41) / (84) and (85), and by using the definition of various

parameters as given in section 2 and Appendix I-A, that at q = qQ

, the matrix elements of M(q) and £(q) reduce to

*AA<V " 2 Re

[

^ + C

A

2s ]

• (82)

Sm<V =

"

2 Re

[ ? 25

U

]

•

< 83 >

5ba<V - " 2 Re

[ C

B

2"

]

- (84)

?WV - " 2 Re " C

B

2^ ]

. (85)

£a<°>= ' (86 >

£B

(0) " " ^ba'^q) • ( 8? )

The quantities Q

1 - Q

IV in eqs (82) -(85) have been defined in eqs

(A. 7) -(A. 10) in Appendix I-A . They should not be confused with Qr

which denotes a particular value of q.

Using eq (A. 25) of II, we find that

- ?WV = ?WV - !?bb<V - °< < 88 >

which shows that

D(q

Q

) = 0. (89)

Since E(q) also has a zero at qQ[see eq (C.17)], we see that the

integrand in eq (54) has a second-order pole at q = Q. = i/2.

28
The contribution of this second order pole has been evaluated

in Appendix I-C. Since this pole is in UHP, it does not contribute

to the integral in the limit mod(z) -» 0. For large values of

mod(z) , we calculate its contribution by using eq (C.29) as

described below.

By comparing eqs (C.16) and (54), we note that in the

calculation of (z) , the function P(q) as defined by eq (55) for

q = qQ

is given by

E

A

C3o» - [£aa<«o> SA (q0 ) +

£AB <V HB <V]- (90>

In order to achieve some simplification, it would be useful

at this point to derive certain relations between the block

elements of the matrix C(q) at q = q Q

, From eqs (49) and

(86)-(87), we see that

Wo' £a<°> " sAB C30 ) £B

(°> = D <V - 0 < 91 >

and

SbA^0» £a< 0 ' " SBB (q 0 ) CB

(0) = 0. (92)

We can derive the following relations by taking the product

of M(q) and C(q) in that order and by using eqs (49) and (88)

Saa^o' + SBA <V = 0 < 93 > and

SAB (q0

) + £BB (q0) = 0. (94)

Now we come back to eq (90) . From the definition of N(q) as

given by eqs (42) and (4 3) , we note that

W =

" £a (0) (95)

29
and

SB

(q0) =

^B (0) * (96 >

From eqs (91) , (95) and (96) we find that eq (90) gives

PA

(q Q

) = 0. (97)

Similarly, we can show that

P

B

(q Q

) = 0. (98)

Hence we need to keep only the first term on the RHS of eq (C.29).

A similar argument applies to the calculation of H (z) . Thus, in

the limit c = 0, we get the result [for mod (z) -» «] :

R^(Z) =

L Z

1

P

A/

(q Q

) (99)

and

D' (q Q

)

l z"

1

..B' R„(z) = — P (q n

) , (100)

where the primes denote differentiation with respect to the

argument

.

Equations (99) and (100) give the contribution of the

second-order pole at qQto the H functions, if there are no other

higher order poles in the integrand of the H-functions. In case

some other values of q and Q coincide, we need to take the limit

c = 0 of those terms in a similar manner, and they would be

included in R (z) . In those terms in eqs (56) and (57) , in which

qn

is not equal to Qr

, c can be put equal to zero immediately.

The pole at q^ as discussed above is effectively a first

order pole because P(qn ) is also zero along with D(qn

) at q = qn

.

30
This is the reason why logarithmic terms cancel out in the final

result given by eqs (99) and (100)

.

In deriving eqs (99) and (100) , we have only considered the

degeneracy of the poles between set I and set II that is between

qn

and Qr

as defined by eqs (C.17) and (C.18). It is obvious from

eq (C.17) that there cannot be a degeneracy among the poles of set

I because no two q 's can be equal.

. n

The poles of set II as defined by eq (C.18) might be

degenerate among themselves. Since the matrix M(q) is 6x6, it

would have six eigenvalues. The determinantal equation (C.18) can,

therefore, have up to 6-fold degenerate roots. If tj is the

degeneracy of a particular root, that is if tj number of Q are

equal, then, as given by the Cauchy integral theorem, their

contribution to the integral * will contain a factor of the type

77-I [ln(z)] . The function H(z) will not diverge at 00 because the

logarithmic term will be multiplied by a negative power of z.

3B. Behavior of H (z) for low values of mod(z)

P

In this subsection we shall discuss the form of the

particular solution for H(z) in the limit mod(z) -> 0. The behavior

of H(z) for low values of mod(z) can be analyzed following the

same steps as in the preceding subsection for large mod(z). The

singularities near mod(z) ~ 0 in the total H(z) , including the

homogeneous solution, will be discussed separately in section 3E.

The expression given in eq (73) for H(z) is also valid in the

limit mod(z) -> 0, but the allowed values of n and Qr

will be those

given by eq (C.23). The first series in eq (73) will contain only

positive integral powers of z. These terms will vanish in the

limit mod(z) -> 0. The second series may contain also fractional

powers. Except for those values of Q which have an imaginary part

31
greater than -0.5, the powers of z in the second series will also

be positive. This shows that, apart from the exception given

above, the terms in both the series in eq (73) will approach 0 as

mod(z) goes to 0.

In the second series, we again write Q in terms of its real

and imaginary parts as given below [c.f. eq (78)]

Q = g - Lk (101) *r ^rp rp v '

with the constraint

k = -Im(Q ) > 0. (102) rp v r' v '

The subscript p on g and k

r

is introduced to identify this term

as coming from the particular solution. The resulting expression

for H ~(z) , which is analogous to that given by eq (80) is

P<£

(k -0.5)

H

p2

(z) = 1 Yp2 (Qr

) [z]

^ exp[Lg

rp

ln(z)] (103)

for mod(z) -> 0.

Now we consider the contribution of the higher order poles

following the same procedure as used in the preceding subsection.

It may be verified that D(q) is zero at q = Q = -i/2. Since E(q)

is also zero at q = q_1 = -i/2, the integrand has at least one

second-order pole at q = Q = -l/2. Unlike the case for mod(z)

-> oo as discussed in the previous section, P(q_

1

) is not zero in

this case. The contribution of this second-order pole which is

represented by the remainder term, has been evaluated in Appendix

I-C [see eq (C.30)] and is given below:

32
D' (q. x

) L

P'(q_

1

) + LP(q_

1

) Ln(z)

-0.5 P(q

-:L

)D" (q_ 1 )/D' (q_ 1

) + 7TP(q_

1

) (104)

Equation (104) shows the ln(z) behavior of H(z) at low values

of mod(z) . In deriving this equation, we have assumed that the

root Q_1of the determinantal equation (C.18) is not degenerate.

As remarked earlier, if eq (C.18) has a 7}-fold degenerate root at

Q

r7J

/ H(z) will have terms of the type,

R

R

(z) * [z]

(kv -0.5)

r7] ' T) TM 1 [ln(z)] 7? ~ exp[Lg

r7)

ln(z) ] , (105)

where gr^ and k

r^ are, respectively, the real and imaginary parts

of Q . If k is half-integer and g =0, then E(q) will also be

0 at q = Q . This is similar to the case previously considered

for which Qri^ = -l/2. This will lead to an additional power of

ln(z) as in eq (104)

.

3C. Contribution of the homogeneous solution

We now consider whether the homogeneous solution H

Q

(z) has to

be included in the total solution or not. For this purpose we need

to calculate the integral of the force functions with the

objective of ensuring that eq (12) is satisfied. The integral

required for calculating the force function from the H function

has been evaluated in Appendix I-C. However, we shall not need the

explicit expressions for the force functions for calculating the

Green's functions; a knowledge of the H-functions will suffice. In

this section, therefore, we need only to consider the space

integral of the force function, which also has been evaluated in

Appendix I-C by using contour integration.

33
As shown in Appendix I-C, the only poles which contribute to

the force integral are those at Qr

, which are the zeroes of D(q)

and which satisfy the constraint given by eq (C.35). Each value of

Qr corresponds to a force term containing LQ

r

~0.5 as the power of

y and a stress term containing the same power of y. Each of these

terms will be the homogeneous solution because, by definition,

D (Qr

) will be 0 for all of them. We, therefore, choose additional

force terms such that their integral is negative of that given by

eq (C.34) so that the total integral of the force function is 0 as

required by eq (12) . The corresponding H-functions which denote

the homogeneous solution are written in the form,

h£(z) = - £ v£ (Qr

) z LQr"

°- 5

, (106)

Qr

*<-/2

5o< z

> = " I

v

o <V zLQr lQ - 0.5

Qr

*L/2

and

(107)

exp (ttQ ) P

B

(Qr ) V°(Q ) = 27T — , (109)

E(Q

r

) D' (Qr

)

with the constraint

0 ^ Im(Q

r

) * 0.5 . (110)

3D. Singularities in H(z) at mod(z) = 0

Now we discuss the singularities in H(z) as mod(z) approaches

34
0. As given by eq (110), H^z) contains two main terms, one is the

particular solution which is an integral over q, and the other is

the homogeneous solution which is a sum over Q , where Qr

is a

root of eq (51) with the constraint given by eq (110)

.

First, we consider the homogeneous solution. As described in

the preceding section, this term will arise only if eq (51) has a

solution QrQ which satisfies the constraint given by eq (109) . The

additional subscript 0 has been introduced to identify this term

as arising from the homogeneous solution. This singularity can be

expressed as

-(k +0.5)

H(z) * (z)

ru exp[(,gr0 ln(z) ] , (111)

where the weight (coefficient) of each singular term is given by

eqs (108) and (109) , and grQ and k

rQ are, respectively, the real

and imaginary parts of Q n ,

Qr0 = gr0 + «.kr0 , (ii2)

where

0 s k

rQ

s 0.5 . (113)

We see from eq (111) that the exponent of the singularity

varies between 0 and -1. In addition, if the real part of QrQ (in

addition to the constraint on its imaginary part) is nonzero, H(z)

will also have a logarithmic—cos or sin[g ln(z)] —type

oscillatory component, which is represented by the exponential

factor in eq (111) . This oscillatory component is of the same type

which arises in the stress and the displacement field near an

interfacial crack in bimaterial composites (see, for example, II)

.

Depending upon the material parameters of the two constituents of

the composite, this term can be significant. The singular term as

35
given by eq (111) will of course be zero if a Qr

does not exist

which satisfies the constraint given by eq (113)

Next, let us identify the singular terms in the particular

solution for H(z) , which is given by an integral over q in eq

(54) . As shown in section 3A, for mod (z) « 0, this integral can

be represented in terms of two series and a remainder term as

given in eq (73) . As explained in section 3B, the first series in

eq (73) has no singularity. The second term in eq (73) is a series

over Qr

, where Qr

are solutions of eq (51) with the constraint

that Im (Qr

) is negative since these terms arise from the poles in

the LHP. This series is given by eq (103) .

We see from eq (103) that only those terms will be singular

at mod(z) = 0 for which k < 0.5. If the real part of Q is

nonzero, the logarithmic oscillatory factor will also be present.

Thus we obtain the following condition for a term in the second

series in eq (73) to be singular:

0 ^ k

rp

s 0.5 . (114)

The corresponding singular terms in H(z) are of the type,

(k -0.5)

H

p2

(z) * (z) rp exp[Lg

rp

ln(z) ] , (115)

where k has to satisfy the constraint given by eq (114) , and the

additional subscript 2 in H identifies the second series in eq

(73) as the origin of this singularity. The exponent of this

singularity varies between 0 and -0.5.

Finally, we consider the singularities in the remainder term

R„(z) in eq (73), which arise if there are higher order poles. The

hi contribution of one second-order pole at q = -i/2, at which both

E(q) and D(q) are 0, is given by eq (104) . This gives a

36
singularity in H(z) which varies as ln(z). This singular term has

(51) has degenerate roots, H(z) will have higher order logarithmic

singularities as given by eq (105) . In general, the logarithmic

oscillatory factor will also be present as shown in eq (105)

.

The singularities in the H-function arise from the

homogeneous part of eq (36) and are, in general, independent of

the inhomogeneity , that is the RHS in eq (36)

.

4. Expressions for the Green's Functions

In this section we shall give the final expressions for the

displacement Green's function and the stress-2 Green's function.

The former gives the displacement and the latter the stress

components t^

2

for a unit line force. The stress-1 Green's

function, which gives the stress components t1for a unit line

force is given by eqs (15) and (16) . We shall give the results

only for the case when the unit line force is in the UHP at the

point (x',y') with p defined as

The corresponding expressions for the case when the unit line

force is in the LHP can be obtained by a simple change of

variables, y' and p.

4A. Displacement Green's Function

no oscillatory factor since the real part of Q is zero. If eq

Pa = x' + p y' . (116)

a

a/3

37
+ L a u I

A

(a)

A

,

-A

f Z + C N -«

a a

+ L

-A

z

- L

-A

z

a/3 L p/3 P/3

(117)

(for 0 s x

2

s oo)

and

G

B

(x,x') = - | Re ^ i

B

(a) Q* 1 ln(z*- p

A

)g

a(3

+ L

a/3

-B

z EBA (a,/3) l^f + ERa (a,/3)

1 A

p

e

; BA

B

Z ^ a

+ L a

B , . TTB

1 (a) y

B^

f z + c

a

B

+ r

B

(«) y

B

-B, V C

-B

- L I

a/3 <-

3

-B

f Z X ERR (a,/3) U" + ERR (a,0) U

B

B

Z

a

BB (118)

(for -oo s x

2

s 0) , where

(119)

38
EAB (a,/3) = y

A

(a) Q

1

/3 (120)

EBA (<*,/3) = i

B

(a) Qj

1

, (121)

(122)

^(z) = -

27TL

00

F (t) ln(z - t)dt

0

J

dz (123)

and

y

B

(z) =

2TTL

F

B

(t) ln(z - t)dt

-oo

J

H

B (z)dz . (124)

In deriving eqs (117) and (118) from eqs (10) and (11) , we

have neglected terms which contain factors of the type ln(p

a

) if

they represent only rigid body displacements. The integrals in eqs

(123) and (124) are indefinite integrals. These integrals can be

evaluated by using either the integral or the series

representation for the H-functions which have been derived in the

previous sections. It can be verified by using eq (6) that eqs

(117) and (118) yield eqs (15) and (16) for stress components t^^-

4B Stress-2 Green's Function

The stress-2 Green's function which can be derived from eqs

(117) and (118) by following the procedure given in section 2 or I

39
and II, are given below:

T

2 (x,x<) = - i Re £ a

A

«x) [z

A

- p

A

]

-1

a

1 _ A /N .I [A -Al - - Re ^ a (a) \z

a

- PfB j

a/3

-1

+ L

(X

, z a

+ c ^

a

-A,

z + e n + f-A, . , -A\ UA( *aT fc i {2 («>/ Pa } «

I -=s-j

a

+ t

-A

z £ Baa^.P) hA

(

-£") + 5AA («^) hA

(

a/3 L P^

- L

a/3 L

-A

DAB (a,^) H

B a

^ P

B

0

+ DAB (a,^) H

B a

-B (125)

and

a0

-1

+ L

a/3 L

-B

DBA (a,^) a

-A

+ (a , /3 ) H'

B

^a

Po '

+ L

B^

( Z + C X a

B

a

+ {?(«)/ P«} H

B

-B,

-B

40
- L I

Ot/3 L

DBB («,0) H

-B

Bf Z

_a_ \

B

P

0

+ DBB (a,^) H

B

B

z a

(126)

where

D^fof, 0) = a

A

(a) / p

-A

/3 (127)

DAB (a,0) = a

A

(ot) Ql

11

s / 0 <3 (128)

DBA («^) -B . . _II / -A

e («) / (129)

and

DBB (<x,/3) = -B. . ^IV / B

2

(a)

2/3 / P/3 (130)

In the next two sections, we shall illustrate the formalism

developed in the preceding sections by applying it to two simple

cases: the antiplane-strain problem and the plane strain problem

in a cubic solid containing a Z-5 tilt grain boundary.

41
5. Application to the Antiplane-Strain Problem

For illustration, we apply the formalism developed in this

paper to the simple case of antiplane-strain in a bimaterial

composite containing a free surface. This case is of little

practical interest but has the important advantage of simplicity.

It can, therefore, serve the purpose of illustrating the

mathematical technique for calculating the SI Green's function

without unnecessarily clouding the essential features of the

technique. The experience gained in treating this simple case

should be useful to readers in solving a realistic problem. The

calculation of the SI Green's function for a more realistic case

has been given in the next section.

In the antiplane-strain case, all the 3x3 matrices such as y,

a, Q and of course the Green's function itself reduce to lxl

matrices, that is, pure numbers. Further, in this case the indices

a, |3, etc., will take only one value, that is, a = j3 =1. We assume

the solids A and B which constitute our composite model to be

cubic. This is the same model which has been considered in I and

II. Values of various parameters for this model are given below.

Some of these have been calculated in I and II; others can be

obtained quite simply by following the same method.

P« = P/3 = L

' (131)

r

A' B

(«) = *g'

B

= l/2dAB , (132)

crA' B(a) = (A B

= i/2 , (133)

Q

1 = Q1 = Q

IV

= QIV= d , (134)

a s a s '

v '

Q

11

= Q

11 = 2dD , (135) a s B

v '

Q

111

= Q

111

= ~2dA , (136)

a s A

v '

n

A,B = 1/2 ; (137)

42
KAA (tt/<3) = Ld/2; KAB (a,<S/ = LdA

'

KBA («,^) = R(a,/3) = Ld

g

,

K

BB (a,^) = P(a,/3) = -0.5 Ld , jA,B = r = _ L/2 ^

EAA

= d / 2dA ' EAB =

- 1

'

EBA = 1 EBB = d/ 2dB

'

a

a/3 = -1 = exp(TTL) ; = -1 = exp(-Tri) , (138)

b

a/3 = d

a(3 = 1

'

d = dA

- d

B

, (139)

d A _ = c

A,B / [ cA

+ c

B

] , (140) A,B 44 '

L 44 44 J

*

v '

dA

+ d

B = 1 . (141)

Substituting the values of the constants in eqs (131) -(141)

into eqs (117) and (118) , we obtain the following expressions for

the displacement Green's function for a unit line force applied at

x' in the UHP :

i. x in UHP

G(x,x') = (-l/2TTdA

) Re[ln(z-p) + d ln(z-p)]

(142)

43
ii. x in LHP

G

B

(x,x') = (-1/7T) Re[ln(z-p)] + t-[u

A

[- \ )

+ U*( f ]]

- S^l !)"*(-? )]

• (143)

Following the procedure given in section 2, or directly from

eqs (15) and (16) , we obtain the following expressions for the

stress (r ) Green's functions:

TA(x,x') = - Re — x

2 71 L (z

+

- p) (z - p )

-|— [

HA[(z+e)/i] - HA[-(z + C)/L

[

HA[-z/L] - HA[z/L]] + dA [H

B [z/L] - HB

[-z/c]

and

T

B

(x,x') = - Re [

1

_^ ] - d

B [

HA[-z/ L ] - HA[z/i]

| [

H

B [-(z+e)/i] - HB [(z+e)/t

E

B

[z/l] - H

B [-z/l]

j

, (144)

44
where, in view of eqs (131) and (A. 5),

z = x12+ l x = x + lv.J (145)

As in section 2, we obtain the following equations for the

two H-functions by using the boundary condition given by eqs (4)

and (5) or directly from eqs (25) and (26)

:

HA

(y + lc) - HA(y - lc)

j =

27TL

k y + Lp y - Lp

and

+ d

L y + tp y - Lp J

J

(146)

[

H

B

(y + lc) - H

B

(y - lc)

L y-Lp y+Lp J

(147)

Following the method given in Appendix I-B and the steps

given in section 2, we obtain the following solution of eqs (146)

and (147) .

V^q) = s(-Lp) - s(Lp) + d [ s(-Lp) - S(L P ) ], (148)

V

B

(q) = 2 d

B

L exp (irq) [ s(-Lp) - s(l P ) ] , (149)

H

A

(z) = f Aq) z Lq" 0 * 5 dq , (150)

—oo

and

45
H*(Z) = f

V

B

(q) z Lq-°' 5 dq , (151)

where

s(p) =

2jf- [exp(-27rq) / E

2

(q)j p

Lq °* 5

. (152)

Now we consider evaluation of the integral in eq (150) with

(q) given by eq (148) . The integrand has second-order poles at

all values of q [see eq (C.17) for definition of q ]. We can use

n n eq (C.24) for evaluating this integral by putting all Q = q and

r n accounting for the second-order poles as described in Appendix

I-C. This would give a series in q which can be easily summed. We

n can also evaluate this integral directly by making the

substitution t = exp(-2irq). The limits of integration over t

become 0 to », The integral can then be carried out by using the

same contour as used in Appendix I-C for evaluating I(z) in eq

(C.l). The integral for H (z) in eq (151) can be evaluated in a

similar manner. Thus we obtain the result

H

p

(z) = ~ 2nT [S^'"^) " S ( Z ' L P)

+ d S(z,-l P ) - S(z,Lp) (153)

and

H

P

(Z) = " WZ [S < Z ' L P) ~ S(z,-Lp)J , (154)

where

S(z,p) = 2ir(z - p ) ln (Z/P) * (155)

46
In eqs (153) and (154) , the subscript p on the H-functions

indicates that they are particular solutions of the Hilbert's

equation. We have to add the homogeneous solution to the

particular solution unless the corresponding force functions

satisfy the condition given by eq (12) . It is therefore necessary,

as discussed in section 2, to evaluate the integral of the force

functions which is given below.

In the present case, in view of the Plemelj relation, as

given by eq (B.42), the force functions are directly given by eqs

(146) and (147) . The integral of the force functions can now be

easily carried out. The result is

^00 0

I

F =

J

F

A

(y)dy + F

B (y)dy = 1 , (156)

^ ^00

where we have to take the limit y = «.

00

In deriving eq (156) , we have accounted for the fact that in

A

the definition of H (z) in eq (146) , the cut is on the positive

real axis. Hence, for the discontinuity in H (y) at the real axis,

we have used the relations

i = exp(TTL/2) , -L = exp (3TTL/2 ) ; (157)

ln(yw + ip) = In y , ln(y - up) = In y + 2ttl; (158) 00 00 00 00

where p has a positive real part and mod(p/yOT

) « 1.

Similarly, for H (z) in eq (146) , the cut is on the negative

side of the real axis. In this case, therefore,

l = exp(TTL/2) , -I = exp(-TTL/2) ; (159)

and

47
ln(-y v J + ip) = In y + ttl , ln(-y - ip) = In y - til; (160) 00 oo

N -*oo J

00 x '

We find from eq (158) that the force functions given by eqs

(148) and (149) do not satisfy the constraint given by eq (12) . We

therefore add the following force functions, which are solutions

of the homogeneous Hilbert equation, that is, eqs (143) and (144)

with their RHS equal to zero.

F^(Y) = -e/7r(y2

+e

2

) (161)

and

Fg(y) = -e/7r(y2

+e

2

) , (162)

where, in both eqs (161) and (162) , the limit e = + 0 has to be

taken. The two force functions have a delta function dependence on

y. In the limit e = 0, they approach zero for all values of y

except for y = 0. Their integral is finite, and they make a

nonvanishing contribution to the H-functions. The corresponding

H-functions which represent the homogeneous solution of the

Hilbert equation, are given below.

H

o<

z

> = 2¥T z_1 < 163 > and

H

o<

z > - WZ z_1 •

< 164 >

The integrals of the force functions as given by eqs (161)

and (162) can be evaluated as follows

I

F(J = Lim [ J

Yo

° F^(y) dy +

J°

F^(y) dy

J.

(165)

Yoo = 00 0 ~yoo

48
Using eqs (157) — (160) , and taking the limit e = + 0, and =

oo, we obtain

In writing the force functions given by eqs (161) and (162)

A B the total integral has been divided equally between F and F but

this division is quite arbitrary. Since both the force functions

are applied at the same point, that is, y = 0, what matters is

their total resultant. Thus, in general, the coefficient in eq

(161) could be taken to be any number—say £ instead of 0.5,

provided the coefficient in eq (162) is taken to be l-£.

From eqs (156) and (166) we see that if the force functions

given by eqs (161) and (162) are added to those given by eqs (143)

and (144) , the constraint given by eq (12) will be satisfied.

These forces will contribute zero stress at the free surface,

since they are solution of the homogeneous equation. Thus the

final solution which will satisfy eq (12) and all the other

boundary conditions is given by

A R A R A R

HA'*(Z) = h£'*(z) + H^' tJ (z) ; (167)

A B where the particular solutions H ' (z) are given by eqs (153) and

P A B

(154) and the homogeneous solutions H

Q

' (z) by eqs (163) and

(164) . Although the two expressions given by eqs (163) and (164)

look similar, the former has to be evaluated by taking a cut on

the positive half of the real axis and the latter by taking the

cut on the negative half of the real axis.

The stress-1 Green's functions are thus given by eq (144) in

the UHP and eq (14 5) in the LHP with the H-functions given by eq

(167). The displacement Green's functions are given by eqs (142)

and (143) in terms of the U-functions which, according to eqs

49
(123) and (124), are obtained by integrating the H functions. The

result is

ljA(z) = " 2Wr[v(z '~ Lp) " V < Z ' L P")

+ d ^v(z,-Lp) - v(z,tp)J] + JL_ln(z) (168)

and

.B

uB(z) = "

t?T [

V ( Z ' L ^) " v(z,-ip)j + _L_ln(z), (169)

where

l r 2 v(z,p) = - 0.5 (In z) + In z ln(z-p) - In p ln(z-p)

" I ~h ( p/

z)n

1

(for l z/p l > 1 ) (170)

n n J

and

v(z,p) = p In p In z + In z ln(p-z) - In p ln(p-z)

+ X

-^2 ( z/p) n

(for |z/p| < 1 ) , n n -1 (171)

and where the sum is over all positive integers n (n—1,2,..). The

last logarithmic term in eqs (170) and (171) corresponds to the

homogeneous solution.

The stress-2 Green's functions can be obtained from the

displacement Green's functions by using the method given in

section 2 or directly from eqs (125) and (126) . These are given

below.

50
i. x in UHP :

T

A

(x,x') ^Re -1 + —*

*

2 TT L ( Z - p) (Z - P )

-I

+

[

HA[(z+e)/L] + HA

[-(z + c)/l]

" -^f- [

HA

[-z/L] + HA[z/L]

+ LdA |^ H

B [Z/L] + H

B

[-z/L]J (172)

ii. x in LHP :

T

B

(x,x') = - •^

B- Re

[ ^_

1

_^ J

+ Ld

B [^ HA

[-z/l] + HA

[z/L]

+ —^-Jh

B [-(z+c) /L] + H

B [(z+e)/L

+ H

B [z/l] - HB [-z/l]

J

(173)

The calculated values of the displacement field u(x)

,

stress-1 component t^^x) an<3 the stress-2 component r

32 (x) have

been shown as 3-D plots in figure 1-2 through 1-4 respectively. In

these calculations d, and d_ have been taken to 0.8 and 0.2 A B

respectively. The unit line force has been applied at x' = (3,

1.5). The vector x' is chosen so that it does not fall at any

point of the mesh where the Green's function is calculated. This

is to avoid the characteristic singularity in the Green's function

at x=x' . We see from figure 1-2 through 1-4 that t 31 is 0 at x =0

and u(x) and x

22 are continuous at x

2

=0 whereas Z21 ^ s

discontinuous at *

2

=0.

51
An important point concerning the numerical calculations: In

deriving the expression for H (z) , we have taken the cut on the

positive half of the real axis. This implies that the argument of

the variable changes by 2ttl when it crosses the real axis on the

positive side. This has to be properly accounted for in the

computer programming. For example the complex log should have the

following values in the limit c = +0 where y is real and positive.

ln(y + lg) = In y ; ln(y -lc) = In y + 27tl (174)

and

ln(-y + lc) = In y + til ; ln(-y - lc) = In y + til. (175)

On the other hand, in normal Fortran on most computers, the branch

values are taken as

ln(y + lc) = In y , ln(y -lc) = In y; (176)

and

ln(-y +lc) = In y + til , ln(-y - lc) = In y - til . (177)

In the derivation of H (z) , the cut has been taken on the

negative side of the real axis. In this case the argument of the

variable would change from til to -7tl as it goes across the cut.

The branch values of the log function for this cut will be the

same as given by eqs (176) and (177) .

In the numerical calculations for the H functions, it is

therefore necessary to program the computer so that the

appropriate branch values are taken, that is, as given by eqs

(174) and (175) for the HA

- function and by eqs (176) and (177) for

the H -function.

52
6. Application to a Cubic Crystal Containing a Z-5 Grain Boundary —Plane-Strain Problem

As a further illustration of the formulation given in section

3, we apply it to calculating the SI Green's function for a cubic

crystal containing a Z-5 tilt boundary. We have chosen this

particular case for illustration because of the current interest

in the structure of grain boundaries and because its high symmetry

simplifies the calculations considerably. The SI Green's function

as calculated in this section can be applied to a variety of

problems concerning the elastic properties of the grain boundaries

such as dislocation pile-up and crack propagation.

A tilt boundary in a solid can be visualized as a composite

solid in which two half solids are welded together along the

boundary and the crystallographic axes of the solids are oriented

with respect to each other. In particular, the S-5 boundary in a

cubic solid can be modelled as follows (see, for example,

[24-25]). Consider two pieces of the same cubic material with

their crystallographic axes parallel to the coordinate axes.

Rotate one of the pieces about the Z-axis by an angle 9 = tan 1

(3/5) while keeping the other piece fixed. Weld the two

pieces together so that they form a planar interface. This

interface will be the S-5 grain boundary.

The cubic solid, which we have chosen for our calculations,

is stainless steel. The Green's functions for this solid without

the free surface have been calculated in I . In this paper we use

the same model as described in I. We label the solid which has its

crystallographic axes parallel to the coordinate axes as the solid

A, that is, in the UHP. The solid, which has been rotated as

described in the preceding paragraph, is labelled as the solid B,

that is, in the LHP.

53
The values [18] of the elastic constants of stainless steel

are taken as follows (in units of c..): c, , = 2.2, c, _ = 1.3, c MM 44' 11 12 44 = 1.0. These will be the elastic constants of the solid A in our

model. The elastic constants of the solid B, that is, the rotated

solid, can be obtained by using the rotation transformation law

for the fourth-rank tensors as given below:

B A

cijkl

= S

ii' S

jj

,Skk/S

ll' c

i'j'k'l' ' (178)

where S is the matrix of rotation. Its elements are

S ll = S 22 = COS 9 < 179 > and

S 12 = - S 21 = Sin 6 , (180)

where, for a Z-5 boundary, 8 = tan 1

(3/5).

The calculation of various parameters for the present model

has been described in I. Here we only give the values of those

parameters which are needed for the present calculation. As given

in I, in the present case, a and /3 which label the roots of eq

(A. 19), take only the values 1 and 2. The third root, a = 3,

corresponds to the antiplane-strain mode which has already been

discussed in the preceding section.

The most basic parameters required for the Green's function

A B calculation are pft

'

, which are the roots of eq (A. 19), and the

matrices y(a) and g;(a) as defined by eqs (A. 17) and (A. 15)

respectively. The expressions for these parameters are given

below. Other matrices can be easily calculated in terms of these

matrices by using various relations given in section 2. In the

expressions given below the superscripts A and B are not shown for

notational brevity.

54
Pl = [(1 - K

2

)/ 2] 1/2

+ L[(l + K

2

)/ 2] 1/2

p2 = - [(1 - K

2

)/ 2] 1/2

+ L[(l + K

2

)/ 2] 1/2,

Z ( 1 + C D

2

) .

k (a) k (a) = — Z 8 d *21^ ' a

p 0 pa

'

*22 (a) z

a < C + Pa

2

) .

a

ll (a) Z

a P«< 1 * p 0

+ t Pa" >

'

a

l2 (a) Z [ ( 1 - /3 n )p

2

+ C ] /

<r 21 (a> Z ( -1 + B„ - C P

2

) ,

a

22 (a) c AA Z p ( C

2 ~ £2

+ /3_ 44 a *a v s 1 0 '0

2

) a

where a = 1 or 2 , and

Z = - ^

[1 - K

4 ]" 1 / 2

, « 4cllPa

K

2 = 1 + (5 Q / C) O0

+ 0.5 5

Q

) ,

^ = C ll /

C

44 '

*0 = < C

12

+ C44^C44 '

and

5

0 = (C11 " C12 " 2C44 ) /C44 V

The point of application of the unit force is taken to be the

55
same as in the previous section, that is, at x' = (3.0,1.5). As

mentioned in the previous section, for numerical convenience, the

value of the vector x' should be chosen so that it does not fall

at any point of the mesh where the Green's function is calculated.

This is to avoid the characteristic singularity in the Green's

function at x = x'

.

The calculated values of the the displacement u^x) and the

stress components z^ix) , "^l^-^

= x

\2

^ ' and T 22^-^ nave been

shown as 3D plots in figures 1-5 through 1-8 respectively. The

plot for u

2 (x) is similar to that of u1

(x) and, therefore, has not

been given.

We see from figure 1-5, 1-7, and 1-8 that u

1 (x), t12 (x) and

t

22 (x) are continuous across the interface at x2 = 0 as required

by the boundary conditions given by eqs (2) and (3) . The stress

component x \\^^ i- s n°t continuous across the interface as shown

in figure 1-6. The components T 21^-^ anc^ T ll^-^ are zero a^ ^e

free surface, that is, at x

1 = 0 as required by the boundary

condition given by eq (4) and (5)

.

The singularities in the stresses will be discussed in the

following paper on the generalized plane-strain problem. As

mentioned in section 3D, the singularities in H(z) arise from the

homogeneous part of eq (36) and, in general, do not depend upon

the RHS of eq (36) . As we shall see in the following paper, the

matrix M(q) is exactly the same in the present plane-strain case

and the generalized plane-strain case.

It may be emphasized again that, in using the formulation of

this paper for numerical calculations, the values of the various

functions such as ln(z) and z

L(^ over the upper and lower branches

of the cut have to evaluated carefully according to the scheme

given by eqs (174)-(177).

56
7. Summary of Main Results

In this paper we have obtained the displacement Green's

functions and the stress Green's functions for a general

anisotropic bimaterial composite solid which has a plane interface

and a free surface normal to the interface. These are applicable

to a variety of plane-strain problems and also account for

displacements and forces normal to the plane.

We find that the stress field is singular near the

intersection of the free surface and the interface. The

singularities are of the type r

11 and [ln(r)] n where 7) is a

positive number between 0 and 1 and n is a positive integer

between 0 and 6. The actual values of 7) and n and the coefficient,

i.e., the weight of the singularities, depend upon the material

properties of the solids. In addition to the singularity, an

oscillatory behavior of cos or sin[ln(r)] type, may also be

present in certain cases.

A method of solution of the generalized inhomogeneous vector

Hilbert problem has been given. The method is simple and

convenient for applications to different problems in the stress

analysis of materials involving cracks and interfaces in which the

inhomogeneity and the nonsingular part of the Hilbert' s kernel are

of the form l/(y - p) where y is a real variable and p is complex.

For this purpose, we have introduced a complex transform which is

of the form y

L<^ 0,5 where l = v'(-l) and q is a real variable. This

transform seems to be particularly interesting, because it is an

eigenfunction of the Hilbert' s kernel and is orthogonal over both

variables y and q in the domain 0 ^ y < « and -« < q < w.

57
Appendix I -A

Elastic Green's Function for an Infinite Bimaterial Composite

In this appendix we quote the expressions for the elastic

Green's function for a bimaterial composite as obtained in I and

used for calculation of Green's function for a cracked solid in

II. The expressions given below have been taken from II and are in

a form which is more convenient for present application. The

coordinate system used for these expressions is the same as shown

in figure 1-1. The indices a and /3 take the values 1, 2, or 3.

The Green's function given below is the displacement Green's

function G(x,x'), which gives the displacement field at the point

x when a unit line force is applied at the point x' . As in I and

II, we distinguish between four different cases corresponding to x

and x' being in UHP and LHP.

1. x and x' both in UHP ( x * 0 ; x ' * 0 ) :

a

(A.l)

a/3

2 . x in LHP, x' in UHP ( x

2

s 0 • x

'

' X

2

2: 0 ) :

Gq (x,x') =

1 ~B, v ^11 -, ,-B -A , ^ 2_ *

(0°

2f?

ln(z oT P

/3

> (A. 2)

a/3

58
3. x in UHP, x' in LHP ( x

2

£ 0 ; x

2

' ^ 0 ) :

GqB

(x,x') = . 1

r

A

(cx) q"1 ln(z^- p

A

) , (A. 3)

a/3

4 . x and x

7 both in LHP ( x

2

^ 0 ; x

2

' £ 0 )

:

GBB (x,x') = - I £ J

B

(a) In (5*- p

B

) a

1 V "B, v ^IV , ,-B B.

+

n ) I ( a ) 2|3 ln ( z a

~ <V ' ( A - 4 )

a/3

where

A, B A.B /T. ... V = X

l

+ V X

2 ' (A» 5)

p

A,B = xl + pA/B xl , (A. 6)

Qa = M

[

a

A

(oc) - e

B

(i

B

) _1

r

A

(a)] , (A. 7)

Q* 1 = N

[

a

A

(a) ~ /(is 1 " 1

I

A

(«)] » (A. 8

)

Ql

11 = M

[

ciB

(cO - (Ig)"

1 *V)] , (A. 9)

9

1V

= n

[ e

B

(«) - ^s^s 1 " 1 *B

< a

>] ' ( A - 10 >

, A. -1 T ~B ,-B. -1 A, A. -ll" 1 „„, M = (I s ) [ 2S ^s* " 2s ( 2S ) J

' (A. 11)

and

59
,-B. -1 T -B ,-B. -1 A, A.-ll

-1

(A. 12)

(Note. The following formulas are valid for either superscript, A

or B, corresponding to the solid in UHP or LHP, respectively.)

Is

= I

' a

(A. 13)

Zs = I 2t«) / (A. 14)

a(a) = L

1Z

(a) r(oc) , (A. 15)

II L

ij«*> =

C i2kl + PaC i2k2 ' (A. 16)

K • • (a)

r

ii<q2

>

aq. <p« - p«' n

a*/3

(for q2 = q pj , (A. 17)

• (g) is the ij cofactor of the Christoffel matrix A defined by

A

ij (3) = cikjl qk qx (A. 18)

qn and q_ are components of the wave vector g and p is obtained

such that q2 = pa q1is a root of the equation,

I |A(q) ||=0, (A. 19)

|A| denotes the determinant of the matrix A and a is the

coefficient of q. in | |A| | , and

Im (p ) > 0 . (A. 20)

60
In analogy with eqs (A. 13) and (A. 14), we define the

following quantities as sums over a:

<4

,11, in, iv = y Qi, ii, in, iv

_ (A 21)

a

Several useful relations between the various parameters of

the Green's function (such as j, cr, Q ) have been given in

Appendix A of II. Two more relations which involve the constants

introduced in this paper are

A

D

r = 0 , (A. 22)

where

n(a)/Pa = -v(cc) = J(a) , (A. 23)

*ik = C ilkl + PcJ C ilk2 + C i2kl>

+ Pa C i2k2 ' (A. 24

)

61
Appendix I-B

Solution of Generalized Inhomogeneous Vector Hilbert Problem

In this appendix we give a solution of the generalized

inhomogeneous vector Hilbert problem. The corresponding scalar

problem, its properties and solution have been discussed in detail

in the excellent treatise by Muskhelishvili [20] . The vector

problem has been discussed by Vekua [21] . Our method of solution

as given below is different from those given in [20,21] and is

based upon the use of a complex transform. This method should be

particularly convenient for those types of kernels and

inhomogeneities which occur in problems of stress analysis in

solids containing surfaces, interfaces, and cracks.

First we consider the scalar problem as defined below:

D H(y+LC) + 6 H(y-LC) + £

C

a

H(a

a

y) = f(y), (B.l)

a

where D, C

a

and a^ are complex numbers, a is a summation index

which takes an arbitrary but finite range of values, t and y are

real variables which vary between 0 and », c is a small positive

number which approaches zero in the limit, f (y) is a known

function, and H(z) for any complex z, is defined as

CO H(z) =

1

2ttl

0

F(t) dt . (B.2) t - z

The integral in eg (B.2) has to be interpreted as the

principal value integral for z on the real axis when the integral

is singular. The first two terms in eq (B.l) correspond to the

62
principal value of the integral in the definition of the Hfunction. The function F(t) in eq (B.2) is the unknown function

which has to be determined from eq (B.l).

In eq (B.l), f (y) is the inhomogeneity , which is zero for the

homogeneous Hilbert problem. If C = 0, this equation will have no

nonsingular term in the kernel and is then called the special (not

generalized) Hilbert problem (the prefix 'special' is usually

omitted) . The particular form of the nonsingular kernel as given

in eq (B.l) is of the type 1/ (t-ay) where a is complex, and is of

special interest in problems concerning stress analysis in solids

containing surfaces, interfaces, and cracks.

We solve eq (B.l) by using a complex transform method. For

this purpose we introduce the following function of two real

variables, q and y (0 ^ y ^ » and -co < q s co) :

T(q,y) = y

Lq-°- 5

. (B.3)

By using the values of the branch points as given in Appendix I-C,

we can verify that T(q,y) is an eigenfunction of the homogeneous

Hilbert problem given by the LHS of eq (B.l). Further, the

functions T(q,y) are orthogonal in q as well as y space as given

below:

.00 T(q,y) T(q',y) dy = 2tt 5 (q-q' ) (B.4)

0

J

and

.00 T(q,y) T(q ,y') dq = 2tt 5 (y-y'

)

(B.5)

-co

63
These orthogonality relations can be verified by using the

following standard definition of the delta function:

5(y) =

2TT exp(iqy) dq (B.6)

The integral in eq (B.4) can be evaluated by using the

substitution x = In y. In writing eq (B.5), we have used the

following property of the delta function of an arbitrary

dif f erentiable function p(y) of y:

5[P(Y)] = S(y-y

0 ) [ Hv^l

1

' (B * 7)

where y = y Q

is a real root of the equation p(y) = 0.

We now write F(t) in terms of its transform v(q) as

n

1 _ _ F(t) = v(q) t iq" dq . (B.8)

— 00

Substituting for F(t) from eq (B.8) in eq (B.2), and carrying out

the integration over y as shown in Appendix I-C (see eq C.7), we

obtain the following result for H(z)

:

H(Z) = f V(q) z iq~°' 5 dq , (B.9)

-co

J

where

V(q) = v(q)/E(q) (B.10)

and

64
E(q) = l/[l+exp(-27rq) ] (B. 11)

As described in Appendix I-C, in deriving eq (C.7), we have

taken a cut on the positive real axis (see figure 1-9) . The values

of argument just above and just below the real axis are therefore

taken to be 0 and 2tt, respectively, as given in eqs (C.3) and

(C.4). Thus, we get the following, in the limit e = +0,

00 , H(y+ie) = V(q) y

lc5" 0 * 5 dq (B.12)

-or

and

00 . H(y-te) = -[ V(q) exp(-27rq) y

lc3" 0 * 5 dq . (B.13)

In order to solve eq (B.l), we use the representation of H(z)

as given by eq (B.9) and use eqs (B.12) and (B.13) for the first

two terms. Then we multiply both sides of eq (B.l) by T(q,y) , i.e., y

LC

* ®'

, and integrate over y from 0 to oo. Using the

orthogonality relation as given by eq (B.4), we obtain the result

V(q) = [M(q)

j" 1 N(q) (B.14)

where

M(q) = [ D - D exp(-27rq) + £

C

ft a^q °* 5

] (B. 15)

and

00

x= t \ -Lq-0.5 , f(y) y dy (B. 16)

0

J

65
In many cases of practical interest such as stress analysis

in solids, the function f (y) is of the form

f (Y) = XR

« / [y " Pa

]

' (B ' 17)

a

where pais complex. In such cases, the integral on the RHS of eq

(B.16) can be easily obtained by using the same contour as used

for deriving eq (C.7) in Appendix I-C. The result, which can be

obtained by replacing q by -q (not i by -l) in eq (C.7) is given

below:

=

L exp(-27iq) j

-iq-0.5

E(q) L a

Vl a' v '

Equation (B.14) along with eqs (B.8), (B.9) and (B.10) gives

the particular solution of eq (B.l). In order to obtain the

general solution, we have to add a solution of the homogeneous

part of eq (B.l) to the particular solution. The homogeneous

solution can be written in terms of Q as

r

lQ -0.5

H (z) = £h(Q ) z r

, (B.19)

r

where h(Q

r

) are arbitrary constants; and Q , which may be complex,

is a solution of the following equation:

M(Q

r

) = 0 . (B.20)

The arbitrary constants can be determined with the help of the

boundary conditions on F(t).

Now we consider the vector problem which arises if there is

more than one unknown function F(t) and, of course, an equal

66
number of equations. In this case we regard eq (B.l) as a matrix

equation. If there are n unknown functions and n equations, then D

and in eq (B.l) will be n x n square matrices, whereas H(z) , f (y) will be n x 1 column matrices or n-component vectors.

Similarly, M(q) and N(q) will be n x n and n x 1 matrices

respectively. Equation (B.20) will become a determinantal equation

and h(Q

r

) will be an eigenvector of M(Q

r

) corresponding to a zero

eigenvalue. One of the components of this eigenvector will remain

arbitrary, which has to be determined from the boundary

conditions. With these modifications, the solution given above

remains valid for the vector Hilbert problem.

As a special case of the general method, we solve the vector

Hilbert problem for two unknown functions. This solution is used

in this paper (section 2) , where the unknown functions are the two

A A

force functions F and F . As in section 2 , we label the unknown

functions with A and B and also label the corresponding matrix

elements of eq (B.l) by the same indices. Equation (B.l) now

becomes a set of two simultaneous equations as given below:

DAA AA

a

(B.21)

a

and

+ DBB BBH (y-LC)

(B. 22)

a

67
The H-functions are now defined as

"a

F

A

(t)

t - z dt (B.23)

and

H

B

(Z) =

1

2TII

-oo"

° F

B

(t)

t - z dt (B.24)

As in eq (B.8), the force functions are written in the

following form:

FA(t) = fvA

(q) t^"0 - 5 dq

and

—00

F

B

(t) = f v

B

(q) |t| i(J- 0 - 5 — 00 dq

(B.25)

(B.26)

The corresponding expressions for the H-functions are

HA

(Z) =-- J

V

ft (q) z iq-0 .

5

dq (B.27)

and

H

B

(Z) = fV

B

(q) z^5 dq . —00 (B.28)

Using eq (B.25), the integration for H (z) in eq (B.23) is

carried out as in eq (C.7). The result is identical to that given

by eqs (B.9) and (B.10) except that the superscript A has to be

added to the relevant matrices. The corresponding integral for

gH (z) has also been evaluated in Appendix I-C [see eq (C.13)].

68
In case of H (z) , we have taken the cut on the negative real

axis. The argument of z is taken to be n just above the real axis

and -7i just below it. The corresponding values of the H-functions,

as obtained by using eqs (C.ll) and (C.12) are given below: (in

the limit c = +0; -« < y < 0) :

00 H

B

(y+LC) = -lJ V

B

(q)exp(-7rq) |y|

lc*- 0 - 5 dq (B.29)

—00

and

H

B

(y-ie) = l\ V

B

(q)exp(7rq) | y | . 5 dg # (B.30)

—00

To solve eqs (B.21) and (B.22), we proceed as before, except

that we multiply eq (B.21) by y~ L(*~ 0 ' 5 and eq (B.22) by |y|" Lq" 0,5

and then integrate over y between the limits given in eqs (B.21)

and (B.22) to ensure the orthogonality of the transform. Thus, we

obtain the following particular solution for the force functions

and the corresponding H-functions:

V

A

(q) = VA

(q)/E(q) = NA

+

^ , (B.31)

V

B

(q) = iv

B (q)exp(-7rq)/E(q) = M*A

NA

+ Mj^ N

B

, (B.32)

where M

1 is the inverse of M. The matrix elements of M and N are

M

AA " DAA " 5AA ««P(-2"q) + £

^

(aj ' 5

, (B.33)

a

M

AB

a

I WViq"°- 5

< (B. 34)

69
m = V c r-c ) i(J-°- 5

BA L aBA * c

a

; a

(B.35)

M

BB

=" L [

DBB exP(-"q) - DBB exp(Trq)] + J

CaBB ) iq"°' 5

f (B.36)

J a

.00 NA (q)

1

27T

_ , . — iq— 0 . 5 , fA

(y) y

H dy (0 * y * oo)

, (B.37)

and

N

B

(q) _1

27T

oo

fB

(y) y|

L(3" 0 - 5 dy (-oo < y < 0) (B.38)

Equations (B.31) and (B.32), when substituted into eqs (B.27)

and (B.28), give the particular solution for the H-functions. The

homogeneous solution is given by eq (B.19), where h(Q

r

) is now

identified as an eigenvector of the matrix M(Q

r

) for a zero

eigenvalue with Qras a solution of the following determinantal

equation,

| |M(Q

r ) ||=0. (B. 39)

The final solution is the sum of the particular solution and

the homogeneous solution as given below:

rA(z) = r

v

A

(q) z±q 0,5 dq + 1 w z r

-oo

J r

iQ -0.5

(B.40)

and

r r°° -irr-n r ^ lQ -0.5

H

B

(z) = V (q) z iq U * D dq + £hg(Q

r

) z r

. (B.41)

-oo

J r

70
One of the two components of the eigenvector h, (Q ) and h_,(Q_ ) is

arbitrary and has to be determined from the boundary conditions or

the constraints imposed on the force functions.

The force function can be obtained directly from the

H-functions by using Plemelj relation [20]

F(y) = H(y+te) - H(y-Le) . (B.42)

In the present case, eq (B.42) may be easily verified by using the

transforms of F and H as given by eqs (B.31) and (B.32).

In many problems concerning stress analysis in solids, the

inhomogeneity f(y) is of the form given by eq (B.17). In the

present case, we write the two components of the vector f(y) in

the form,

fA

(y) = I

R

a / [y " pa

] (0 " y "

°°

) (B * 43)

a

and

fB

(y) = £ R* / [y -

pjj|] (-00 < y < 0) . (B.44)

a

The values of the integrals for N A (q) and N„(q) in eqs (B.37)

and (B.38) can be obtained from eqs (C.7) and (C.13) by replacing

q by -q and not l by -l. The result for NA

(q) is identical to that

given by eq (B.18) except that the superscript A has to be added

to R and p on the RHS. The result for ND (q) is

V«> -- SSE^X R

2

(Pa

>- tq -°- 5

• (B.45,

71
Appendix I-C

Evaluation of Integrals

In this appendix, we shall evaluate certain integrals by

using contour integration. These integrals are required for

deriving various formulas in the paper. We shall also obtain the

force function and its integral to ensure that the constraint

given by eg (12) is satisfied.

required for calculating the H function from the force function

First, we shall evaluate the following integral, which is

.red for calculate

in eq (B. 2) or (B. 23) .

00 . _ t iq-0 .

5

I

A

(q,P) =

t - p

dt , (C.l)

where p is complex and q is real.

To evaluate I (q,p) , consider the integral of the following

complex function of variable z :

iq-0.5

Z = — , (C.2)

z - p

where z is complex. We integrate Z over the contour shown in

figure 1-9. The contour consists of an outer circle C with a cut

on the positive real axis. The radius of the circle approaches

infinity in the limit. The argument of z is taken to be 0 at the

upper branch (just above the cut) and 2n at the lower branch (just

below the cut) ; that is,

t = t + ie = t (C.3)

72
and

t_ = t - lc = t exp(27TL) , (C.4)

where t is the real part of z.

Thus we obtain from Cauchy's theorem,

Z iq-°- 5

i> . .A

C Z - p

dz + E(q) I

A

(q,p) = 2ttl Residues (C.5)

where

E(q) = 1 + exp(-27rq) . (C.6)

The first term in eq (C.5), denoting the contribution of the

outer circle, vanishes in the limit when C extends to infinity.

Since the integrand has only one pole inside the contour at z = p,

we obtain from eq (C.5)

i.q-0.5

TA, , 2ni p ^ -7\ 1 (q ' p) =

Efq) * (C ' 7)

In order to calculate the value of H (y+te) and H (y-te) , we

need the value of I(q,p) just above and just below the real axis.

Using eq (C.3) and (C.4), we obtain the following values of the

integral over the upper and the lower branches of the cut:

I (q,y+Le) =

2ttl y

Lq-0 .

5

E(q) (C.8)

and

I (q,y-te) =

2?Tiexp (-27rq) y

E(q)

(.q-0 . 5

(C.9)

73
We now evaluate the following integral, which is required for

g calculating H (z) in eq (B.24):

I

B

(q,P) =

-co

°

|+-l ic3- 0 - 5

-L^J dt . (C.10)

t - p

g To evaluate I (q,p) , we consider the same function as defined by

eq (C.2) and evaluate its integral on the contour given in figure

1-10, which has a cut on the negative real axis. The argument of z is taken to be n on the upper branch and -n on the lower branch

so for -co < t ^ 0,

t = t + lc = |t|exp(7Ti) (C.ll)

and

t_ = t — LC = 1 1 1 exp (-7TL ) . (C.12)

Proceeding as before, we obtain

I

B

(q,p) = -

2 "eXP( l(V

tq~ 0 " 5

• (C.13)

g The values of I (q,p) just above and below the real axis as

calculated by using eqs (C.ll)- (C.13) are

I

B (q,y+tc) =2TIcexp(-21Tq) | y | '-q- 0 - 5

(c 14)

and

i

i Lq-0 .

5

T B. , 27TL y ^ ,~ i<=\ I (q,y-LC) = . (C.15)

74
We now evaluate the integral over q required for the

calculation of the H function. As we see from eq (54) , we need to

evaluate an integral of the form,

CO P(q) dq z Lq~ 0 ' 5

, (C.16)

E(q)D(q)

— 00

where D(q) is the determinant of the matrix M(q) and P(q) is an

analytical function of q with no poles.

To evaluate the integral in eq (C.16), we choose a

semicircular contour either in the UHP or in the LHP as shown in

figure 1-11. In the present case, the contribution of the integral

over the large semicircle vanishes in the UHP if mod(z) > 1 and in

the LHP if mod(z) < 1. Accordingly, we choose the contour in the

UHP if mod(z) > 1 and in the LHP if mod(z) < 1.

Let q

n

and Q

r

denote the zeroes of E(q) and D(q)

respectively. The integrand in eq (C.16) has two sets of poles;

those at q will be called Set I and those at Q Set II. Thus

r

qn = (n + 0.5)i , (C.17)

where n is an integer- which can be positive, negative or zero and

are the roots of the determinantal equation,

D(Q

r

) = 0 . (C.18)

We assume that all Qr

are distinct, which should be generally

true. When they are not, it would be a higher order pole which can

be treated by the standard methods of contour integration.

Moreover, in general, some q^ will be equal to some Q^.. These

second-order poles can also be treated by standard methods.

75
However, presently our object is to write a general expression for

the integral with the view to identify its singularities and to

obtain an expression for the force function and its integral. For

this purpose, it would be convenient to introduce a small positive

imaginary part lc in q in the denominator of the RHS of eq (C.16)

as follows and, in the end, take the limit c = 0.

00

IH(z) = P(q)dq z iq 0,5. (C.19)

E(q-ie) D(q+ie)

—00

The two sets of poles in the integrand in eq (C.19) can now

be written as follows:

Set I. Zeroes of E(q)

.

q = qn

+ ic . (C.20)

Set II. Zeroes of D(q)

.

q = Qr

-ic . (C.21)

The allowed values of qnand Qr depend upon whether the contour of

integration (see figure 1-11) has been chosen to be in the UHP or

LHP and are given below:

I. Mod(z) > 1 (Contour in UHP).

n = 0 , 1 , 2 , . . . and

Im(Q

r

) > 0. (C.22)

II. Mod(z) < 1 (Contour in LHP).

76
and

Im(Q

r

) < 0 (C.23)

The residue of E(q-ie) at qnis 2tt and that of D(q+ie) at Qr

is D' (Q )

, where the prime denotes differentiation with respect to

the argument. Thus, we obtain the following value of the integral

by using Cauchy's integration formula:

IH

(z) ^

D(q

P(q +ie) z-( n+1 >-e

_ (q + 2ic)

n

v ^n '

iQ -0.5+c

P(Q

r

-ie )z r

*

2 " L Y. E(Q

r

-2; g )D'(Q

r

) * R < Z ' C

> <

( C « 24 ) r

where the sums in eq (C.24) are over those values of qnand Qr

which are not degenerate, that is, which represent only a

single-order pole; and R(z,e) denotes the contribution of those

values of n and r which represent a pole of second or higher order

in the integrand. The plus signs in eq (C.24) are to be taken when

mod(z) > 1 and negative when mod(z) < 1.

In the cases considered in section 3A and 3B, D(q) is zero at

the following two values of q which are part of Set I:

qo = Qo = L/2 (C.25)

and

q_1 = Q_± = -l/2 , (C.26)

which correspond respectively, to n = 0 and -1 in eq (C.17). These

values therefore represent second-order poles which will

contribute, respectively, when the contour is chosen in UHP and

77
LHP. From eqs (C.22) and (C.23) the pole at qQ

will contribute

when mod(z) > 1 and the one at will contribute when mod(z) <

1. We shall now evaluate R^(z,c) at these poles.

First we consider the former case, when mod(z) > 1. We use

the expansion formulas,

E(q

Q

-2L£) =-47TLe(l + 2nic) (C.27)

and

D(q

o+ 2 lc) = 2 lc D'(q

Q

) [1 + LCD" (qQ) / D' (qQ) ] (C.28)

Thus, we obtain the following for mod(z) > 1:

-1

Lim R (z,e) =

e->o

u

D' (qQ

) |_

P' (q0

) + LP(q

Q

) Ln(z)

^

P(q

0 )D" (qQ

)/D' (qQ ) +TrP(q

o

) (C.29)

The corresponding expression for R(z,e) when mod (z) < 1 is

formally similar to eq (C.29) and is written as

Lim R, (z , e) = —-

e->o

1

D' (qQ

)

P'(q_

1

) + LP(q_

1

) Ln(z)

1

2 P(q_1 )D" (q^/D' (q_ 1 ) +7TP(q_

1

) (C. 30)

We now calculate the integral required for the force

functions. An explicit knowledge of these functions is not

required for calculating the Green's functions for which, as we

78
see from eqs (116) and (117) , a knowledge of H-functions is

sufficient. However, we do need to calculate the integral of the

force functions in order to satisfy the constraint given by

equation (12)

.

Using the Plemelj relation, as given by eq (B.42), we see

that the two force functions can be written in terms of an

integral of the form,

CO E(q) P(q)dq | y

|

iq °* 5

. (C.31)

E(q-ie)D(q+ie)

— 00

In case of the force function in the LHP, the numerator of

the integrand in eq (C.31) will have a factor exp(Trq) which can be

easily included in the result derived below. The integral in eq

(C.31) can be evaluated by using the same method as described

earlier. The result is

P(q

n+ie)|yf (n+1) "C

I (y)= -2TTG

D(q + 2ic)

+ 2711

E(Q

r

-Lc)P(Q

r

-ic ) |y|

E(Q

r

-2(.e)D' (Qr

)

iQ -0.5+e

r

+ R(y/p,e) , (C.32)

where the allowed values of n and r are defined by eqs (C.22) and

(C.23), and R(y/p,e) accounts for any second order poles as the

corresponding R term in eq (C.24). The complex variable p is

contained in the function P(q)

.

The first term on the RHS of eq (C.32) (along with the last

term where it exists) gives the particular solution of eqs (25)

79
and (26) . The second term is a solution of the homogeneous

equation since D(Q

r

) is zero for all Qr

» The first term on the RHS

of eq (C.32) is vanishingly small in the limit e = 0. However,

this term can-not be neglected because it gives a finite

contribution to H(z) , as we can see from eqs (C.19) and (C.24).

The 1/y term which corresponds to n = 0 or q = i/2 , is actually a

delta-function term which is zero for all values of y except y=0.

This is similar to the term in eqs (158) and (159) for the

antiplane-strain mode.

The integral of the force function in the q representation

can be directly obtained from eq (C.32). We write the total

integral of the force function as a sum of those for the UHP and

the LHP as

it = -

+ I

B

(C. 33)

where

ip(y) <Jy

oo

—00

P

A

(q)dq

D(q) (iq+0.5)

iq+0 .

5

00

iq+0.

5

(C. 34)

.B with a similar expression for 1 . The corresponding limits on the

-B

y integration in I will be from -y to 0. In the end we shall

take the limits y = 0 and y = oo. J

0 oo

80
Now we evaluate the q-integral by using the contour

integration as in the case of integrals in eqs (C.19) and (C.31).

First, consider the integral of the lower y limit, the term

containing y . In this limit, as discussed earlier in this

appendix, we have to take the contour in the LHP as shown in

figure 1-11. Thus we find that the integral will contain only

positive powers of yQand therefore will be zero in the limit yQ =

0.

Now consider the integral of the term containing y . In this

case we have to choose the contour in the UHP. We have to include

the residue at Q , the zeroes of D(q), and at q Q = l/2, which

arises from the factor (q-L/2) in the denominator. For the moment,

ignore the pole at qQ = l/2 The integral will then contain terms

with (t.Q

r

+0.5) as powers of y^. In this case, as given by eqs

(C.17) and (C.18), Im(Q

r

) ^ 0. Obviously the terms for which Im

(Qr

) - 0.5 will also vanish in the limit y^ = ». Hence the only

terms which will contribute to the integral are those for which

0 ^ Im(Q

r

) < 0.5. (C.35)

Now we consider the special term for which q = l/2 for which

the exponent of y^ is 0. As shown in section 3A, the value of D(q)

is also zero at q = l/2. However it does not represent a

second-order pole because, as given by eq (97), P(q) is also zero

at q = l/2. Therefore, the integrand in eq (C.34) has only a

first-order pole at qQ = l/2. The contribution of this pole to the

integral can be written as

I

A

Q = 2n Lim [ -4$-]- < c - 36 > q=L/2 L J

Similarly, we obtain the following expression for the

integral of the force on the negative Y axis:

81
I

B

0 = 2tt Lim

q=L/2

(q) (C. 37)

Adding eqs (C.36) and (C.37), and using eqs (93) and (94) in the

definition of P(q) as given by eq (55) , we find that the numerator

of the limit will approach zero faster than the denominator. The

total contribution to the integral from the pole at qQ = i/2 will,

therefore, be zero.

For satisfying eq (12) , the only contribution to 1^. can

therefore come from the poles at q = Qr

, where Qr

is subject to

the constraint given by eq (C.35). Hence, we include only these

contributions, so

with a similar expression for I . The total integral of the force

function is given by eq (C.33).

If a Qr

does not exist which satisfies (C.34), then 1^ is

zero and eq (12) is satisfied. If 1^ is nonzero, then additional

force functions have to be applied to the solid corresponding to

the same Qrwhich give a finite contribution in eq (C.38) so that

eq (12) is satisfied. These terms will give zero at x = 0,

since they are the solution of the homogeneous Hilbert problem.

The expression on the RHS of eq (C.38) diverges for y^ = ».

The singularity will cancel from the force integral calculated

from the final solution. The force integral will be zero, since we

have ensured that eq (12) is satisfied by the final solution.

Qr

*L/2

0.5

(C. 38)

82
Free Surface

f

A

(y)

/

//

/

/

/

Solid A (UHP)

Interface

/

/

f

B

(y) /

/

/

Solid B (LHP)

Figure 1-1: A bimaterial composite containing a free surface normal to the

interface and the coordinate system used in these calculations. The force

functions FA

(y) and F

B

(y) are applied just outside the free surface.

83
S-i O

S6 '0- 03 '/ - '<5'- \?- '€- 60'*>-

84
85

LO

+ 1 '0 Z0'0 00 '0 90 '0- CI '0- 02 '0-

87
T3H

i—

I

O

M

O

t-^'SSl 06' 1 1 ZS'Cir88
ooHo

COo

4-1C

CD

4-1Oco

•H

+JoC3

4-1

ctf

co

CN

CU

co co

CO Cfl CU o

CO -i-l

cu

Cfl

5-1

co

cu •• c

—i

M CX

a; cu

fc'- H

C

89
10 '0 *0'0- 0t '0- Si '0- 0Z'0~

90
Imaginary Axis

i

Figure 1-9: The contour

function HA (z) . The cut

angle is 0 on the branch

used for evaluation of the

PQRS is along the positive

PQ and 2tt on the branch RS

integral for the

real axis. The phase

91
Figure I- 10: The contour used for evaluation of the integral for the

function H

B (z) . The cut PQRS is along the negative real axis. The phasf

angle is n on the branch PQ and -tt on the branch RS

.

92
Figure I -11: The semicircular contours used for evaluation of the

q-integrals in the UHP (solid line; for mod z > 1) and in the LHP (dotted

line; for mod z < 1).

93
!

I
Part II

Generalized Plane-Strain Analysis of a Bimaterial Composite

Containing a Free Surface Normal to the Interface

by

V.K.Tewary and R.D. Kriz

95

1. Introduction

In this part, we extend the elastic plane-strain Green's

function calculated in Part I to account for the generalized plane

strain. We apply it to calculate the stress and the displacement

field in a bimaterial composite containing a free surface normal

to the interface which is subjected to generalized plane strain.

The generalized plane strain is induced by subjecting the solid to

an out-of-plane load such that the strain component c 33 is

constant.

The result is obtained in terms of a closed integral

representation, which is evaluated analytically as well as

numerically. As in Part I, the analytical evaluation of the

integral leads to a series representation—not necessarily

convergent—of the stress from which the singularities in the

stress can be identified precisely. The numerical calculation of

the integral representation does not suffer from convergence

problems. It contains singular as well as nonsingular terms and

gives reliable values of the stress and the displacement field in

the entire solid. The method is applied to a cubic solid

containing a Z-5 grain boundary and to fiber-reinforced laminated

composite with different relative fiber orientations.

In general, the nature of singularities in the generalized

plane-strain problem is similar to that obtained in Part I. The

power of a singularity depends upon the elastic constants of the

constituent materials, whereas its weight depends upon the nature

of the loading as well as the elastic constants. In some cases the

singular term may also have a logarithmic oscillatory factor.

The plan of this part is as follows. In section 2, we

reformulate the Green's function method for generalized

plane-strain problems and apply it to calculate the displacement

97
and stresses in an infinite (without a free surface) bimaterial

composite subjected to an out-of-plane load. The displacement and

the stress components t . thus calculated satisfy the required

continuity conditions at the interface. We then use these results

and the technique qiven in Part I to calculate the stress and the

displacement field in a similarly loaded composite which contains

a free surface. These calculations are given in section 3. The

analysis is applied to a Z-5 tilt grain boundary in a cubic solid

in section 4 and to laminated composites in section 5. A summary

and a discussion of the results are presented in section 6.

2. Displacement and the Stress Field in an Infinite Composite

Subjected to Generalized Plane Strain

In this section we calculate the displacement and the stress

field in an infinite bimaterial composite which has no free

surfaces. A schematic representation of the model solid considered

in this section is given in figure II-l, which also shows the

coordinate axes. This model is same as used in Paper I. The

notation and the coordinate axes used in this paper are the same

as in Part I. The notation in this paper is mostly the same as

used in Paper I

.

The displacements in both UHP and LHP have to satisfy the

following equation of elastic equilibrium (add superscripts A or B

corresponding to UHP or LHP respectively)

:

a

2

u.

c ikjl 3_ = o. (1)

dXkSX-j^

Equation (1) is the homogeneous part of the elastic

equilibrium equation valid in the absence of any body forces. As

98
shown in Paper I, its solution in a particular region can be

written using the elastic Green's function as

u(x) = G(x,x')f(x') dx' , (2)

L

J

where f(x') is an arbitrary function which is applied in a region

L such that the region L is outside the region of solution of eq

(1) . The integration in eq (2) is carried out over the entire

region L. The function f (x' ) in eq (2) is usually referred to as

the force function. The force function in eq (2) does not

represent a physical force on the solid. It is a mathematical

artifact which is used to satisfy the prescribed boundary

conditions

.

In the plane-strain problem, u(x), G(x,x') and f(x') are

functions of only the x

1 ~ and x^- coordinates and are independent

of the x.j- or the Z- coordinate. Even in plane-strain problems,

the Z- components of these quantities, in general, will not be

zero. The solution of eq (1) using eq (2) for the plane-strain

case has been discussed in Paper I.

In the generalized plane-strain problem, the solid is

subjected to a load in the Z- direction such that the e 33 is

constant everywhere in the solid and has a prescribed value, e

being the strain tensor. The displacement field in the solid can

therefore be written in the form

Vl^) = ^ * l3 x

3

+ u

p

.(x), (3)

where t) is the prescribed value of e 33 so

H = c 33 (4)

99
is constant and u (x) is a vector which is a function of only x„ -p - J 1

and x

2

and not of x

3

. The subscript « on u indicates that it

denotes the displacement field in an infinite composite with no

free surfaces. Similarly, the subscript p on u on the RHS

indicates that it is a solution of the plane-strain problem as

described below.

Following the notation used in Part I, x, x' will denote a 2D

(2 dimensional) vector on the XY- plane with Cartesian components

x^ f x

2

and x' ^, x'

2 respectively. The 3D vectors will be denoted

by the upper case variables X, X' . For example, the position

vector X in eq (3) denotes a 3D vector with Cartesian components

x, , x0 and x_ . Our object in this section is to determine u (x) 1' 2 3

J -p v -' such that u

m

(x) is a solution of eq (1) and it satisfies the

following boundary conditions,

u*. [x

1 ,0] = u*. [x1 ,0] (-co * Xl s .) (5)

and

T«2i r x

i '°] = T

!2 i

t x

i '°3 ("M * xi * < 6 >

where T

m2 is a vector which is defined in terms of the components

t . of the stress tensor (for superscripts A and B) as

T

»2i = Ti2 • < 7 '

The stress components can be calculated from the displacement

field by using the formulas,

Su, flu. 5u,

T il = C illk — + C il2k — + C il3k — ' < 8 > SX

1

9X

2

aX

3

100
Suk

Suk

8uk

Ti2 = Ci21k ^

+

° i22k + Ci23k 17

3

(9)

and

T

i3 = C i31k — + C i32k — + C i33k — ' 9X

1

ax

2

5X

3

The first term in eq (3) is a constant with respect to the

differential operators in eq (1) . It would therefore automatically

satisfy eq (1). Since u (x) depends only on x and x , we can use

the plane-strain Green's function obtained in Paper I to determine

this part of the displacement field.

First we consider the displacement field in the solid A in

the UHP. Following the method given in Paper I, we apply

hypothetical forces just outside the region A at the set of points

(x 1# -e) where x^j^ varies from -co to + oo and c is a small positive

constant which will be taken to be zero in the limit. The positive

value of e ensures that the hypothetical forces are applied just

outside the region of solution. We denote this force function by

A

f (x

1

) . By using eq (2) and as given in Paper I, we can write

Up(x) in the form,

,00 Hp

A

(£) = GA (x;x^-c)fA

(x^) dx^ , (11)

—CO

where G (x;x') denotes the Green's function for the solid A when

the solid B is not attached to it.

As described in the previous paragraph, the Y component of x'

has been taken to be -e to ensure that the force f (x^) has been

applied just outside the region A. Equation (11) then gives a

solution of eq (1) in the entire space of A.

101
The Green's function for the solid A, in the absence of solid

B, can be written as follows( see, for example, Paper I):

GA(x;x') « - JL £ y

A

(a) ln(zA - p

A

) , (12)

where

za - x

l

+ P« x

2 ' < 13 >

P* " X

i

+ P«

X

2 • < 14 >

A A

and r (a) and p

a

have been defined by eqs (A. 17) and (A. 19)

respectively in Appendix I-A of Part I. From eqs (3) and (12) we

obtain

u

A

(X) = u

Q - -A- £ y

A

(a) ln(zA - t-Le)

f

A

(t) dt, (15)

a -co

J

where it is understood that only the real part of the RHS in

eq (15) represents the displacement field and where

U

0i = 77 5 i3 X3 * (16)

In the arqument of the loqarithmic function on the RHS of eq

A

(15) we have written p e = te. This is because we have to put e =

• A • 0 in the limit and therefore the maqnitude of p^ is not

siqnificant. The only thinq which could possibly be siqnificant in

the loqarithmic term, as in Part I, is the siqn of the imaqinary

part of pa

which, as qiven in Appendix I-A, is positive. However,

in the present case, we do not need to introduce a cut in the

complex plane so that the limit e = 0 can be taken immediately.

102
Following the method given in Papers I and II and also Part

I, and using eq (9), we obtain the following expression for the

stress

£2 <*> - 202 - '4- I 2 <«> J A ~ ' < 17 > a -oo J Z - t+LC

a

where the vector T

Q

is defined as

T

02i " C i233 * • < 18 >

In a similar fashion we obtain the following expressions for

isplac

in the LHP:

the displacement field and the stress component t . 9 in the solid B

00 (X) = u

Q - ][ i

B

(oc) ln(z* - t+te)fB (t)dt, (19)

a -oo

J

f

B (t)dt

=^02 -^rl?w I -b - > < 20 > a -oo

J z - t+LC

a

where

S02i " C!233 * • < 21 >

103
Now we apply the boundary condition given by eq (5) on eqs

(15) and (19) . This gives the following relation between the two

force functions:

Is £A(t) = *s £B(t) * (22)

A B For convenience, we write f (t) and f (t) in terms of another

force function f

Q

(t) as

^(t) = [Ig]"

1

f

Q

(t) (23)

and

f

B

(t) = cig]"

1 f

0 (t) (24)

A B We now substitute for f (t) and f (t) from eqs (23) and (24)

in eqs (17) and (21) and apply the boundary condition given by eq

(6) . After rearranging the terms, we obtain the integral equation

TT -co x

i " t = 5 t

0

, (25)

where

D = -B -B -1 A . A -1a [If ]

-1

(26)

T = T

B

-0 -02 TA

-02 ' (27)

104
and the integral on the LHS of eq (25) is a principal value

integral

.

Equation (25) is in the standard form of a Hilbert transform.

Its solution is given by

-00 dt

, (28)

—00

x

l -

1

where we have used the fact that the RHS of eq (25) is a constant.

The integral on the RHS of eq (28) can be easily carried out. For

this purpose, it would be convenient to evaluate the integral

taking the upper and the lower limits to be t and -t

respectively and then taking the limit t - oo. The result is

f

Q

(x

1

) = i D Tq , (29)

where we have taken the following limiting values of the

logarithmic function:

ln(x. x + t ) = ln(t ) (30) 1 00 x 03 '

x '

and

ln(x, x - t ) = ln( t ) + til (31) 1 oo'

x oo ' x/

in the limit t -» oo.

00

From eqs (23) , (24) and (29) , we obtain the following

solution for the two force functions:

fA(t) = -L M T

Q (32)

105
and

f

B

(t) = (33)

where the matrices M and N have been defined in eqs (A. 11) and

(A. 12) in Part I.

We see from eqs (32) and (33) that the two force functions

are constants. Using these results in eqs (15) and (19) gives the

required expressions for the displacement field in the solid. The

logarithmic integrals in eqs (15) and (19) can be easily

evaluated. Thus we obtain the following result for the

displacement field in the UHP and the LHP:

u

A

(X) = u0 - Y «r A

(a)M Tn zA

-oo v -0 L - v ' - -0 a

a

and

u

B

(X) = un - V y

B (a)N T„ z B

, -oo x -0 L - v

' - -0 a

a

(34)

(35)

where we have used eqs (30) and (31) and, as is usual in continuum

mechanics, ignored the terms which correspond to the rigid body

displacements in the solid.

For later use, the expressions for the stress components

and whicn can be obtained from eqs (34) and (35) by using eqs

(8) and (9) , are given here.

ti - loi - 3

T

o ' (36 >

T

B = TB

- Q

B N T (37}

-ool -01 -s - 0 '

{ '

106
and

^2 " ^02 " 4 S T

0

< 38 >

£2 m

So2 " 5 H T

o ' < 39 >

A B A B where the vectors and are defined in analogy with eqs

(7) and (18), as

£li«> = Tii

B

< 40 >

and

mA , B A , B , . _ . 2dii = c il 33 71

• < 41 >

Equations (34), (35), (38) and (39) satisfy the two boundary

conditions given by eqs (5) and (6) . Moreover, the second term in

both the eqs (34) and (35) is a function of x

1

and and not of

x

3

, whereas the first term is a function of x

3 only. Hence the

value of the strain component e 33 will be 17 as prescribed by eq

(4) . Equations (34) and (35) are, therefore, the desired solution

for the displacement in an infinite composite solid. It may be

noted that x . _ is constant in the solid. 12

3. Displacement and the Stress Field in a Composite Solid

Containing a Free Surface

In this section we shall apply the results obtained in the

last section to calculate the displacement field and the stress

distribution in composite containing a free surface. The model

solid and the coordinate axes are shown in figure 1-1 of Part I.

Mathematically, the difference between the calculations of Part I

and the present calculations arises because in Part I a line load

107
was specified at a given point whereas in the present case the

strain e 33 is specified.

We start by writing the displacement fields in A and B in the

following form, which is analogous to that given in eq (3)

.

u

A

(X) = u

A

.(X) + u

A

.(x) (42) l

v -' ool v - 7 Sl v -' v '

and

u

B

(X) = u

B

.(X) + u

B

.(x) , (43) 1 — col v -' SI v -' ' v/

A B where u ' (X) are the displacement field vectors in an infinite

composite with no free surfaces as given by eqs (25) and (36) , and

A B

u ' (x) are functions of only x. and x_ , and not of x_

.

-s — 12 3

In view of the remark made at the end of the preceding

section, the strain tensor calculated from displacement field as

given by eqs (42) and (43) will satisfy eq (4) . Our task,

A B therefore, is to determine u ' (x) such that eqs (42) and (43) are TM* s

^TM solutions of the elastic equilibrium equations [eq (1) with

appropriate superscripts] in the HHP and LHP respectively and

satisfy the boundary conditions,

u

A

j,[x1 ,0] = u

B

.[x

1 ,0] (-co < Xl * »), (44)

x

A

i2 [x

1 ,0] = z

B

Ll [x1 ,0] (-co < x± * co), (45)

and

T il C °' X

2

] =

° ( ° " X

2 "

M)

' (46)

T

B

±1 [0,x2

] =0 (0 * x

2 ^ -oo). (47)

108
We can determine u ' (x) by using the plane-strain Green's

function as quoted in Appendix I-A. This Green's function is a

solution of eq (1) and satisfies the boundary conditions given by

eqs (44) and (45) . In order to satisfy the boundary condition

given by eqs (46) and (47) , we shall use the same technique as

used in Part I

.

3A. Solution in Terms of an Integral Representation

As in Part I, we apply a hypothetical distribution of line

A

forces just outside the free surface denoted by F (t) in the UHP

and F (t) in the LHP. As shown in figure II-2, these forces are

applied just outside the free surface at continuous set of points

x

1

=-e and x

2 =t, where e is a small positive number which is taken

to be zero in the limit. The range of the real variable t is

(0 * t * oo) for F

A

(t) and (-« < t ^ 0) for F

B

(t) in the LHP.

Thus, using the definition of the Green's function and as

given by eq (2) , we can write the following expressions for

u

A ' B

(x) in UHP and LHP

s -

oo 0

Ug.(x) = I G

Q

AA (x;t)FA

(t)dt + G

Q

AB (x;t)FB (t)dt (48)

0

and

CO 0

u

B

.(x) = ' G

Q

BA (x;t)FA

(t)dt + G

Q

BB (x;t)FB (t)dt / (49)

~00 0

where G

Q

denotes the Green's function for an infinite bimaterial

composite which is given by eqs (A.l) — (A. 4) in Appendix I-A.

Using these equations and eqs (117) and (118) of Part I, we obtain

the following expressions for the displacement field from eqs

(42) , (43) , (48) and (49)

109
u

a

(x) = u^(x) + l I

A

(a)

r Z + C ^ a

zA

+ e

a (X

+ L

-A

z

Sp L

1 Pp J 1 Pp J

a0 l

*AB< a '^ 5

B a

K P

B

$

+ EAB (a # /3) U

B

(50)

and

(for 0 ^ x

2 ^ oo)

u (x)

-B

z = H*(X) + t £

EBA (a^) UA

[ ^] + ( a , /3

)

a/3 L

B

f z~ \ a

K P

+ L a u

B / x TTB

y (a) y

' z a

+ g

- B

+ r

B

(a) u

B

-B

,

Z

a

+ C

-B

- L

a/3 L

EBB (a,P) U

B

-B

Z _. \

a

v P

B

0

+ EBB (a,0) H

B

B

z a

(51)

(for -oo < x

2 ^ 0) ,

where, as in Part I, the U- functions are defined as the

indefinite integrals,

U^z) = HA

(z)dz , (52)

110
U

B

(z) = H

B (z)dz , (53)

HA

(z) =

F

A

(t)

27TL

0

J

and

H

B

(Z) =

t - z

-° F

B

(t)

dt

27TL dt,

(54)

(55) t - z

and the other symbols have been defined in Part I.

For later use, we also write below the stress components

'

' 1' w^^c^ can ke easi]

given in eqs (50) and (51)

.

and t • «i which can be easily obtained from the displacement field

To(x) = TA

2oo

+ I ]T {**(«)/ p

a

} 5

A

a

f

Z " 1 A v P

* +

j

2

A

(a)/ p

A

| H

A

-A

.

Z

cc + c

-A

a/3

DAA (a,<3) H

A

-A

-A

z a

^ P

A

'

- L

a/3 l

DAB (a,P) H

B

A

' P,

B

+ DAB (a,f3) H

B

-A

z a

(56)

111
(x) = T

B

2 00

+ L

+ L

r -B B

L eBA<«^) s

A

( ^-) + eBA(a , P) „M !«_]

I [ p£} s

B

( {?(->/ Pa

B

} b

b( !llLj

Pa/a

- L

r -B B

(57)

ool

«

L Pa ^

1 Pa J

+ i

a/3

.A'

2

A

a

Saa'"'*3

) !T 4~ + K (a,P) HA

f -I-'

112
- L

a|3 L

K

BB («,0) H

B

-B

k P

B

0

+ K

BB (a,/3) H

B

B( fa_ 1 -B 1 P

0

(59)

where the vectors and T

m2 are given by eqs (3 6) -(39) and, in

an analogous manner, we have defined the vectors T. (x) and T_ (x)

th

— such that their i components are equal to, respectively, z ±i^^

and r

i2 (x)

.

As mentioned earlier, eqs (50) and (51) will satisfy eq (1)

everywhere in the region of solution and also the continuity

conditions at the interface as given by eqs (44) and (45) for all

A B values of F (t) and F (t) . This is because the forces have been

applied outside the region of solution and the Green's function

used in writing eqs (50) and (51) satisfies eqs (44) and (45) . As

in Part I, we have to determine these two force functions so that

the free surface boundary conditions as given by eqs (46) and (47)

are also satisfied. In addition, the force functions have to

satisfy the following condition required for the displacements to

be single-valued:

00 (J F

A

(t) dt + F

B

(t) dt = 0 (60)

—00

The H- functions are determined in the same way as in Part I.

We write these in the form,

HA(z) =

—00

TTA, . iq-0.5 , \f (q) z dq

113
and

Y

B (q)z^-°- 5 dq , (62)

—00

where (q) and V

B

(q) are the unknown functions of the real

variable q, which have to be determined.

Following exactly the same procedure as given in Part I and

using the same block representation, . we obtain the following

result for the vectors V^q) and V

B

(q),

H

B

(z) =

M(q) V(q) = N(q) (63)

where M(q) is exactly the same as given by eqs (38) -(41) of Part I

and the block elements of N(q) are

N a (q) =

-

and

N

B

(q> = -

-03 1

2m

-col

2TTL

00 -Lq-0.5 , y

H dy (64)

0

J

-00

I y I * dy

.

(65)

The integrals in eqs (64) and (65) can be evaluated in an

elementary manner. However, in order to be able to handle the

apparent divergence at oo, we evaluate the integral in eq (64)

between the limits 0 and y^ and that in eq (65) between -y^ and 0,

and in the end take the limit y = «. Thus we obtain the following

expressions for the block elements of N:

114
NA

(q) = t£(q)7(-Lq+0.5) (66)

and

N

B

(q) = t®(q)/(-Lq+0.5)

,

(67)

where

and

t

A

(q) =

—oo

-1

t

B

(q) =

— 00

x

-oo 1 -Lq+0.5

^00

2TTL

-oo 1 -Lq+0.5

^00

27TL

(68)

(69)

The transform of a constant in q-space as given by eqs (68)

and (69) is singular in the limit = oo. However, as shown in

Appendix II-A, this singularity cancels out when the inverse

transform is taken.

The solution of the matrix equation (63) can be easily

written as

Y(q) = [M(q)] 1

t

M (q)/(-Lq+0.5)

.

(70)

For later use, we write the following block structure of eqs (63)

and (70) as in Part I:

-AA

M ^ -BA

M

-AB

M

-BB

1 Y*(q)' 1

J , Y

B

(q).

(-Lq+0. 5)

A 1 — oo

.B

•oo

(71)

115
and

I Y

B

(q)J

(-Lq+0.5) D(q)

C v -BA

r a i tA

Sab — 00

c t

B

-BB J < -00 '

(72)

where C is the matrix of cofactors of M.

As remarked in Part I, the representation of V(q) in terms of

the cofactor matrix as given by eq (72) is useful mainly for the

purpose of a qualitative analysis of the nature of H(z) . For

numerical calculations, it is better to use eq (70) , which

involves only the inversion of a 6 x 6 matrix.

Equation (70) gives the particular solution of the Hilbert's

equation. The general solution is obtained by adding the

homogeneous and the particular solutions and is written as

H(z) =

( A }

3p< z > +

r a i

. 3o< z >.

(73)

where the subscript p denotes the particular solution which is

given by eqs (61) and (62) along with eq (70) . The homogeneous

part of the solution in eq (73) , labelled with the subscript 0, is

H

0

(z) =£h(Qr

) z LQr-°- 5

, (74)

where the vector h is a solution of the matrix equation,

M(Q

r

)h(Q

r

) = 0 (75)

116
and Qr

is a root of the determinantal equation,

D(Q

r

) = 0 . (76)

One of the components of the vector h for each Q will be

arbitrary. It has to be chosen so that eq (60) is satisfied. These

constants can also be chosen to satisfy any other boundary

condition on the stress field, such as those which may be required

at the outer surfaces in a finite solid.

3B. Series Representation for the Particular Solution

In this subsection, we shall obtain the series representation

for the particular solution for H(z) . It may be remarked that for

an actual calculation of H(z) , it would be generally more

convenient to evaluate numerically the integral over q in eqs (61)

and (62). The series representation for H(z) may not give

numerically accurate values of H(z) except for limiting values of

mod(z), because the convergence of the series is not guaranteed.

However, the series representation for H(z) serves two useful

purposes: (i) it can be used as a numerical check on the direct

calculation of H(z) obtained by numerical integration in eqs (61)

and (62) , and (ii) it helps in identifying the singularities in

H(z) at mod(z) = 0, which provides a better physical understanding

of the nature of the stress distribution in the solid. Moreover,

as mentioned in section 1, a knowledge of the singularities in

H(z) is essential for a proper numerical modelling of the solid by

using finite element or other such techniques.

First, we shall consider H (z) . We see from eqs (61) and

(72) that the particular solution for H (z) can be written as

117
?

A

(q) [z/yj iq-0.5 dq

(77)

00 (-Lq+0.5) D(q)

where

P

A

(q) = AB (78)

with a similar expression for H (z)

.

We evaluate the integral in eq (77) by following the same

procedure as used in section 3 in Part I. Comparing eq (77) with

eq (C.16) of Part I, we find that the only poles which can

contribute to the integral in eq (77) are those arising from the

zeroes of D(q) and the single pole at q__ = -i/2. The set-I poles

in eq (C.16) of Part I arising from the function E(q) are not

present in eq (77) except for the single pole at q = CL-^*

Moreover, in the two examples discussed in section 4 and section

5, we shall see that D(q_

1

) is also zero. Hence, the pole at q =

is a second-order pole in the lower half of the complex plane.

Since y^ approaches oo in the limit, we note that mod(z/y

w

) in

eq (77) will approach 0. Thus, as indicated in section 3B and

Appendix I-C in Part I, we have to choose the contour in the LHP

for evaluating the integral over q in eq (77) . Thus, by using eqs

(C.34) and (73) of Part I, we obtain the following series

representation for H(z)

:

lQ -0.5+c

r - RH

(z)

,

(79)

r

where Q is a root of eq (7 6) with the constraint,

118
Im (Qr

) < 0 , (80)

and the other quantities in eq (79) have been defined in Part I.

In eq (79) and in what follows in this section, z is expressed in

units of y for notational convenience. In actual numerical 00 calculations y serves as a scale factor.

00

In analogy with eq (103) of Part I, we write eq (79) in the

form,

(k -0.5)

H

p

(z) = 1 ^p (Qr

) [z] r exp[Lg

r

ln(z) ] + R

R

(z) , (81)

where gr

and -k

r

are the real and imaginary parts of Q^.,

Qr

= gr

- *.k

r

(82)

and

Y

p

(Qr

) = - *(Qr ) g(i<Q

r

) • (83)

We see from eq (81) that H(z) will have the same kind of

singularity and the logarithmic oscillations as discussed in

section 3B of Part I. The singularity will exist only if a Qr

exists such that its imaginary part satisfies the constraint that

0 < k < 0.5 . (84) r

Similarly, the logarithmic oscillatory factor will exist only if

the real part of Qris nonzero.

In eq (81) , the first term contains sum over only those

values of Q which are not higher order poles. The contribution of

119
higher order poles is included in Rjj(z) . In particular, the

contribution R^

Q

(z) of the second-order pole at

q_i = Q_i = -l/2 (85)

to Rjj(z) ^as been given below. Following the method leading to eq

(104) in Part I, we obtain

d' (q_ x

)

E' (q_ x

) + ln ( z ) (86)

Equation (86) shows the ln(z) behavior of H(z) in agreement

with some of the earlier work (see, for example, [10] and other

references given there) . Further, if the root Qris degenerate, we

shall get logarithmic terms with higher powers. For example, if

for some Qr

* -l/2 which satisfies the constraint given by eq

(84), D(Q

r

) has a factor (q - Q^.)^ , then by using the Cauchy's

theorem we can easily show that R„(z) will have the form,

Rjj(z) « [z]

(k -0.5) „ ~ 1 [ln(z)]^" exp[Lg

r7?

ln(z) ]

,

(87)

where is written in terms of its real and imaginary parts as

Q = g - t.k (88)

with the constraint

k > 0 . (89)

The subscript tj in eq (88) identifies the particular value of

Qr

as a 7] —fold degenerate root.

120
3C. Force Functions and the Homogeneous Solution

In this subsection we calculate the force functions and their

integral by using the same technique as in Part I. By using the

Plemelj relation as given by eq (B.42) of Part I, we obtain the

following expression for the two force functions:

ip(y) = E(q)PA

(q)dq

(iq-0.5) D(q)

iq-0 .

5

— 00

(90)

and

i?(y)

exp(irq)E(q)PB (q)dq

(iq-0.5) D(q)

'y/y.

iq-0 .

5

(91)

—00

where P (q) has been defined in eq (78) , and

P

B

(q) = (92)

The integrals over q in eqs (89) and (90) can be evaluated by

using the same method as given in Appendix I-C of Part I for

calculating the integral in eq (C.31). Unlike eq (C.31) in Part I,

the integrands in eqs (89) and (90) do not have the factor E(q) in

the denominator but have the factor (-iq + 0.5). The singularity

arising from this factor at q = -t/2 cancels out with the zero in

E(q) in the numerator. The integrands, therefore, have only a

simple pole at q = -l/2 which is a zero of D(q)

.

We need only to calculate the total integral of the force

functions in order to satisfy eq (60) . The total integral of the

121
force functions has also been evaluated in Appendix I-C of Part I.

From eq (C.33) of Part I and eqs (88) and (89), we obtain the

following expression for the integral of the force in UHP:

Y

CO ip(y) dy

,0° Yn E(q) PA

(q)dq

. D(q) (iq+0.5) (-iq+0.5)

— CO

We take the limit as Y > y -> co. As explained in the CO CO Appendix I-C of Part I, only the poles in the LHP contribute to

the q-integral at the lower limit of y, and their contribution is

zero. At the upper limit of y, the q-integral is shown in eq (93)

.

In this case the poles in the UHP will contribute in the limit

Y =00.

CO

First we consider the pole at q = l/2. In contrast to eq

(C.33) of Part I, the integrand in eq (93) has the E(q) term in

the numerator, which cancels the singularity arising from the

factor (-Lq+0.5) in the denominator. However, as given in eq

(C.25) of Part I, D(q) is also zero at q = l/2. Thus the integrand

in eq (93) has a simple pole at q = l/2 which arises from the zero

of D(q). The contribution of this pole is given by

J T0 ** n\ /

D '<V ' < 94 >

where the prime denotes differentiation with respect to q and

q0 = l/2 . (95)

(93)

122
Similarly we obtain

Using eqs (78), (92), (94), (96) and eqs (93) and (94) of Part I,

we see that the contribution of the pole at qQ

to the total force

integral is zero.

Further, as in Part I, we note from eq (93) that the

contribution of a pole for which ImfQ^J > 0.5, will be zero in the

limit = oo, because it would lead to a negative power of Y^.

Thus, we see that only those poles will contribute to the integral

which satisfy the following constraint:

0 < Im(Q

r

) < 0.5 . (97)

This leads to the following for the total force integral:

lQ + 0.5

i = _ x E (Qr

)yw( Y

M A. ) r

:, = 2TTL ) —

J-^ D'' <V < Qr

+ x /4)

r

1

[

P

A

(Qr

) - Lexp(TTQr ) P

B

(Qr

) 1, (98)

where Q is subject to the constraint given by eq (97) . If a Qr

which satisfies this constraint does not exist, then 1^. is 0.

As remarked in Part I, the apparent singularities in 1^ are

not relevant. When 1^. is not zero, we have to include additional

force terms corresponding to the solution of the homogeneous

equation. The H-functions as obtained by these terms are

123
= £ v£(Q

r

)z LQr-°- 5

(99)

and

H*(z) = £ V* (Qr

) z LQr" °' 5

f (100)

where

Qr

*t-/2

3£(QJ = 2TTL

0 (Qr> =(-LQ 4-0.5)0- (Q ) < 101 >

and

b E

B

(Qr W ) -

2l" (-.Qr+0.5)D'(Q

r

)

' < 102 >

If we separate Q into its real and imaginary parts as

Qr0 = ^ro

+ Lkr0 ' < 103 >

then we see that the H—functions will have the following

dependence on z

:

-(k +0.5)

H(z) * (z)

ru exp[cgr0 ln(z)]. (104)

3D. Singularities in H(z) at z = 0

In this subsection we identify and briefly discuss the

singularities in H(z) . The nature and origin of these

singularities are the same as in Part I, and the discussion given

there is also generally applicable to the present case. The only

difference between the plane-strain case as discussed in Part I

and the present generalized plane-strain case arises because in

124
the former case the integrand for H(z) [see eq (54) of Part I] has

E(q) in the denominator, whereas in the present case it has the

factor (-Lq+0.5).

The singularities in H(z) arise from (i) the homogeneous

solution and (ii) the particular solution. These are discussed

below:

(i) Singularities in the Homogeneous Solution

The form of the singularity is given by eq (104) , where k

rQ

satisfies the constraint given by eq (97) . The exponent of q will

have a value between 0 and -1. This term will, in general, have a

logarithmic oscillatory factor unless is zero. This

singularity will be present only if a Qr[a solution of eq (76) ] exists which satisfies the constraint given by eq (97)

.

(ii) Singularities in the Particular Solution

The singularities in the particular solution can be

classified in two categories according to the nature of poles in

the integrand for H(z) as:

(a) Singularities arising from simple poles. These are given by

the first term in eq (81) . The exponent of a singularity in this

category can have a value between 0 and -0.5. These terms will

also have an oscillatory log factor unless Real (Qr

) = 0.

(b) Singularities arising from higher order poles. These are given

by the second term in eq (81) . There is at least one second-order

pole at q = -l/2 arising from the factor (-tq + 0.5) and a root of

D(q) , which gives rise to a logarithmic singularity as given by eq

(86) . This singularity has no oscillatory factor.

125
If D(q) has other degenerate roots, H(z) will have

singularities of the form given by eq (87) . The exponent of z in

this class of singularity may have a value between 0 and -0.5. The

logarithmic singularity may, in principle, have a power of up to

5. There will, in general, be a logarithmic oscillatory factor

unless the real part of the root Qris zero.

The exponents of all the singular terms are independent of

the loading of the solid and depend only upon the elastic

constants and the geometry of the solid. This is because the

exponents of the singularities are characteristics of the matrix

M, which does not depend upon the type of the loading. The loading

only affects the RHS of eq (71) . The weights of the singular terms

in the homogeneous solution are also independent of the loading,

as is apparent from eq (75) . However, the weights of the

singularities arising in the particular solution will depend upon

the nature of the loading.

4. Application to Z-5 Grain Boundary in a Cubic Crystal

In this section, as an example of the application of the

formalism developed in this paper, we shall apply it to calculate

the stress distribution in stainless steel containing a E-5 tilt

grain boundary. We consider the same model which has been

described in section 6 of Part I. We consider the solid subjected

to an out-of-plane load such that e 33 is equal to a constant tj.

First we consider the solid to be infinite, having no free

surfaces. The stress and the displacement field for this solid

have been given in section 2. The expressions for various

parameters such as y and cr have been given in section 6 of Part I.

As in Part I, we choose the units such that c. . = 1 for solid

A in the UHP and tj is taken to be unity. The calculated values of

126
u

l(£)/ T n(-) ' T

12 ^-^ anc* T22^-^ can ^e eas ^^y obtained from eqs

(34)-(39). We notice from these equations that the displacement

field is a linear function of the distance from the origin, and

the stresses are constants.

The calculated values of the components of the vectors T . —ool

and as defined by eqs (36)-(39), and which are required for

eqs (56) -(59), are given below in units of t) [see eq (41)]:

TA

, = (-1.3, 0, 0) and T

B

. = (-1.3, 0, 0) — ool — ool

and

TA

0 = (0, -1.3, 0) and T

B

0 = (0, -1.3, 0). — 03 2 — 00 2

We note from the values given above that T _ is continuous -oo 2

across the interface, having the same values for solids A and B.

This is required by the boundary condition given by eq (6) . In

this particular case, because of high symmetry, is also

continuous and equal in magnitude to T^

.

Now we consider the effect of a free surface in this model

solid using the formalism presented in the preceding section. The

particular solution for H(z) is given by eq (70) and the

homogeneous solution by eq (74) . In order to see whether the

homogeneous solution is to be included or not, we have to solve eq

(76) for its complex roots and see whether it has any roots which

satisfy the constraint given by eq (97) . The numerical method of

calculation of the roots of eq (76) has been described in the next

section. The two roots of eq (76) which have a positive imaginary

part less than 0.5 are

Q _ = 0.0272 + 0.4976 i (105) ^r0

and

127
QrQ

= 0.0118 + 0.4994 i . (106)

The corresponding singular terms in H(z) are given by eqs

(99) and (100) . The z-dependence of these singular terms is

H(Z) * z" 0,9976cos(0.0272 In z) (107)

and

H(z) « z~

0,9994 cos(0.0118 In z) . (108)

We notice from eqs (107) -(108) that these singular terms in the

stress also have a logarithmic oscillatory factor but with a small

magnitude.

The integrals for the particular solution in eqs (61) and

(62) are calculated numerically. In the numerical calculations we

have taken ywto be 10 and the limits of the q—integral from -10

to +10. For testing the convergence we follow the procedure given

in Appendix I

I -A.

The calculated values of the displacement field u

1

(x) have

been shown as function of x.^ and x

2

in a 3-D graph in figure II-2.

We notice from this figure that, as required by the boundary

condition given by eq (44), the u

1 (x) is continuous across the

interface. The behavior of other displacement components is

similar to that of u^.

The calculated values of the stress components "t^fx) and

r 22 (x) have been shown as functions of x.^ and x

2

in 3D graphs in

Figs II-3 and II-4 respectively. We notice from these figures

that, as required by the boundary condition given by eq (45) , these stress components are continuous across the interface. A

similar behavior is shown by the stress component T 22' In

addition, the stress component x

12

is zero at the free surface,

128
which satisfies the boundary condition given by eqs (46) and (47)

.

The stress component t-^Cx) is shown in figure II-5 as a function

of x^ and x

2

. It is also seen to be zero at the free surface,

which satisfies the boundary given by eqs (46) and (47) . It is, of

course, not continuous at the interface. A similar behavior is

shown by the stress component ^ 31 »

As is apparent from Figs II-3 through II-5, the stress

components are singular at the origin. The nature of the

singularities has been discussed in the preceding section. The

exponents of the two singular terms which arise from the

homogeneous solution are determined from the values of the roots

in accordance with the discussion given in section 3D. The

behavior of these singular terms is given by eqs (107) and (108)

.

In order to identify the singularities in the particular

integral, we have to calculate the roots of eq (76) in the lower

half of the complex q-plane. As remarked in the previous section,

eq (76) has one root at q = -i/2. The corresponding singularity in

the H-function will therefore be a log singularity as given by eq

(86) . The other root of eq (76) with negative imaginary parts and

which satisfy the constraint given by eq (84) is

Q = 0.0163 - 0.3342 i. (109) pr v '

The corresponding singular terms in H(z) will be given by eq

(81) . The nature of this singularity is

H(z) « z~°* 1658 cos(0.0163 In z) . (110)

We find that this singularity also has a logarithmic oscillatory

factor with a small coefficient. As remarked in section 3D, the

exponents of these singularities are independent of the loading

129
and are the same as obtained for the plane-strain problem in

section 6 of Part I.

5. Application to Fiber-Reinforced Laminated Composites

In this section we further apply the formalism developed in

this paper to calculate the displacement field and the stress

distribution in a fiber-reinforced laminated composite containing

a free surface and subjected to an out-of-plane load. We shall

consider a planar interface separating two different

fiber-reinforced layers, terminating at a stress-free edge, which

is taken to be normal to the interface. As mentioned in the

introduction, this problem is important for studying the failure

of composite materials.

A real fiber-reinforced laminated composite contains several

layers and interfaces separating them. In general, the interfaces

are not planar. However, we make the usual simplifying assumptions

about the material by considering it as a bimaterial composite

with a single planar interface. We assume that both the layers are

homogeneous and anisotropic, and have no boundaries except for the

interface separating them and the stress-free edge. In spite of

various approximations, the analytical model calculations are very

useful because they provide a good understanding of the underlying

physical principles of failure of composite materials. Moreover,

they provide a basis for an efficient numerical approach for

calculation of the stress distribution in a real material [14].

Thus we are considering a model in which the two solids, A

and B, in figure II-l are the same material except for the

orientation of the fibers. As in the earlier papers, we specify

the fiber orientations in terms of 0, the angle of rotation about

the Y-axis (x

2 ~axis) in figure II-l. The angle 9 is measured from

the Z-axis as shown in figure II-6. The elastic constants of the

130
two solids, A and B, are therefore related by a rotation

transformation given by the orthonormal matrix,

S =

cos 0 0 sin 0 ]

0 10

-sin 0 0 cos 0

(111)

The bimaterial composite in figure II-l is usually specified

in terms of the fiber orientations as (8,0'), which means that the

fibers in solids A and B are at angles 0 and 0', respectively,

from the Z axis as shown in figure II-6. The orientation 0=0

refers to the case when the fibers are parallel to the Z axis.

In this paper we shall consider composites with the

orientations: (0,90), (30,-30) and (45,-45) where the angles are

given in degrees. There is a strong interest in the engineering

properties of composites with these orientations. We shall

consider a typical high-modulus graphite/epoxy material. The set

of engineering elastic moduli for this material for 0=0, as given

in [10], is

l

= E

2 = 1.447; E3 = 13.79; G

12 = = = 0.586;

"21 = "31 = "32 = °- 21

where E. and G. . are the Young's and shear moduli respectively in

1 _ 1 j

units of 10 kPa, and are the Poisson's ratio. The above

constants are based upon the assumption of hexagonal symmetry

(transverse isotropy) and are the same as used in [10]. The

calculated values of the corresponding elastic constants, which we

2 shall refer to as Set I, are given below (in units of MN/m ):

Set I

131
C ll

= C22 = 15 * 25; C

12 = °' 329; C13 = C23 = °' 389;

c

33 = 139.5; c44 = c55 = c66 = 0.586.

In addition to the assumption of transverse isotropy, the

elastic constants of Set I are based upon the following

assumption, referred to as the Pipes-Pagano approximation:

G

12 = G31 and V21 = V32 * ( 112i

The validity of the assumption in eq (112) has been questioned by

Kriz [13], who has obtained a set of elastic constants without

making this assumption. This set, which will be referred to as Set

2

II, is given below in units of MN/m along with the corresponding

. . . 7 set of engineering elastic moduli in units of 10 kPa:

Set II

C ll = C

22 = 1 ' 347; C12 = °' 6747 C13 = C23 = °' 618;

c__ = 14.17; c.. = c,, = 0.579; o ce. = 0.336; 33 44 55 66

E

1 = E2 = 1.002; E3 = 13.79; G±2 = 0.336; G23 = G31 = 0.579;

v21

= 0.491; v 21 = v

32 = 0.306.

In our calculations, we have included both the sets of

elastic constants given above. The calculations are carried out

exactly along the lines described in section 3 and section 4. The

main parameters of the results are given in table 2 . This table

gives the values of pa

, the roots of the Stroh determinant, as

defined by eqs (A. 19) and (A. 20) of Part I, for a = 1, 2 and 3 and

for solids A and B. These roots are purely imaginary and agree

with those calculated in [11].

132
TABLE 2

Results for composite materials with different fiber orientations

Parameters Set I Set II

(e,e')-> (0,90) (30,-30) (45,-45) (0,90) (30,-30) (45,-45)

A

Pi

A

P2

A

P 3

B

Pi

B

P2

B

P3

0. 8940L

1.1185L

1. 0l

0. 6345L

1. 0l

4.7669L

1. 0692L

0 . 8626L

2 . 53621

Same

above

1. 0348L

0.80931

3 .44421

Same

above

0. 96881

1. 0322L

1. 0L

0 . 6820L

1. 3126L

4.75551

1. 0915L

0.9265L

2 . 5306L

Same

above

1. 1686L

0.85461

3.4362L

Same

as

above

TA 1

col 2

3

0. 657

0

0

4.617

0

6.873

5.934

0

5. 301

1. 067

0

0

4 .883

0

6. 948

6. 159

0

5.471

T

B 1

ool 2

3

u . / uy

0

0

4

a

. cnbl / 0

-6.873

o . y j 4

0

-5. 301

n on

0

0

4 . 883

0

-6. 948

6. 159

0

-5.471

TA,B 1

oo 2 2

3

0

0.614

0

0

0.639

0

0

0.613

0

0

1.093

0

0

1.078

0

0

1.102

0

Roots con1:ributing to the heDmogeneous5 solutioi

0. 00008

0. 49900

0

0. 49992

0

-0.46563

0

-0. 49999

I

Re

Im

II

Re

Im

Roots con1

0. 00175

0.49978

-0. 00011

0.49993

:ributina

-0. 00351

0.49522

-0. 00097

0.49905

to the pj

-0. 00588

0.49812

-0. 00019

0.49981

articular

-0. 00181

0. 49981

0

0. 49995

solution

-0. 00008

0.49900

0

0.49992

0

-0.48102

0

-0.49999

Re

Im

II

Re

Im

0

-0.46662

0. 00043

-0.49928

0

-0.48839

-0. 00069

-0.49976

0

-0.47443

-0. 00019

-0.49981

0

-0. 38765

0

-0.49999

133
Further, Table 2 gives T - and the corresponding values of

T _ as defined in eqs (36) - (40) .The value of T , is needed for

calculating the RHS of eqs (78) and (92) . We immediately notice

that both the stress components are constants as shown in section

2. In addition, T

m2 has the same value in UHP and LHP as required

by the continuity condition given by eq (6)

.

Table 2 also gives the permissible values of Q in the

exponent iQ-0.5 of z in the definition of H(z) . These values of Q

are the roots of eq (76) . The values of Q which contribute to the

homogeneous solution given by eqs (99) and (100) are those which

satisfy the constraint given by eq (97). The values of Q which

arise in the particular solution are those which satisfy the

constraint given by eq (84)

.

The roots arising in the particular solution are not

explicitly needed for the calculation of the stress and the

displacement field since, as explained in the text, they are

included in the integral over q in eq (77) . However, they

represent the singularities in H(z) and the stress field at the

origin as represented by eq (81) and discussed in section 3D (ii)

.

Table 2 shows that there are two roots which contribute to

the singularity in H(z) as described in section 3D(ii) . Out of

these, the exponent of the singularity arising from the first root

agrees with that obtained by Zwiers, et al., [10] and Ting and

Chou [11] for the same composite. This root has no real part;

therefore, as described in section 3 D(ii) , it will have no

logarithmic oscillatory factor in agreement with [10] and [11].

The second root is very close to -l/2. As described in the

section 3, it can be shown analytically that eq (76) has a root at

-l/2. The numerical method used in the present calculations seems

to give one more root at about the same value which shows that the

134
root at -i/2 is doubly degenerate. This degeneracy has not been

reported by earlier authors. Since this root of eq (76) is itself

doubly degenerate, the integrand on the RHS of eq (77) will have a

third order pole. As discussed in section 3 D(ii) , therefore, this

2 implies that H(z) has a [ln(z)] singularity at z=0 and not just

ln(z) . The weight of the singularity is given by eq (87)

.

The roots which contribute to the homogeneous solution are

explicitly introduced in the total solution for H(z) , as given by

eqs (99) and (100) . From Table 2 we see that in addition to the

singularity in the particular solution, H(z) will have two more

-1+e singular terms of the type z where e is very small—of the

-3 -5

order of 10 for one root and 10 for the other root. The two

roots are apparently degenerate, which indicates a higher order

pole at i/2. However, as described in section 3, the root at l/2

does not contribute to the stress.

The complex roots of eq (76) , as quoted in Table 2 and in

section 4 for the cubic case, were obtained by the Muller's method

by using a library program on the mainframe computer at NIST.

Equation (76) has roots at ±l/2 which are known analytically.

These were given as known roots in the Muller's method and were

not calculated numerically.

The numerical method, however, gave two apparently new roots

close to this value as shown in Table 2. Since these roots just

satisfy the constraint given by eq (97) , they were included in the

calculation. We could not ascertain with more than an accuracy of

_ g 10 that these are independent roots (possibly degenerate with

the one at l/2) or the same root at l/2 which is being picked up

by the numerical method with a round-off error.

In order to verify the independence of these roots, we first

2 divided the LHS of eq (76) by (z +0.25), which should remove the

135
known roots at l/2 and -i/2. Then we applied the Muller's method

again to obtain the remaining roots. The method still gave the

roots as given in the table, which suggested, within an error of

-5

10 , that these are indeed independent and degenerate roots.

However, the uncertainty is not serious because, as shown in

section 3, the contribution of the root at l/2 to the stress is

zero. The contribution of the roots in Table 2, which are close to

this value, is also quite insignificant.

We calculated the integral over q in eq (77) numerically

between the limits -8 to 8. The convergence of this integral has

been discussed in Appendix II-A. The value of y was taken to be

10. The calculated values of u^ T

22 and T ll ^or (°' 9 °)

composite, using Set I elastic constants, have been shown as

functions of x.^ and x

2

as 3D graphs in figure II-7 through figure

11-10 respectively. The same quantities for (30,-30) and (45,-45)

composites using Set I elastic constants and for (0,90), (30,-30)

and (45,-45) composites using Set II elastic constants have been

shown in figure 11-11 through 11-3 0.

We have given the above results for composites with different

orientations because of a strong engineering interest in these

materials. These curves should be useful in detailed numerical

studies of the stress distribution in real composites.

We see from the various figures that, as required by the

boundary conditions, the stress components t

12 and t22 and the

displacement field are continuous across the interface, whereas,

in general, is not. However, in case of (30/-30) and (45/-45)

composites, x 11 is also continuous across the interface. This is a

consequence of the mirror symmetry across the interface.

Apparently for the same reason, u

1

has an extremum (maximum or

minimum) at the interface along the Y- axis in these materials.

136
Although the displacement field is continuous across the

interface, in case of the (0/90) composite u.^ has a very large

gradient at the interface along the Y-axis. This is a consequence

of the strong elastic anisotropy of the solids in UHP and LHP in

the composite. This can be understood qualitatively as follows.

The stiffness is very large in a direction parallel to the

fibers and very low in a direction perpendicular to the fibers. In

case of the (0/90) composite, the fibers in the UHP and the LHP

are parallel to the Z- and the X-axes respectively. In the UHP u

1

is the displacement component in a direction perpendicular to the

fibers, which is against relatively weak elastic forces. On the

other hand, in the LHP, u.^ is in a direction parallel to the

fibers, which is against relatively strong elastic forces. The two

are therefore very different. The same argument explains the large

difference between the values of in the UHP and LHP in the

(0/90) composite.

The stress components and t.. are zero at the free

surface, as required by the boundary conditions given by eqs (46)

and (47) . All the stresses are singular at the origin, as shown by

the peaks in the figures. These peaks appear to be finite in the

figures, because in actual calculations we calculate the stress

close to the origin but not exactly at the origin for obvious

numerical reasons.

6. Summary of Results and Discussion

In this part we have generalized the elastic plane-strain

Green's function method given in Paper I and used the method of

Part I to calculate the stress distribution and the displacement

field in a bimaterial composite subjected to a generalized

plane-strain loading. The stress components are expressed in terms

137
of a function H(z) , which is given by the solution of a

generalized inhomogeneous Hilbert equation. The main effort is in

the determination of the function H(z) . The final result contains

two parts: a particular solution which is expressed in the form of

a closed integral representation, and a homogeneous solution which

consists of a finite number of terms.

The integral for the particular solution can be evaluated

analytically by using contour integration. This gives the result

in the form of a series, with each term representing the

contribution of the poles in the integrand. The series contain

singular as well as nonsingular terms and may not be a convergent

series. The advantage of the series representation is that it

helps in precisely identifying the singularities in the

stress—their weights (the coefficient of a singular term) as well

as exponents. The earlier work in this field, in general, only

attempted to identify the singular terms and calculate their

exponents. As suggested in [10] and [11]/ the weight of the

singularity has to be obtained by a numerical solution of the

boundary value problem.

In these calculations we have considered only one free

surface in the solid and assumed the solid to extend to infinity

in other directions. The effect of outer boundaries of the solid,

which in real cases would have irregular shapes, is very difficult

to include in an analytical calculation. This can be easily done

by using a numerical technique such as the finite element method.

However, a purely numerical technique becomes inefficient and

unreliable in a region where the solution has singularities and

discontinuities. In such cases, the proper procedure would be to

combine the analytical and the numerical techniques by using

enriched finite element method [14] or taking the exact analytical

solution near the origin and joining it by numerical solution in

138
the outer regions. The formulation given in the present paper

should be particularly suitable for this purpose.

In this paper we have shown that the stress will have

—5 singularities of the form z (0 < 8 < 0.5) and ln(z) (as well as

— 1+5 its higher powers) arising from the particular solution and z arising from the homogeneous solution. The homogeneous solution

will contribute only in certain cases as given in section 3 D (i)

.

In addition, depending upon the material constants, the

singularity may have a logarithmic oscillatory factor which is

similar to that associated with interfacial cracks (see, for

example, Paper II)

.

We did not find the logarithmic oscillatory factor in case of

the laminated composites which is in agreement with the earlier

work on these materials. However, in case of the cubic solids

described in the preceding section, we did find this factor,

although with a small magnitude. Our result on the nature of the

singularities of the type z

5 and ln(z) in laminated composites

agrees with that obtained in some earlier papers [10,11]. However,

the singularities with exponents between -0.5 and -1.0 have not

been identified in the earlier papers.

The existence of the homogeneous solution in our result is

due to a boundary condition which we impose on the solution of the

Hilbert eguation. This condition requires the displacements to be

single-valued. It is obviously a physical requirement which comes

naturally in our method of solution because otherwise the Hilbert

equation will not have a unique solution. This seems to have been

missed in the earlier published work. This may be the reason why

the singularities arising in the homogeneous solution with

exponents between -0.5 and -1.0 have not been identified so far.

It can be easily seen analytically that if the exponent of the

139
singularity in the homogeneous solution is exactly -1, it will not

contribute to the stress. In the case of the laminated composites

which was considered in section 5, the only terms which we found

in the homogeneous solution have exponents very close to -1

(imaginary part of Q close to 0.5). These terms do not make a

significant contribution to the stress and therefore do not

invalidate the results published in the earlier papers on such

laminated composites. Even in the cubic case which was discussed

in section 4, the exponents of these terms are quite close to -1.

Regarding the singularities in the particular solution, we

seem to obtain a degenerate root at Q = -i/2 as shown in Table 2.

As remarked in the preceding section, this implies a singularity

2 ... . . . of the type [In z] . The possibility of such singularities has

been acknowledged in [10] but not actually included in the earlier

calculations. In our method, we do not need to isolate the

singularities in the particular solution which are all

automatically included in the integral representation.

One of the advantages of our method is that we do not have to

calculate the exact roots of eq (76) except for the homogeneous

solution when it exists. In these cases, we have only to search

for the roots in a finite range as given by eq (97) . Moreover,

even when we have to calculate the roots, we need to solve only a

6x6 determinantal equation. In the earlier methods, [10,11] one

has to solve a 12 x 12 determinantal equation. The main

disadvantage of this method is that it cannot account for the

other surfaces and interfaces which exist in a real solid. In real

cases, a combination of a numerical and the present technique, as

suggested earlier, should provide an efficient method for

analyzing the stress distribution.

140
Appendix I

I -A

Transform of a Constant and Some Numerical Considerations

In this appendix, we describe the transform of a constant

which is required for calculating the RHS of eq (63) . In this

context we shall also discuss some aspects of the numerical

convergence which are relevant to the calculations in this paper.

Let C denote a constant, which is expressed in terms of its

transform c(q) as

00 C =

—or

c(q)y Lq °* 5dq . (A.l)

Following the method given in Appendix I-B of Part I, we obtain

the following inverse relation for c(q)

C y" 1-^0 - 5

, ( A . 2)

0

J

which is an equation of the same form as eq (64) or eq (65)

.

We evaluate the integral in eq (A.l) in the limits 0 and y00 and take the limit as y goes to infinity. The integral is then

00 elementary, and we obtain the following expression for c(q)

* x CI , ,-Lq+0.5 -i

C «3) = " JWZ (q+0 .5L )

* (A * 3)

In order to verify that this gives the correct transform of the

constant C, we substitute for c(q) from eq (A. 3) in eq (A.l).

Obviously the inverse relation is valid only if the following

141
integral is independent of y in the limit y -» oo and is equal to

unity; that is,

00

27TL

-00' (q+0.5L) w(YjY) oo Lq+ °' 5 dq = 1 . (A. 4)

The integral on the LHS of eq (A. 4) can be easily evaluated

by taking a semicircular contour in the lower half plane as shown

in figure 1-4 of Part I. In the limit of large Ymthe contribution

of the semicircle will vanish, and the integral can be seen to be

unity as required by eq (A. 4) for all values of y.

Equation (A. 4) can be used for testing the convergence of the

numerical integration over q in eqs (61) and (62) . In order to

carry out the numerical integration for H(z) , we have to choose

finite but large values of y^ and the limits of the integration

for q. We also evaluate numerically the integral in eq (A. 4) with

the same values and check whether the integral is acceptably close

to unity.

One problem which arises in this context is that if the

limits of integration for q are taken to be too large, the

integrand in eq (61) may become numerically unstable. Some matrix

elements of the matrix M(q) contain factors like exp(±7rq), which

will become too large near the upper or the lower limit of q in

the integration. The matrix will then become ill-conditioned and

its calculated inverse may contain large numerical errors. On the

other hand if the limits of q are taken to be too small, the

integral may not converge in the sense required by eq (A. 4).

The stress T

1 (x) as calculated from eqs (58) and (59) will be

zero at x

1 = 0 only if eq (A. 4) is satisfied. It is therefore

necessary to choose optimum values for the limits of q so that eq

142
(A. 4) is satisfied without making the matrix M(q) ill-conditioned.

To make the matters worse, the convergence also depends upon

values of y in eq (A. 4). This would suggest that different limits

for the q-integral have to be chosen for different values of y.

This is not advisable because it might introduce artificial

discontinuities in the result.

One way out of this difficulty, which we have used in our

calculations, is the following. We choose reasonable values for y^

and the limits of q and evaluate the integral over q for H(z) in

eqs (61) and (62) using these constants. Then we calculate the

integral in eq (A. 4) by using the same values of y, y , limits of

q and the integration interval. The calculated value of the

integral in eq (A. 4) is taken to be the effective value of unity.

We then divide H(z) by this effective value of unity. The errors

due to lack of convergence are then largely cancelled out from the

result. This is indicated by the fact that T^x) as calculated

from eqs (58) and (59) will become zero at X- = 0 in accordance

with the free-surface boundary condition.

143

f

B

(x) //////

Solid A (UHP)

//////

fw \\\\\\° www

Solid A (LHP)

Figure II-l: An infinite bimaterial composite containing a plane interface and the coordinate system used in these calculations. The force functions

f (x) and f (x) are applied just outside the UHP and LHP respectively.

145
u

o

c

CM

i

i—

i

i—O

i—

Ll

26 '2 90 '2 02"I CC'0 CS'0-

CA 'X )L n

146
o

n

i—

i

o

i—

Ll

8

o

CO

o

4-1

n3

<4-l o

c

o

•H

4-> O

I

cOW CD

cO

^—' cd

CN S-l rH 4J

t- w

W cd w C

a> co

5-l <-h

J-> P- w

CD CDN H »H

r1

CO •• S-l <d

• C

m a; i—i W)

a> a)

hi

•H C

fa •--(

90 '0 *0'0 20 '0 00 '0- 30 '0-

a 'x

147
o

CO

o

i—

i

o

cOU

o

4-1

nj

4-1oS

o

•H

4JO

C

tH

Cfl

CO <u

cfl CO ^ o

- C

X tH ^—< cfl M !-l CN +J K CO

CO CD

co pj

CD Cfl

CO

CD CD H "H•H

cfl •• U

<t CD

• c M CD M M

(!) CD

•H C

fa »H

0£ m0 6t'0 60 '0 10^0- il-0-

148
LO

O

i—

i

Ll

T3

•r-l o

coo

•H

•9o

rrt CO

O/ 1 i

•0aa

CMO

c

o

•Ho

s

CoW CD

co

CO

CJ

•rHCO

tH,

4—* M

CO W CD

S-i ill rrt p-H

1 1

to T3

cu x: N

H •H

r-l

CO m CU

i c M CD M hO

0) CU H rC

4-> bO c

•rH

TTTTTTJ-r7-TTTTJ-r7-rTTTJ

ii '0 90 '0 10 '0

TTTTrTTTTTT-r

fi>0 '0- 60 '0-

149
Figure II -6: The rotation angle of the fibers with respect to the Z or the

X axis. The rotation is about the Y axis which is normal to the plane of

the paper. At 0 orientation, the fibers are parallel to the Z axis. At

orientation 9, the fibers are parallel to OZ'

.

150
O

i—i

£9 '0 * I > > ) ) ) I > i ) i > i ) ) 1 1 )

61 '0 >0'0- 9Z'0- 6*

m0-

'xn

n

u

o

o «H

*4—I 4J

CO

Ctj >—I

T3 <D C

Cfl <—I X 4-1 CD

<4-l CO o c

•I-l w

0)

to

o

•H

4-> Oc

4-1 cfl o

CO

c

10 «H

ccj ctj

5-1 4-> >-> w

CD 3 ' 1 0 CD •H N

4-1 «rH —I

4-> cfl £ !-l CD CD S C

CD CD O bO

ca •—I CD

CO -u

X)

c

•H

CD

CU

H 4->

•--( to .. o ' E • MOW M O 4-> c

cu ^-n nj ^ o u 3 ON 10

•H

O

O

151
00

I

o

I—

I

Ll

o

\o^ o

•rH

flj 4-> CO

5-1 ClJ O •—

1

U-t CD

>*, «-H

i 4-> £ <D

bO

Pi

<4-4 «H o w

c

o cd •H CO

4-> nj o o

C

<4-l -i-l

n3 5-i 4-> CO CO

CD

x a.

*»—

'

CM T)

CO

k- N

•rH 10 i—

I

W Ctf 0

CD 5-1

5-1 CD

4-1 C

t

M CDbO

CD ,C 0) H ^4->

•• c

00 «H

> > > 1 1 > I

) ) ) I )

1 1 ) } 1 1 ) } I

) ) ) I ) } 1 1 ) t 32 '0 02 '0 Ct'0 S0'0 20 m0-

I—I CD CO M 4J 4-> •H C

a) w cd o +j &, CO ME C

o o

u o

5-1 3

152
O

h—

f

Ll

o

\o

O

H

P

CO

S-i cfl o i-H 4-1 <D

U

c CU CO

txO

4-1 •r-( O CO C

o

•H co

4-> O o c

C

4-1 • H

cd »

co CO CO

CUC

cC —1

&,

cs Q)

t- N

•t-l co <—

1

CO cd

CD

J-l 0)C

CO CUCxO <D

<U H

c

<n •r-l

i (—

i

0) M 4J

•i-l CU CO

5-1 o

Qu

W) E

•t-l O

ft. CJ

1 1 f ) /

1 1 / / 1 1 /1 J / / / / / 1 f

I / f > 1 1 J 1 1 / 1 / / CB'0 S9'0 l+'0 02 '0 10 '0-

153
I

o

I—

I

Ll

1 %

%/1

^—* <ys \

(~\ \—* o

ni •HU

rH W /svj td —

1

(D H

a 1 •U

rrt <D

t/3

fi i H

W

M P

i CD

i i CO

r \ \J to r* (J ^

ii

j C

•H

crt 03u

W U

03 W

a;C

03

rS r-l

<uN

w •H

co 1 <u 03

!-l 5-1 4-> Qi W C

0)

(1) W)

H 0) ,C

4J O (3

I •H

i M <D w M 4-) 4J

•iH c

0> W 0) O 4J (X M

w> E c

•-i o o o 0

1 1 ) i 1 1 1 n 1 1 1 1 1 i 1 1 1 1 1 1 ) 1 1 i 1 1

1 1 1 1 1 1

1

91 '0 i»t '0- £> '0- CZ '0- S0 • i - SC 'I a 'x n

154
o

I—

I

Ll

1 1 1 1 j J

I I } I i I J ) I i i DJ

) ) ) t 1 1 J i I ) t f ) 7 62 '0 Z*'0- 21 'I- CB'la 'X ) L n

o

5-1 »H O 4->

Q>

4->

4-1o

c

o

•Ho

c

<4-l

c

•I—

I

w

<uw

o

•rH

cfl S-i

nj ^ —

i

—I 0)

0) N

•H >H

4-> !-i C 0) a) C

e <u o

-I 42 w

0 a) H T-lw

o

a.6

»—

1

Oo

1—

1

l-l 4-> o C

a) cD

S-l 1 4J 3 \ W

W) O C

•H CO o

fn o

155
o

CM

i—

i

O

i—

Ll

} l I l J I 1 I f l J 1 t l I l JV I I I t J t l 1 I

! J ) l ! I I 29 '0 CC '0 £0 '0 92 '0- S3 '0- SB '0-

o

CO o

CO M

S-l W

o m

O

O

CO

bOC

•HW

3

C

CO

oC

o

•r-t O

U to o to C o

d

•H

to

5-1 tO 4-> tO to

>> c

- to X M ^ ex

cm

b- CDN

tO M

tO M

O tO

5-1 S-l •U (U

to C

O

O bO

H CD

1-1

CM C

MOW M -U 4-> •H C

o w to

5-1 O -U 3 ft M

HOE C

•H O O

156
o

COo

CO

o

cO .HU

J-l

w

om

4-1 i-ha;i • C W

cc a; oo X

4-1Oco

txO c

•l-l w3 c

4-1

O

cfl CO

CtJ

W

/-n a; >^ C - cfl M

CN Tj

t- <uN W «r-l W —I ^ S-i p cu w C

CD oo u

t-l >H i • HUB H 4J <J •H C

CD CO cfl

5-1

O -P

d cx

w ME C

•H

O

O

tn

O

U

6< *X)Z3*

157
o

I—

I

CO

O

v - •HU

CO

ca i—

i

r\<j a)

q_|

cu

r!

txO rN c

•H

CO O

C CD o CO •H cd

4-> o

o

c c

•H

ca

ca 4J

CO

10

CO CDc

caa

tj

rH 0)

iH N

•H

rH

co ca 03 u

CD CD U G

4J CD

CO bO

<u

CD H rC

<t c

rH •H

1 M CD w

1—

1

4-) •p •H c

CD CO u O 4J o- co W> E C

•H o o

fa o 0

7 / 7-7 I I J fi t ) ) > J

}) ) ) , ) I

, J }i ; jj , , f } , f Si-0 IZ'0- 93'1- BS'C158
LDOi—i

L.

99 '0- ZClo

S-i •>-< O 4J

4-1W

C

cd • AJ

CO

moco

•Ho

4-1

bO

r-l

CDw

cOo

ca

C

•-I to cfl co ^P

/-s to X C

^_-/ c0 rH

T3 T3

,—1 QJ CD N

03U

4-1

4->E

CD bO O

03 CD 1—I ^

Pu -P W

•r-l

C

CD CD

10O

P-6oo • M W M ^ J-> m C

>-l 4-> 3 \ W

W)

m

C

•H <}•O

fa ^ o

a 'x ) i n

159
'A

I—

I

o

I_L

\U 1 <t

r J

crt »H

jj

CO

M-l

tiC 40

0)

C/3 rS uuhri

4-1 r| o

w

r! -j

\j *H

ij 4-^ rnWi

crt

«-t 0-1 (-1 >rj

crt rrt ycH W 4J

crt W

<De

cd X

M

*H T)

t» <UN

w •H w •—

1

0) cS

)-l 5-( QJ tn c

0)b0 H 0)

4J

M 4-> 4-> •H C

(U M tt)

i-l O 4J M E C

•rH O O

} 1 1 1 1

1 1 1 1 1 1 1) j ) i in ) }

1 1 > 1 1 / }I } 1 1 1 1

)

l* '0 >l '0 CI '0- 0* '0- 49 '0-

160
7TTTTTJTTTTTT1

9fS- 0Z'3- CS'C161
GO

i

O

i—i

Ll

m

i \m

o

4-> U w

O (3

<4-l •—!

CO

>->

-0 •

cfl CO CO

O i-l

o

•-1 CD •U w

s °

i-ln

CO 4->

rt) to ^ s

X •—

i

1-»

fc- CON

CO t-i

co —

i

CD <fl

S-l 5-1

4-1 co

co C

0)

CD bO

H CO

*->

oo C

M O co

4-> i-l c

CD CO U O 4-> 3 a. co oo e C

•H O o

fa O o

) 1 1 J 1 1 / 1 1 1 1 1

1 1 } 1 1 1 1 ) 1 1 1 ) 1 1 } 1 1 1 1 i i I 6S'0- te'l- 0C '€- 99'*- t0'9~

162
cO

% 1

o

) ) ) I ) ) I

) } 1 1 1 ) I} }1 1 1

1 1 1 i ) ) } 1 1 } } ) } }

)

69-0 10 '0- IZ'0~ t*'l- IfZu

o oU

>^ M

cO X) —I C <u cO CN X ' 4J

4-1 Q) O tSi C bO O C

•H «H

4-1 05 c

3

U-t Cfl tO

cO O

cfl C

cO -i-i

cfl

4-) - Cfl

(J) - G

3 cO —

i

—I

CD X)

CD CW N

4-> C

CD6

CDO

CO bO —

I

a* a)

Cfl ,c

•H 4->

a) H

Cfl o

' E • MOW M O 4-> c

CU ^-n cO

5-i O 4->

•h o o

tin ^' CJ

163
CM

i

O

I—

I

ONV

O

trt •HU

l_i CO

vj cd 0 I

CD

CM

*d i f-M •U

crt 0) CO

li i HHw

i

rft WJ

rj r-t P r\\J

/1\

(L*

r

)

u

'

rrt

f l_i C

•H

»u

(/] •u

cd to

(U >> C—

1

CN

i—

1

X)

^, 0)N

CO •H

CO >—

1

CD cS u M

p 0)

CO c

CO

CD bO H CO

4J O a

CNJ

1 •H M CO to M u

H C

cu to u O •p a to E C

•H o o

fa o

i 1 1 ) JI ) 1 1

1 1 1 1 1 1

1

82"I C9'0 10'0~ 99'0~ 0C'l164
CM

i

O

i—

i

Ll

o

\o

O

co •HP

n CO o cO ^

0)

CM c P

co cu CO

bO

tw c

o

CO c

o

•H CD •P co O cfl c O

C

•H

cd cO

}-l to P

cO co

cuc

cO

>—

/

fc- 0)N

01

co ^

cu cO P. U P CU

co

cu CU bO H CU rCP

<—

i

c

CM •H M CU co M P P

•i-l c

<U CO CO

V-i O 4-> 3 & CO W> B c

•H o o

<j o

f 1 I i I

I I I I ) I

} I I t 1 1 I } I I 1 1 } I I I ! I

t i I ) I 96 'I S9'l 92 • I 16 '0 99 '0 61 '0

165
o

ON

is-*- &

.

s. ^ .c. ^.g_ si !jr^Jrg.

c

•HW

o

CO «H

4->

$-1 CO O ct)

<+-l H

CD >•, CM

T3 i C w

cO <D 00

<wO

c

o

•H (J) •P W

O cfl §

°

•H

cO cfl

cO co

^-s CDc

- co X •—

i

ex

iH K 4)N

CO -H

CO

<U co

CDC

CD

CD bO

CD o

CM C

CM «H

t H CD CO t) iJ •H C

CD CO CO

S-l O ±J 3 (X CO W) E C

•H o O

fci o o

166
o

U -H o v

4-1 U5

ctj '

on

4-1o

n

CMi Oi—i

Ll

co

•H

•POC

4-1

(A 'XH

n - 9*> 2~ C2'Cc

lt-1 tno

•H

w to

CO )-i 4-> 5 §

.—i 3 p-

-I <D

CD

N

4-1 r-l V u C o)

0)

C 6 a>

CD b0 o

n) a) —i ^

Pu 4->

»ri

C

t3 >H

(1) CDo

..

S

CM

O i O • M W M ^ •»-> O C

5-1 i 4-> 3 \ W

6DO

C

COo

167
C\J

i

o

I—

I

1 1 1 1 / 1 1 1 1 1 / 1 1 1 1 tyi 1 i i i j. i i i / / / 1 1 i i

1

1 21 '0 90 '0 10 '0 S0'0~ 11 '0- 91 "0-

o

rO

i \o

a") ^-^ o

co •i-l 4J

5-i w o cO I—

1

<D

CM

i 4_) cO cu CO X

o

o

•H 0) 4J w o CO c o

'•1-1 •H

co cOu

w 4-1 cO w

^ N CDb

cO .—

1

M

T>

k- CDN

W •H w ^

CD CO

5-( >-i P CD w b

CD

0) U)

x:H CD

CM

CDU

B

CD W

4-> 4J

.--i B

w <0 o 4-> p« w

E B

o o

o o

168
o

o

CO o

cO •H

4J U 03 O CtJ 4-1 .—

1

CM

i C 4J

cd CD

C/}

faO 4-1o •H

03 c

o

•H d)

4-> 01 O CTJ C O

t+-i C

•H

cd cdH

w 4-1 cd 03

CDC

cd f—

1

v—s p- M

<N

CDN

CO •r-l 03 i—

4

0) cd

S-i 4-> CD

03 c

CD H <u

4-»

tn c

CM •1-1

i M CD 03 1—

1

4-1 4-> •H c

<D 03 cd

S-i O 4-) CL. bO 6 C

•H o o o o

1 1 1 1 1 1 f

1 1 1 1 1 f I ZS ' I ~ 96 • I- 3C m2-

TTTTJ~T7~r7~rT~TTTTTTT

169
\0

CMi

CD

i—i

Ll

o

CO1 \o

CO

o

4J

V-l t/3 O nj

CO

CM

t3 i c w

03 CO * .

<4-lC

O "Hw o

*r-l CO

4->

W

O 0) c o

C|_|

c

•H

nj u) -v CO

>->

C - cd ^ &,

k. CON

CA >i-l W .—1

CO tfl P S-l •U CO W C

CO

CO W) H co p

CM1 M CO

W

M

U 4J •H c

CO t/1 P. 0 P 3 a- to &o e e

•H

O 0

u+ 0 0 1 1 } / 1 } 1 1 1 1 i } 1 1 1 1 1 1 1 1 1 ? ? I } 1 > J I } }

170
o

5-t «H O 4-> IW CO

CO

CD

P CM

ct) '

4-1O

(Mi O

Ll

o

•HUO

4-1

C

cn

COoP

9i 'Z- - I •t~1

t/1 c3

CO r-t 1 \ £— 5 

—1 5 p- ^ <D

<D

N

•HcwP Q>

<u

P

S <u

<D bO O

OS

O ^ .C

P- ^

in •H

P

Q) CD ^ *i H 05O

re

CM

O i o

0) u

inP

cO

5-t ' 4-> 3 > CO bO lO C

•rl <!

o

171 I
in

00

CM

i

i—

i

o

h-H

Ll

TTTTTTTTTT1 ) TTT7 T7 J

} //////// I f / / / I

21 '0 P-0'0 >0'0~ ll'0- 61 '0~

txO c

•H

003

S* o

Cd >H

•U

5-1 CO o cd tW r—

I

(1) >^

CM

-O i C 4J

CE) 0) CO X

m

o

c

o

•T-l 0) u co o cd C o ^ C

•H

nj cd

co ±j

cd co

C

- cd

CM-Xl

t- 0)N

00 t-I

00 —

I

<u cd

to C

0)

co C

CM >H

I • H dl « M 4J -U

a; co cd

>-l O -P 3 ex oo M E C

•H O O

172
) 1 1 1 1

1 J 1 1 J 1 J I P-C't- 99' i1 1 1 r t f r 1 1 1 j ) 1 1 1 1 1 1 1 r i

IT) <f

i \

cd •H

4JW

Q

4-1

CN) TiM 4J

bO

til r

- »-i

'—J w

o

-J f\\ •U w

o

c O

<4-l c

•I-l

w 4-> to

(Dc

CSM Tj

<UN w •H w —1

<u « U (1) w C

0)

<D bO

43H a)

4Jc

CM

i M a> M 4J

•r-( 0) W

c

3

o -u ME C

•-I O o

173
n

i

o

Ll

7

90 '0

TTTTTTJ~TTTTTTJTTTTTT'J7H~rTTTi 23 '0 '0 01 '0- C3'0~

m

\in

<J- o

cO •HU

U 05 O cO r-l

CD >•> CM

i M 4J

co CD CO

txO cw C o •H

03 c

o

•1-1 (1) +J 0) o CO £

O

Cw C

•r-l co CO

CO U

CO

CD >^> C

cO X r-|P1-1

r-t T)

0)N

w •r-l CO i—l

<d COu

CD CO C

CD

<D W> H CD rC O c

co •H

i M CD CO M •P 4J •H C

CD CO CO U O P

ex CO W) 6 C

•r-l O o

rH o o

174
Acknowledgments

The authors thank Dr. H.I. McHenry for his constant encouragement

and support. The authors are also grateful to Drs. Robb Thomson

and Ed Fuller, Jr. , for reading the manuscript and many important

suggestions. This work was carried out in part under a research

contract between the Ohio State University and the National

Institute of Standards and Technology.

References (For parts I and II)

[1] R.B. Pipes and N.J. Pagano, J. Composite Materials, 4, 538

(1970)

.

[2] P.W. Hsu and C.T. Herakovich, J. Composite Materials, 1JL, 422

(1977) .

[3] A.S.D. Wang and F.W. Crossman, J. Composite Materials, 11, 92

(1977) .

[4] N.J. Pagano, Int. J. Solids and Structures, 14., 385 (1978).

[5] S.S. Wang and I. Choi, "Modern Developments in Composite

Materials and Structure," (ASME Publication; Editor J.L. Vinson),

p-315 (1979)

.

[6] R.L. Spilker and S.C. Chou, J. Composite Materials, 14, 2

(1980)

.

[7] S.G. Lekhnitskii, "Theory of Elasticity of an Anisotropic

Elastic Body," (translated by P. Fern), Holden-Day, San Francisco

(1963)

.

[8] Dempsey, J. P. and G.B. Sinclair, J. Elasticity, 9, 373

(1979) .

[9] T.C.T. Ting and S.C. Chou, Int. J. Solids and Structures, 17,

1057 (1981).

[10] R.I. Zwiers, T.C.T. Ting, and R.L. Spilker, J. Appl.

Mechanics, 49, 561 (1982).

175
[11] T.C.T. Ting and S.C. Chou, "Fracture of Composite

Materials"; (Proceedings of the II U. S . A. -U. S . S .R. Symposium,

Lehigh University, Bethlehen, PA, March 9-12, 1981; Edited by G.C.

Sih and V.P. Tamuzs) ; p265 (1982).

[12] N.J. Pagano and R.B. Pipes, J. Comp. Mat.; 7, 679 (1973).

[13] R.D. Kriz, "Effect of Material Properties on Interlaminar

Stresses in Angle-Ply Composite Laminates"; Report No. VPI-77-16,

Virginia Polytechnic Institute and State Univirsity, Blacksburg,

VA (1977) . (A copy of this report can be obtained from the

authors .

)

[14] H.K. Stolarski and M.Y.M. Chiang, Int. J. Solids Structures;

25, 75 (1989)

.

[15] C.E. Bakis and W.W. Stinchcomb; Composite Materials: Fatigue

and Fracture; ASTM STP-907 (Editor H.T. Hahn; American Society for

Testing and Materials) ; p 314 (1986)

.

[16] A.H. England, J. Appl . Mechanics, 32, 400 (1965).

[17] J.R. Rice and G.C. Sih, J. Appl. Mechanics, 3_2_, 418 (1965).

[18] V.K. Tewary, R.H. Wagoner and J. P. Hirth, J. Materials

Research, 4, 113 (1989)

.

[19] V.K. Tewary, R.H. Wagoner and J. P. Hirth, J. Materials

Research, 4, 124 (1989)

.

[20] N.I. Muskhelishvili, "Singular Integral Equations,"

(Noordhoff , Groningen) (1977)

.

[21] N.P. Vekua, "Systems of Singular Integral Equations,"

(Noordhoff, Groningen) (1977)

.

[22] A.N. Stroh; J. Math. Phys. ; 41, 77 (1962).

[23] J.E. Sinclair and J. P. Hirth, J. Phys. F (Metal Physics), 5,

236 (1975)

.

[24] W. Bollman, "Crystal Defects and Crystalline Interfaces,"

(Springer- Verlag, New York) (1970)

.

[25] V.K. Tewary, R.M. Thomson and E.R. Fuller, Jr, J. Materials

Research, 4, 309 (1989)

.

176
BL-114A

(5-90)

U.S. DEPARTMENT OF COMMERCE

NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY

BIBLIOGRAPHIC DATA SHEET

1. PUBLICATION OR REPORT NUMBER

NIST SP-802

2. PERFORMING ORGANIZATION REPORT NUMBER

B90-0273

3. PUBUCATION DATE

March 1991

TITLE AND SUBTITLE

Effect of a Free Surface on Stress Distribution in a Bimaterial Composite

5. AUTHOR(S)

V.K. Tewary, R.D. Kriz

6. PERFORMING ORGANIZATION (IF JOINT OR OTHER THAN NIST, SEE INSTRUCTIONS)

U.S. DEPARTMENT OF COMMERCE

NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY

BOULDER, COLORADO 80303-3328

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (STREET, CITY, STATE, ZIP)

Same as item #6

7. CONTRACT/GRANT NUMBER

8. TYPE OF REPORT AND PERIOD COVERED Final

10. SUPPLEMENTARY NOTES

11 UTERATURE

( SU^EY°MENt"o^?T HEReV

AL SUMMARY OF M0ST SIGNIFICANT INFORMATION. IF DOCUMENT INCLUDES A SIGNIFICANT BIBLIOGRAPHY OR

This paper is in two parts. In part I, the elastic plane-strain Green's function is calculated for an anisotropic

bimaterial composite solid containing a free surface normal to the interface. An exact integral representation is

obtained for the Green's function, which is evaluated numerically. The integral is also evaluated analytically,

which gives a series representation for the Green's function. For illustration, the formalism is applied to the

antipiane-strain problem and the plane-strain problem in a cubic solid containing a E-5 grain boundary. In part

II, the Green's function is extended to the case of generalized plane strain and applied to calculate the stress

field in a cubic solid containing a F-5 boundary and to fiber-reinforced composites with various fiber orientations

subjected to an out-of-plane load.

It is found that, as predicted by earlier authors, the stress is singular at the intersection of the free surface

sad the interface in both plane strain and generalized plane-strain case. These singularities in the stress field

associated with the presence o f the fret surface are identified and discussed. The singularies in the stress field

are found to be of the type r~& , ln(r) as well as containing higher powers of ln(r), where 6 is between 0 and 1 and

r is the radial distance from the intersection of the free surface and the interface. In addition, it is found

that, in general, the stress field will also contain an oscillatory factor of the type exp [ (.gln( r) ] , where g is a constant that depends upon the material parameters cf the two solids. However, for the case of the generalized

plane strain in fiber-reinforced composites, the factor g is almost 0 so that the oscillatory behavior is absent.

Part I also gives a simple and convenient method for solving an inhomogeneous generalized vector Hilbert problem,

which contains a nonsingular kernel in addition to the usual singular Hilbert kernel. The solution of this Hilbert

problem is required for obtaining the Green's function in the present case, as well as many other cases dealing with

the stress analysis of composite materials containing surfaces, interfaces, and cracks.

The solution of the Hilbert problem is obtained by using a complex transform of the type y ^ where y and q are

variables on the real axis. It is shown that this function is an eigenfunction of the Hilbert operator and is orthogonal for y, varying in the range between 0 and * oo and q between -<» to + <».

12. KEY WORDS (6 TO 12 ENTRIES; ALPHABETICAL ORDER; CAPITAUZE ONLY PROPER NAMES; AND SEPARATE KEY WORDS BY SEMICOLONS)

antiplane strain; bimaterial composites; fiber reinforced composites; free edge effect;

free surface; generalized plane strain; Green's function; Hilbert problem; plane strain;

singularities in stress.

13. AVAILABILITY

UNLIMITED

FOR OFFICIAL DISTRIBUTION. DO NOT RELEASE TO NATIONAL TECHNICAL INFORMATION SERVICE (NTIS).

X

ORDER FROM SUPERINTENDENT OF DOCUMENTS, U.S. GOVERNMENT PRINTING OFFICE,

WASHINGTON, DC 20402. ORDER FROM NATIONAL TECHNICAL INFORMATION SERVICE (NTIS), SPRINGFIELD, VA 22161.

14. NUMBER OF PRINTED PAGES

185

15. PRICE

ELECTRONIC FORM ft U.S. GOVERNMENT HUNTING OFFICE: 1990 - 2 8 1 -5 57 ' "t 0

l' * 3



1 t JL k3 M. Technical Publications

Periodical

Journal of Research of the National Institute of Standards and Technology—Reports NIST research

and development in those disciplines of the physical and engineering sciences in which the Institute

is active. These include physics, chemistry, engineering, mathematics, and computer sciences.

Papers cover a broad range of subjects, with major emphasis on measurement methodology and

the basic technology underlying standardization. Also included from time to time are survey articles on topics closely related to the Institute's technical and scientific programs. Issued six times a year.

Nonperiodicals

Monographs—Major contributions to the technical literature on various subjects related to the

Institute's scientific and technical activities.

Handbooks—Recommended codes of engineering and industrial practice (including safety codes) developed in cooperation with interested industries, professional organizations, and regulatory bodies.

Special Publications—Include proceedings of conferences sponsored by NIST, NIST annual reports,

and other special publications appropriate to this grouping such as wall charts, pocket cards, and

bibliographies.

Applied Mathematics Series—Mathematical tables, manuals, and studies of special interest to physicists, engineers, chemists, biologists, mathematicians, computer programmers, and others engaged in

scientific and technical work.

National Standard Reference Data Series—Provides quantitative data on the physical and chemical

properties of materials, compiled from the world's literature and critically evaluated. Developed un- der a worldwide program coordinated by NIST under the authority of the National Standard Data

Act (Public Law 90-396). NOTE: The Journal of Physical and Chemical Reference Data (JPCRD)

is published quarterly for NIST by the American Chemical Society (ACS) and the American Insti- tute of Physics (AIP). Subscriptions, reprints, and supplements are available from ACS, 1155 Sixteenth St., NW, Washington, DC 20056.

Building Science Series—Disseminates technical information developed at the Institute on building

materials, components, systems, and whole structures. The series presents research results, test methods, and performance criteria related to the structural and environmental functions and the

durability and safety characteristics of building elements and systems.

Technical Notes—Studies or reports which are complete in themselves but restrictive in their treat- ment of a subject. Analogous to monographs but not so comprehensive in scope or definitive in

treatment of the subject area. Often serve as a vehicle for final reports of work performed at NIST

under the sponsorship of other government agencies.

Voluntary Product Standards—Developed under procedures published by the Department of Commerce in Part 10, Title 15, of the Code of Federal Regulations. The standards establish nationally

recognized requirements for products, and provide all concerned interests with a basis for common

understanding of the characteristics of the products. NIST administers this program as a supplement

to the activities of the private sector standardizing organizations.

Consumer Information Series—Practical information, based on NIST research and experience, covering areas of interest to the consumer. Easily understandable language and illustrations provide useful background knowledge for shopping in today's technological marketplace.

Order the above NIST publications from: Superintendent of Documents, Government Printing Office,

Washington, DC 20402.

Order the following NIST publications—FIPS and NISTIRs—from the National Technical Information

Service, Springfield, VA 22161.

Federal Information Processing Standards Publications (FIPS PUB)—Publications in this series collectively constitute the Federal Information Processing Standards Register. The Register serves as the official source of information in the Federal Government regarding standards issued by NIST

pursuant to the Federal Property and Administrative Services Act of 1949 as amended, Public Law

89-306 (79 Stat. 1127), and as implemented by Executive Order 11717 (38 FR 12315, dated May 11,

1973) and Part 6 of Title 15 CFR (Code of Federal Regulations).

NIST Interagency Reports (NISTIR)—A special series of interim or final reports on work performed

by NIST for outside sponsors (both government and non-government). In general, initial distribution is handled by the sponsor; public distribution is by the National Technical Information Service,

Springfield, VA 22161, in paper copy or microfiche form.
U.S. Department of Commerce

National Institute of Standards and Technology

Gaithersburg, MD 20899

Official Business

Penalty for Private Use $300