United States Department of Commerce

Technology Administration

National Institute of Standards and Technology

NIST Special Publication 910

Green 's Functions and Boundary Element Analysis for

Modeling of Mechanical Behavior of Advanced Materials

Edited by

J. R. Berger

V. K. Tewary

Proceedings of a workshop held

August 14-17, 1996 at the

National Institute of Standards and Technology

Boulder, Colorado

NIST Special Publication 910

Green's Functions and Boundary Element Analysis for

Modeling of Mechanical Behavior of Advanced Materials

Edited by

J. R. Berger

Division of Engineering

Colorado School of Mines

Golden, CO 80401

V. K. Tewary

Materials Reliability Division

National Institute of Standards and Technology

Boulder, CO 80303

Proceedings of a workshop held

August 14-17, 1996 at the

National Institute of Standards and Technology

Boulder, Colorado

December 1996

U.S. DEPARTMENT OF COMMERCE, Michael Kantor, Secretary

TECHNOLOGY ADMINISTRATION, Mary L. Good, Under Secretary for Technology

NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, Arati Prabhakar, Director
National Institute of Standards U.S. Government Printing Office

and Technology Washington: 1996

Special Publication 910

Natl. Inst. Stand. Technol. Spec. Publ. 910

160 pages (December 1996)

CODEN: NSPUE2

For sale by the Superintendent

of Documents

U.S. Government Printing Office

Washington, DC 20402
Contents

Foreword v Workshop Summary 1

Papers 5 Accurate Discretization of Integral Operators 7

Bradley K. Alpert

Efficient Modelling of Two and Three Dimensional Crack Growth

Using the Surface Integral and Finite Element Hybrid Method 17

Balkrishna S. Annigeri and William D. Keat

Boundary Element Analysis of Bimaterials Using

Anisotropic Elastic Green's Functions 29

John R. Berger

Green's Functions in Elastic Fracture Mechanics 43

T. A. Cruse

Mechanical Properties of Metal-Matrix Composites 51

L. C. Davis

Green's Function BEM Research at Rutgers CMAS and

its Applications to the Modeling of Mechanical and

Piezoelectrical Behaviors of Advanced Materials 59

Mitsunori Denda

Approximate Operators for Boundary Integral Equations

in Transient Elastodynamics 77

Thomas L. Geers and Brett A. Lewis

Simulation of the Electrochemical Machining Process

Using a 2D Fundamental Singular Solution 85

H. A. Nied

Boundary Elements and a Green's Function Library 101

F. J. Rizzo, P. A. Martin, and R. A. Roberts

Fracture and Crack Growth of Polymeric Composites

for Use in Dynamic Applications 113

Ronald A. L. Rorrer and Will Fehringer

Elastic Green's Functions for Anisotropic Solids 123

V. K. Tewary

Appendix A. Organizing Committee 129

Appendix B. Program 131

Appendix C. Participants 135

Appendix D. Abstracts 139

iii
Except where attributed to authors affiliated with the National Institute of

Standards and Technology (NIST), the content of individual sections of this volume has not been reviewed or edited by NIST; NIST, therefore, accepts no

responsibility for comments or recommendations therein. The mention of trade names in this volume implies no endorsement or recommendation by NIST.

N

iv
Foreword

The 1994 workshop Green's Functions and Boundary Element Analysis for

Modeling of Mechanical Behavior of Advanced Materials was organized by Dr. Vinod

Tewary and Prof. John Berger to demonstrate the potential of Green's functions and

boundary element methods in solving a broad range of practical materials science

problems. The workshop assembled an outstanding and diverse group of researchers

from the National Institute of Standards and Technology and other government

laboratories, universities, and industry to assess the current state of the art and to identify research opportunities. The technical presentations by the researchers were an

excellent review of recent progress in the field. The contributions of the participants

from industry helped all of the workshop participants to understand the technical areas

in which GF/BEM research can contribute significantly to improved materials and

materials processing methods for commercial use. They also defined specific problems, which will be the basis for follow-on research projects.

From the beginning, Dr. Tewary and Prof. Berger designed the workshop to nucleate these ongoing research projects involving teams of workshop participants

and their colleagues. Several focused research projects, most of them built around

a specific industrial technology need, were started shortly after the workshop. The

Green's Functions/Boundary Element Methods Consortium that encompasses these

projects recently created its own web site. For further information on the projects and

the members of the research teams, see their web site at http://www.ctcms.nist.gov/

—jberger/greens/.

The research projects continue under the Center for Theoretical and Computational Materials Science (CTCMS), which was established in 1994 to develop and

apply state-of-the-art theoretical and computational materials science techniques and

to help industry integrate them into technology development. The center's technical

projects emphasize the solution of problems in materials design, processing, and

application and the transfer of the results to industry through research collaborations.

Each of the center's active research projects typically involves a team of re- searchers from several institutions focused on an industrially important technical

problem defined by the industrial participants. Projects may also be affiliated with

established research groups. The CTCMS has an active web site that is used for col- laboration among center researchers and to provide general and scientific information

to the technical community. For more information on the CTCMS, see the web site at http : //www . ctcms . nist . gov

The project on Green's functions and boundary element analysis began, as most

projects begin, with a workshop to define a focused research area in which the center can make a meaningful contribution. Collaboration continues, with project team members working at their own institutions; in future workshops and meetings, the research

team members will review progress on achieving the goals of this workshop.

Dale Hall

Deputy Director (Acting)

Materials Science and Engineering Laboratory

National Institute of Standards and Technology

V

Green's Functions and Boundary Element Analysis

for Modeling of Mechanical Behavior of Advanced Materials

The use of Green's functions and boundary element analysis for modeling

the mechanical behavior of advanced materials was explored by representatives

of universities, industries, and national laboratories in a workshop. Discussion

groups following the presentation of papers led to the identification of research

topics of industrial interest: fracture and damage in heterogeneous and layered

materials, characterization of multiphase materials, nondestructive measurement

techniques, electrochemical machining technology, and analysis of the growth of

thin films.

Key words: boundary element analysis; computer-aided engineering; elec- trochemical machining technology; fracture; Green's functions;

nondestructive evaluation

Workshop Summary

On August 14 through 17, 1994, the National Institute of Standards and Technology

(NIST) in Boulder, Colorado hosted a workshop on Green's Functions and Boundary

Element Analysis for Modeling the Mechanical Behavior of Advanced Materials. The

organizing committee and program are given in Appendixes A and B. The workshop was

sponsored by the NIST Center for Theoretical and Computational Materials Science. Its goals were

1 . To further explore the use of elastic Green's functions in conjunction with boundary

element analysis for modeling the mechanical behavior of advanced materials;

2. To identify areas for research collaboration and technology transfer between NIST

and other organizations;

3. To plan a cooperative research and development program to develop efficient software

tools for modeling the mechanical performance of advanced materials.

The quality of presentations and subsequent discussions demonstrated that the workshop

attracted the top researchers in the field of Green's functions and boundary element analysis

as well as a collection of individuals with a technical need for advanced analysis. Of the

29 participants, 11 were from universities, 8 from industries, and 10 from national laboratories (see Appendix C). Thus, the presentations and discussions can be considered representative of the thinking in these three sectors. Abstracts of the presentations are included as Appendix D.

The discussion groups were based on subject matter and interest: fracture, quantitative

nondestructive evaluation (NDE), and computer-aided engineering (CAE). Each group was charged with identifying the key technology issues, establishing the relevance of these

1
issues to Green's functions and boundary element analysis, and identifying specific projects.

Ideally, each of the identified project teams included a representative from an industry, a

university, and a national laboratory. The key issue was the relevance of the technology

discussed at the workshop to specific problems in industry. The methods of Green's functions and boundary element analysis are useful tools, provided the problems are amenable

to the method. At the conclusion of the workshop, several projects were identified that precisely fit the proposed model.

The Green's function methods discussed at the workshop are particularly powerful

when applied to the problem of fracture. Since the singular behavior of stress is incorporated directly into the model, no special techniques are needed for simulating the singularity.

Therefore, the power of the Green's function method is in the ability to extract analytic

insights concerning the mechanics of material failure. The discussion group on fracture

identified the key technical issues as grain boundary and size effects, interfacial fracture,

diffuse cracking and damage in heterogeneous materials, near-interface cracking in thermal

barrier coatings, and process-induced damage during manufacturing. The group considered

which of these topical areas were amenable to analysis by Green's function methods and

identified two potential areas for developing specific projects:

1 . Fracture and damage in heterogeneous media, which involves consideration of multiple

crack sites, damage accumulation, the effect of rigid or ductile particles, the effect of

creep and stress-strain behavior, and the development of reliability models.

2. Fracture and fatigue of layered media, which involves the effects of ductility or brittle- ness of the bonding layer(s), high-heat fluxes (such as those experienced by thermal

barrier coatings), and interface geometry.

In both the project areas, the need for observation and experiments to guide and verify the

theoretical developments was considered crucial.

The Green's function methods have a long history in quantitative NDE because they can

represent point sources in infinite, semi-infinite, and finite domains. The quantitative NDE

working group identified two key technology issues:

1. Characterization of multiphase materials, which includes the forward problem as well

as the inverse problem. The inverse problem is especially difficult since the goal is the

detection of microfeatures in the heterogeneous material. This problem is relevant to modern materials since almost all advanced materials contain multiple phases and,

therefore, interfaces. The multiple phases and microstructure directly affect the macroscopic properties of the material. Therefore, the ability to sense and control microstructure is essential for advanced material applications. The specific problem identified under

this problem statement was the development of time-domain boundary elements for

materials characterization.

2. Modeling and optimization of measurement techniques, which relates to probe and sensor

design, sensitivity analysis, and optimization. Additionally, issues such as the probability

of detection and characterization of subsurface features could be addressed. This topic is relevant to the design of experiments, validation of experimental techniques, and the con- gruence of models and experiments. Two specific projects were identified that would be

helpful in solving this particular problem: The first was a library of Green's functions

that would give the experimentalist access to the full power of various Green's functions

for predicting the experimental response of a material under varying initial conditions;

thus, the experimentalist could play a number of "what if" games before beginning an
experimental procedure. The second project selected was the development of the tools necessary to link macroscopic properties to microstructural analysis through models

based on Green's functions and nondestructive measurements.

The CAE group identified process optimization and engineering of materials as the main

technical areas. Several candidate projects were proposed in which Green's function methods

might be useful:

1 . The further development of electrochemical machining technology requires increasingly

sophisticated models to predict the final shape of the machined part. Similar projects in traditional machining, grinding, and profiling of rolled gears were also discussed.

2. The analysis of the growth of thin films is well-suited to boundary element investigation

because the surface conditions of thin films significantly affect film growth. Both of these

projects could be approached through an integrated software product, which would

contain solid modeling capabilities, analysis based on Green's functions, adaptive meshing, optimal design routines, and visualization capability.

The papers submitted for publication follow, alphabetically by first author. Abstracts for

all presentations are included in Appendix D.

The editors thank the Center for Theoretical and Computational Materials Science for sponsoring the workshop. We also acknowledge the help and advice of Dale Hall, Harry

McHenry, and the members of the organizing committee in ensuring the success of the

workshop.

3

Papers

5

ACCURATE DISCRETIZATION OF INTEGRAL OPERATORS

BRADLEY K. ALPERT

Applied and Computational Mathematics Division National Institute of Standards and Technology, Boulder, Colorado 80303

Abstract. The numerical solution of integral equations involves reduction of equations on a function space to a finite system of equations, typically over a chosen set of basis functions (the boundary

element method). The accuracy of the computed solution is determined by the accuracy with which

the true solution can be represented by the finite basis, as well as the accuracy of projections onto the

basis. An alternative discretization technique, the Nystrom method, which reduces the integral equations to a finite system of equations by replacing each integral with a quadrature, can often be designed

to obtain better convergence. We present some quadrature techniques for handling kernel and solution

singularities arising in integral equations and give numerical examples of their accuracy. Key words, numerical integration, singular kernels, quadrature rules

A wide variety of problems in acoustics and structural analysis, electromagnetic scattering and

propagation, and other applications are conveniently formulated as integral equations. Such formulations are particularly appropriate when the problem domain is homogeneous or piecewise homogeneous and has complicated boundaries or is unbounded, where differential equation methods often

present difficulties. An integral equation to be solved numerically must be reduced to a finite-dimensional problem.

This reduction is equivalent to expressing the solution as a function of a finite set of unknowns, and

expressing the integral equation as a set of equations governing the unknowns.

In the mathematically simplest setting, a domain D is specified by its boundary F = 9£), a boundary function g and Green's function K are given, and an unknown boundary function / is to be determined according to a second-kind integral equation

Here x, y e and F is a closed curve in the plane (for a two-dimensional domain) or x, y e

and F is a closed surface in space. Under broad regularity conditions on K and F, and assuming that

(1) with g replaced by 0 has only the trivial solution, the theory of integral equations developed by

Fredholm assures us that (1) has a unique solution /.

1.1. Boundary Element Method. The boundary element method reduces ( 1

) to a numerical problem

by projection to the space spanned by a finite set of functions {(f)\, . . . ,</>„}. If this set is orthonormal,

E-mail address: alpert@boulder.nist.gov.

1. Introduction

(1)

7
the approximate solution

f(x) = J2fi(l>iix) (2)

1=1

is obtained by solving the system of equations

where gi is given by

fi + J2^'jfj=Si^ i = h...,n, (3)

7 =1

8i = j six) (piix) ds(x), i = 1,... ,n, (4)

and Kij is given by

' ^0 y)<f>iM^j(y)ds{x)ds{y), ij = \,..., n. (5)

The quality of the approximation f ^ f depends on the size of the differences

n Eg(x) = g{x) -^gi(pi(x), (6)

EKix,y) = K{x,y)-Y,Y^Kij(l>iix)(l>j{y). (7)

: i=l 7 = 1

In the process of computing the numerical solution, additional errors are typically introduced in the

computation of the integrals in (4) and (5). In most implementations of the boundary element method, the functions </>, are taken to be piecewise constant or piecewise linear functions on P. This choice, while making the computation of (4)

and (5) feasible, results in estimates of the norm of Eg of order 0(n'^), at least. The norm of Efc will

be even larger, in the typical case where K is singular at x = y. These inaccuracies will introduce

corresponding errors in /. If, furthermore, the boundary F has comers, the true solution / will be

singular there, and the norm of / — / will be larger than order

1 .2. Nystrom Method. The Nystrom method represents functions / and g by their values at selected

points {jf 1, . . . , x„ } on r , and replaces the integral of ( 1 ) by a quadrature using those points. We obtain

the system of equations

n fi+^ Wij fj = g(xi), / = 1, . . .

, n, (8)

7 = 1 where the Wij depend on the kernel K. The true solution / is approximated by /, where /(Xi) = fi,

i = 1 , . . . , n. If necessary, / can be determined at points other than x\, . . . ,Xnby interpolation.

The theory of convergence for Nystrom discretizations is well developed (see, for example, Mikhlin

[1]). It can be summarized as follows: the error of / is proportional to the quadrature error

f]wijf{xj)-fK{Xi,y)fiy)dsiy), (9)

8
where / is the true solution. The task in applying Nystrom's method, then, is the construction of

quadratures that are accurate for a class of functions broad enough to contain the solution or closely

approximate the solution. In the rest of the paper we introduce quadrature methods appropriate for

integral equations that arise in physical applications and provide examples for two dimensional domains. These methods have been developed by many authors (Alpert [2], [3], Kapur and Rokhlin [4], Kussmaul [5], Kress [6], Martensen [7], Rokhlin [8], Starr [9]).

2. Quadrature Techniques

In a situation in which both the Green's function K and solution / of (1) are smooth (i.e., possess

several continuous derivatives), classical quadrature techniques such as Simpson's rule and Gaussian

quadrature are quite adequate. In the usual setting, however, the Green's function is singular at x =

y and the solution is singular at comers of the domain. Good convergence and accuracy depend on

correct treatment of these singularities. In what follows, we shall assume that K takes one of the following forms:

K{x, y) = \og(\x - y\) R{x, y) + S{x, y), (10)

v(y) • (x — y) K{x,y)= 7 R{x,y) + S{x,y), i=2,3, (11)

\x - y\'

where v{y) is a unit vector, typically normal or tangent to the boundary at y, and R and S are functions

possessing several continuous derivatives with respect to x and y. Here |;c — y | is the distance between

points X and y, which in two dimensions is given by the formula

\x-y\ = y/{x,-yi)^ + {x2-y2)^. (12)

Equation (11) includes cases (depending on / and the number of dimensions) where the integral must

be interpreted as a Cauchy principal value integral. These forms occur for both equiUbrium and time

dependent problems (Laplace operators) and time harmonic problems (Helmholtz operators) in two

and three dimensions.

2.1. Smooth Boundaries. When the domain boundary F is smooth, the solution / of (1) can be

shown to be smooth. In this case the necessary integrals can be split into two parts, one smooth and

one containing a known singular function multiplied with a smooth function, and each part can be

treated with an appropriate quadrature. Given a smooth parametric representation of F,

xit) = (xiit),X2{t)), re[0,2jr], (13)

with the squared differential arc length a(r)- = x[{t)^ +x'2{tf positive for all t e [0, lit], the integral

of (1) becomes ^

Kix, y) f{y)ds{y) = ^

K(x(t),x{r)) f{x{x))a{T)dT. (14)

With this parametrization, by Taylor expansion we obtain

x{x) -x{t) = (x[ (r), x'^it)) 2 sin ^

+ O (sin-

^

j

, (15)

9
and

logd^ - y\) = log4sin2

^

+ \ogy/a(t)^ + 0(sm^[(r - 0/2]),

(x-y) (x[{t),x'it)) /x[it)x'{it)+x'it)x''(t) 1 T-t — Ill z ^

(16)

|jc-j|2 a(f)2 \ a(f)2 2 2

The expressions

n T — f T — r log4sin2^-, cot-^ (18)

are singular atr = t and smooth elsewhere; the remaining expressions on the right-hand side of (16)

and (17) are smooth everywhere. We remark that if v(y) is normal to F at y, then v(x)-{x[(t), ^2(0) =

0 and the cotangent singularity does not occur. For a smooth kernel K, the integral of (14) is well approximated by the trapezoidal rule, which has

convergence of order corresponding to the number of continuous derivatives of the integrand (spectral

convergence). This convergence follows from the fact that the exponentials

g-((2n-l)r g-!(2n-2)T

| ^ir ^i(2«-1)t

are integrated exactly by the 2n-point trapezoidal rule. For a kernel whose singularity can be isolated, as in (16) and (17), the smooth part can be integrated with the trapezoidal rule and the singular part can be integrated with a 2n-point quadrature rule that is exact for the functions

^T(T)e-'"^ a(T)e-'("-^>\ . . .

, a(r), a(T)e", . . .

, a{r)e'^"-^^\

where cr(T) is the known singular function. We choose equispaced nodes = Ttk/n, k = 0, . . .2n —

1 and obtain the quadrature rules

I.

2n ^ 2n-\

/(r)^ir^- V/(r,"), (19)

/ "

In Usix^ f{z)dr ^-T ql{t) fix;:), (20)

f\oi(^-=A fiT)dr^-Tr'^(t)f(r','), (21)

^0 V ^ / «

where

1 1 ^^(0 = - - cosn(r - tO - 2 J] y

cos/(r - rO, (22)

1=1 n-\

r"^it) = sinn(r - r^) + 2^sin/(? - r^), (23)

1=1 for A: = 0, . . . , 2n - 1 . We remark that if t is one of the quadrature nodes, as in the usual case, the

expression for simplifies; in particular, t" — r"

10
which corresponds to the trapezoidal rule with n nodes such that t = x" lies halfway between two

nodes.

These quadratures are exceedingly powerful where they are appropriate, yielding exponential con- vergence for boundaries that are analytic closed curves. As we shall see in the numerical examples,

other methods do not perform nearly as well. In fact, depending on the wavelength that characterizes a problem and on what quantities must be computed, it may introduce negligible modeUng error

to round the comers so that these quadratures can be exploited. For situations where the boundary

cannot be suitably approximated by such a smooth curve, however, we turn to quadratures tailored to boundaries with comers.

2.2. Boundaries with Corners. As is well known, integral equations on domains with comers generally have solutions that are singular at the comers. We again use the parametric representation of T

given by

xit) = X2{t)), t € [0, 27r], (24)

which we now assume to be smooth (piecewise) only on the intervals t € f,], / = 1, . . . , m,

with 0 = ro < ^1 < < tm = In . The solution / will have the form

f{x{t)) = it- r,_i)"' Qi(t) + iu - r)^' Riit) + Siit), t e [r,_i, f,], (25)

where Qi, /?, , and are smooth functions, for / = 1, . . . ,m. The exponents a, and depend on the

kernel K and the angles of the comers of F at and r, , respectively. For a variety of situations in acoustics and electromagnetics, these exponents are known analytically.

2.2.1. Known solution singularities. The weights of the trapezoidal mle can be 'corrected' to accurately integrate functions with endpoint singularities [8]. The idea is to adjust a few quadrature weights

at the singular end so that the integration is exact for functions

f{t) = ait)t' +tJ, (26)

for small /, j, where a is the singular part. In particular, the linear system of equations

1 -.2k

fit) dt = Zif) + -T fii/n) w), (27)

7 = 1

for /(O = ait) t' and fit) = x',i = 1, . . . , ^, is solved for u;^, . . . , Wj^. Here r„'(/) is the trapezoidal mle with the left-end value /(O) omitted and with the right end corrected to high order for

smooth functions. It has been shown [8] that this 'corrected' trapezoidal rule has kth order conver- gence for functions fit) — ait) Rit) + Sit), for smooth R and S. This means that quadrature error

is of order Oin"^). Furthermore, the weights if;" , . . . , itij^ approach hmiting values wi, . . . , W2k as n -> oo and can be replaced by these values for large n. For functions such as / in (25) that are singular at both ends of the integration interval, weights for the appropriate singularity can be applied at each end.

In solving integral equations on domains with comers, we must contend with the singularities of

both the kemel and the solution. For the integral

f

Kixit),xix))fixix))aix)dx, (28)

the kemel K is split into a singular part K\ and smooth part K2 (as described in section 2. 1 above). The

integral involving K2 is computed as just described. If t G , r, ], however, the integral involving

11
K\ will have a singularity in the middle of the interval. The appropriate procedure here is to split

the interval into t] and [t, ti] and apply the two-ended singularity quadrature on each interval. Likewise, if t is outside , but close to the interval, the correction weights w", . . . , lUj;^ should

be obtained by replacing f by K2 f in (27).

This procedure requires the solution of a rather large number of small scale (i.e., 2k x 2k) linear

systems of equations for the quadrature weights. It does, however, avoid inaccuracies associated with

an uncorrected uniform discretization of the boundary. It also avoids clustering nodes near comers,

which would increase the number of unknowns to be solved in obtaining the solution /.

2.2.2. Unknown solution singularities. For some integral equations it is not possible to obtain the

form of the singularities in (25). This situation is common in three-dimensional problems, where cor- ners are formed by the intersection of more than two smooth surfaces. Accurate solutions can then be

obtained only by clustering nodes near comers. We demonstrate in two dimensions that good conver- gence can be obtained for these problems.

The idea is to make a change of variables for r to eliminate the endpoint singularities. We define

(2^-i-l)! f , Ms)='—7-r- [ril-r)tdr, (29)

(^!)^ Jo where the constant is chosen so that Ak(\) = 1, and for r e [ti-i, ti] we make the change of variables

r =ti-i+iti -ti_i)Akis). (30)

Then (28) becomes

(2k + IV

K(x(t),x(T(s))) f(x(r(s)))(ti -ti-0 ^, ,^/ [s(\-s)rds. (31)

Now, assuming / is given by (25), the integrand is of order O(5"' ) for s near 0, and of order

0((1 — 5)^'*^+'^+'''

) for s near 1. A function that is / times continuously differentiable and vanishes

along with / — 1 of its derivatives at the endpoints of an interval can be integrated on the interval with

the simple trapezoidal rule, with an error of order 0{n~'~^). In our case, / = a, A: + k -|- a,, which

by proper choice of k can be made arbitrarily large since a, > — 1 . In practice, this procedure works

well for k as large as 6 or 8; for larger values of k, round-off error causes quadrature nodes near the

endpoints to become indistinguishable. The procedure, which for = 8 requires more than three

times as many points as the trapezoidal mle with the same point spacing in the center of the interval, can be improved by applying the change of variables only near the interval ends. For example, we

define

1 (2^ + 1)! , , . 52^+' (1-5) ()t!)2 ffJo Jo [r(8-r)Ydrdt, 0<s<8,

s-8/2 . 1 . (32)

, 5 < s < I — 8, \-8

l-AksO-s), l-8<s<l,

where the chosen constants yield continuity of A^s and its first k + 1 derivatives. The corresponding

change of variables yields a very effective quadrature technique, particularly for large problems.

Once again, the singularity of the kernel K must be handled separately. After the change of vari- ables eliminates the endpoint singularities, the integration interval must be divided in two and correction made for the singularity of the kernel on one end of each subinterval. This can be accomplished

with the method of section 2.2.1.

12
The evaluation near an edge of the integral of sources on that edge, as occurs near a comer, must be

treated with some care. Although the field point will not coincide with any of the source points (which

would produce a singularity), the integrand will be nearly singular. Without special provisions, the

quadrature error would be large. We have successfully addressed this issue by effectively applying the

quadratures on a finer grid, where the values of the sources at the additional locations are determined

by local interpolation. By this procedure, the sharp behavior of the kernel is captured directly, while

exploiting the smoothness of the solution, without increasing the number of unknowns.

The treatment of surfaces is not fundamentally different from the methods established for curves.

Integrals on rectangles can be treated one dimension at a time, which requires the edge and comer

techniques discussed above. For a smooth surface topologically equivalent to a rectangle, a smooth

map between the surface and a rectangle leads to a quadrature for the surface. Similarly, for a topological cylinder, a combination of the quadratures for smooth closed curves and those for segments

is appropriate. While the actual implementation of this method is somewhat involved, no additional

mathematics is required.

In this section we present a several numerical examples of the convergence of the quadratures dis- cussed above. Our first example is for an integral equation over a (fairly) smooth boundary approximating a rectangle, with a Cauchy kemel. The equation is (1), with

where x = ix[,X2) and v(y) is the unit tangent to the boundary at y. The boundary F consists of a

rectangle with sides xi =0,xi = 2, ^2 = 0, = 0.5, except that the comers are replaced by quarter

circles of radius s = 0.1, 0.2. The boundary is discretized with n nodes equispaced in arc length.

The quadratures in this case are based simply on the staggered trapezoidal mle for both the Cauchy

singularity and the smooth part, where the two need not be separated. In table 1 we observe roughly

second-order convergence, anticipated because F and its first derivative are continuous. Nevertheless,

we also observe that as the comer gets sharper (e = 0. 1), the accuracy deteriorates. This suggests that the comer should be treated explicitly. The computed solutions for various n are shown in figure 1

Table 1 . The solution of integral equation (1) with kemel defined in (33) and righthand side of (34) is computed by a staggered trapezoidal mle. The relative error of

the computed solution for various e and n is obtained by comparison with the solution

3. Quadratures on Surfaces

4. Numerical Examples

gix) = cos{\Ox lit) + \0x2(t))

(33)

(34)

forn = 1200.

n £ = 0.2 £ =0.1

200

300

400

500

600

1.14e-01

5.26e-03

1.41e-02

3.14e-03

7.68e-04

3.58e+00

1.73e-01

2.84e-01

l.Ole-01

2.87e-F00

13
f(t) X

10"^ Cauchy Singularity Solution (eps = 0.2)

700.00

600.00

500.00

400.00

300.00

200.00

100.00

-0.00

-100.00

-200.00

-300.00

-400.00

-500.00

0.00

—11

H

II

IIII —h

II

IIII

II

1

Iri A

)1

1 1

i U

H—

n

vJl

1 1

If

1.00 2.00 3.00 4.00

n= 200

n=400"

n=666'

Perimeter

Figure 1 . The computed solutions for the "rounded rectangle" with kernel defined

in (33) and right-hand side of (34) for 3 discretizations are shown.

In our second example we examine the convergence associated with the quadratures for unknown

singularities. Integrands with two different singularities at 0 were integrated with the change of vari- ables in (32), and compared with exact values, as shown in table 2. In our final example we solve a first-kind integral equation on a rectangular domain. The equation

is ^

log(|;c - y\) f{y) dl(y) = sin(xi + X2), (35)

14
Logarithmic Singularity Solution

Figure 2. Two computed solutions, with rectangle perimeter discretizations of 52

and 104 points, for (35) are shown.

where T is the rectangle with sides = 0, xi = 10, X2 = 0, = 1- For each side of the rectangle,

the change of variables (30) was used with k — 4. The solution, which develops singularities at the

comers, is shown in figure 2.

5. Summary

This paper had the goal of illustrating some of the quadrature techniques that can be apphed within

the framework of the Nystrom method. We argue that the numerical difficulties that arise in any in- tegral equation method, particularly the singularities of the kernel and solution, can be handled di- rectly and effectively with the Nystrom method, yielding good convergence and excellent accuracy.

15
Table 2. The integral /q f(x) dx for two different integrands / is computed with

the trapezoidal rule with the change of variables in (32) as well as evaluated analytically. Note that each integrand has 500/7r ^ 159 oscillations in the integration

interval. The relative errors of the quadratures are shown.

/(x) =cos(1000x)log(jc)

k == 2 k == 6 n 8 = 0.1 <5 = 0.2 8 =0.1 5 = 0.2

200 2.81e+00 1.36e+01 1.90e+00 5.28e+01

250 5.04e-03 4.39e-03 5.31e-04 4.89e-04

300 3.04e-04 1.55e-04 7.80e-07 7.10e-08

350 2.25e-06 1.66e-05 3.71e-08 3.03e-09

400 1.21e-05 l.OOe-06 4.72e-09 1.90e-ll

fix) = cos(1000;c)x-'/^

k == 2 k == 6

n 8 =0.1 <5 =0.2 8 =0.1 5 = 0.2

200 1.71e-01 4.09e+00 1.25e-01 1.80e+00

250 7.21e-03 2.83e-03 6.31e-05 2.28e-05

300 4.84e-03 1.82e-03 9.89e-06 9.28e-07

350 3.55e-03 1.33e-03 5.35e-06 5.02e-07

400 2.71e-03 1.02e-03 3.14e-06 2.95e-07

We presented some numerical examples to demonstrate these techniques. More detail on the theory

and convergence of the techniques can be obtained by consulting the references.

' ^ • References

[1] Mikhlin, S. G. Integral Equations and Their Applications to Certain Problems in Mechanics, Mathematical Physics,

and Technology. Pergamon Press, New York, 1957.

[2] Alpert, B. Rapidly-convergent quadratures for integral operators with singular kernels. Technical report LBL-30091,

Lawrence Berkeley Laboratory, University of California, Berkeley, CA, 1990.

[3] Alpert, B. High-order quadratures for integral operators with singular kernels. J. Comput. Appl. Math., 60:367-378,

1995.

[4] Kapur, S. and RokhUn, V. High-order corrected trapezoidal rules for singular functions. SIAM J. Numer Anal, to appear, 1996.

[5] Kussmaul, R. Ein numerisches Verfahren zur Losung des Neumannschen Aussenraumproblems fiir die Helmholtzsche

Schwingungsgleichung. Computing, 4:246-273, 1969.

[6] Kress, R. A Nystrom method for boundary integral equations in domains with corners. Numer Math., 58:145-161,

1990.

[7] Martensen, E. Uber eine Methode zum raumlichen Neumannschen Problem mit einer Anwendung fiir torusartige Berandungen. Acta Math., 109:75-135, 1963.

[8] Rokhlin, V. End-point corrected trapezoidal quadrature rules for singular functions. Comput. Math. w. Appl, 20:51-

62, 1990.

[9] Starr, H. P. On the Numerical Solution of One-Dimensional Integral and Differential Equations. Ph.D. thesis, Yale

University, New Haven, CT, 1991.

16
EFFICIENT MODELLING OF TWO AND THREE DIMENSIONAL

CRACK GROWTH USING THE SURFACE INTEGRAL

AND FINITE ELEMENT HYBRID METHOD

Balkrishna S. Annigeri

United Technologies Research Center

East Hartford, Connecticut 06108

U.S.A.

William D. Keat

Clarkson University

Potsdam, New York, 13699

U.S.A.

SUMMARY

The need for modelling crack growth in materials is important for assessing the fatigue life of

structural components. In a variety of industries: aerospace, nuclear, automobile and others; fatigue

crack growth can be a significant source of failure of structural components. This paper describes

recent advances in modelling crack growth using the Surface- Integral and Finite Element Hybrid

method. Cracks in an infinite or semi-infinite domain are modelled using a continuous distribution

of dislocations in 2D or force multipoles in 3D resulting in a singular integral formulation. The

uncracked body is modelled using finite elements. The coupling of the surface-integral and finite element models is obtained by traction and displacement matching on the external and internal

boundaries of the finite plate with crack(s). A significant advantage is that the finite element mesh

remains fixed as the crack propagates. The remeshing then is required only on the surface integral

discretization which is accomplished automatically without user intervention. This paper describes

the recent effort of modelling through cracks in layered materials and modelling of embedded and

surface cracks in 3D. Numerical predictions have shown good correlation with experimental data

and other analytical solutions. The research efforts have resulted in the development of the SAFE2D and SAFE-3D codes for effective modelling of crack growth.

INTRODUCTION

Structural components can fail due to a variety of reasons such as: excessive deflections, excessive

yielding or inelastic deformation of the material, fatigue fracture initiation and growth, chemical

degradation such as oxidation, sulphidation etc. A significant portion of structural failures ,how17
ever, are due to fatigue loading which can result in initiation of cracks, subcritical crack growth and

eventually failure due to dynamic crack growth. Fatigue is a process wherein the repeated loading

causes growth of voids in the material and their coalescence which leads to the formation of microcracks. These microcracks can link with each other to form macrocracks which can be observed in an engineering sense and can subsequently grow. The useful life N of a structural component can be defined as the sum of the number of cycles (Ni) required for crack initiation and the number for

cycles (Np) required for crack propagation. The relative magnitudes of Ni and Np depend on the

material under consideration, the ductility of the material, flaw distribution, surface finish, residual

stresses, etc. Materials inherently have flaws such as inclusions, voids, porosities, microcracks etc. These flaws are usually the result of the processes used for manufacturing structural components.

Determining the service life of structural components thus requires the ability to efficiently model

crack propagation in complex geometries. This capability is utilized for predicting the total useful

life of structural components and for determining inspection intervals.

There are several methods that can be used to predict crack propagation. For simple geometries of

the structural part and for simple shapes such as an edge through crack or embedded through crack

in 2D, or semi- circular, semi-elliptical surface cracks and circular or elliptical embedded cracks in

3D, use can be made of analytical solutions for the Stress Intensity Factors (Ref. 1). For complex

geometries and loading and for crack propagation, recourse has to be made to numerical solutions

for accurately determining the stress intensity factors and predicting fatigue life. Several methods

exist for modelling crack propagation such as the finite element, boundary element and surface

integral methods. The surface integral method which is used here is based on modelling cracks as a displacement discontinuity distribution using dislocations in 2D or force dipoles in 3D. The surface

integral method has been coupled with the finite element method resulting in the S.urface-Integral

And Finite Element (SAFE) hybrid method. The formulation of the governing equations is given

in the next section. Results for layered materials is presented for 2D through cracks and for 3D

part through cracks.

FORMULATION OF THE GOVERNING EQUATIONS

The development of the Surface Integral and Finite Element Hybrid Method has been presented in

earlier papers[2-10]. A summary of the of the formulation is given in the following. Consider the

problem of a plate with a crack as shown in Fig. 1. The plate is loaded by a load R along its external

boundary and by a traction T along the crack. The objective is to determine the Stress Intensity

Factor (SIF) at the crack tips. The plate with the crack can be represented as the sum of a plate

without the crack and an infinite plate with a crack; ensuring that tractions and displacements are matched at the external (and internal) boundaries. The plate without the crack does not contain

any singularities and hence can be effectively modelled using a finite element model (a boundary

element model can also be used) and the infinite plate with the crack is modelled using a surface

integral model. The infinite plate with the crack, however, needs a correction along the boundaries

of the actual plate. This is achieved by an appropriate boundary load correction vector R'^ which is computed along the boundaries of the finite plate using the surface integral model and subtracted

from the load in the finite element model to ensure that the correct total load is applied to the

boundary. Similarly, the traction vector acting along the crack line which is due to the forces

18
R — acting on the uncracked body is computed using the finite element model and applied to the surface integral model. The governing equation for the hybrid method is obtained as follows.

The finite element equation for the uncracked plate is given by:

[K] {C/^^} = {R} - {R^) (1)

[K] — finite element stiffness matrix of the uncracked body; {C/^^} = vector of unknown finite element nodal displacements for the uncracked body; \R\ = vector of applied loads; {R'^} = load

correction vector for the boundary due to the presence of the crack.

The crack in an infinite domain is represented by using a continuous distribution of dislocations

for 2D applications leading to a singular integral equation formulation [3]. For 3D applications,

the crack is represented by a continuous distribution of force dipoles [7]. The dislocations or force

dipoles are events that cause a displacement discontinuity. There is thus no need for modelling

the two surfaces of the cracks separately as is done in a conventional finite element or a boundary

element model. An analogy of this in aerodynamics is the use of fluid source and sink pairs for

modelling flow over thin airfoils in potential flow theory, commonly called a lifting surface model.

Only the mid surface of the airfoil is modelled for thin airfoils.

The 2D integral equation for an infinite domain with a crack is given by:

ti{xo) = i'th component of traction at xq due to a distribution of dislocations along the crack surface;

Uj = the j'th component of the unit normal to the crack surface at xq; T\j{x,xo) = fundamental

stress solution for a dislocation(kernel function); fii{x) = /'th component of the dislocation density

at x; Sc = surface area of the crack.

The 3D integral equation for an infinite domain is given by:

where 5i{x) denotes the /'th component of displacement discontinuity on the crack surface at location X. The singular integral equation is integrated numerically in the Cauchy principal value sense and the discretized form of equations (2) and (3) are both represented by:

[C] = coefficient matrix for the singular integral equation; {F) = vector of dislocation density

amplitudes at the interpolation points along the crack surface (for 2D) or the vector of nodal

displacement discontinuities (for 3D); {T} = vector of applied tractions along the crack; {T'^} =

traction vector for the crack surfaces due to the applied load R- R"" on the finite element model.

Equations (1) and (4) are fully coupled through the boundary force matrix G {{R^} = [G]{F})

and the stress feedback matrix S {{T^} = [S]{U^^}) which leads to the following representation

(2)

'c (3)

[C] {F} = {T} - {T=} (4)

19
for the hybrid formulation:

' K r; 1 f TI^E 1 f 1 (5)

' K

S ?1 {?}

In equation (5), the variable U^^ represents the continuous displacement field in the uncracked

body (the finite element model). To determine the total displacement field, the displacement in

the surface integral model is added so that the displacement boundary conditions can be applied,

as follows:

{U} = {U^^} + {U^'} (6)

where {U^^} is the vector of unknown surface integral displacements occurring at the finite element

nodal positions. {U^^} is obtained by evaluating the integral equation for displacements [3]:

{U''} = [L] {F} (7)

where [L] = is the displacement matrix. Thus equations(5), (6) and (7) can be combined to form

the governing equation for the hybrid formulation:

K G-KL

S C-SL (8)

In equation (8) arbitrary force, traction and displacement boundary conditions can be imposed.

The sparsity and symmetry of the stiffness matrix is taken fuU advantage of using a special solution

scheme which has been reported in [4]. The different matrices K, G, S and C are stored independently with only the upper half of the K matrix being stored. The K matrix is factored only once

using the LDL^ transformation; the results for various load/time steps are obtained by solving the

system of equations [4]. This scheme has been implemented in the SAFE-2D and SAFE-3D Hybrid

computer codes.

SAFE Codes- Capabilities

The SAFE-2D and SAFE-3D codes at present, have the following capabilities:

20
SAFE-2D

.single and multiple cracks

.surface and embedded cracks

.curved cracks

.branched cracks

.multiple crack origins

.cracks emanating from holes

.shear bands (non-opening cracks)

.isothermal isotropy

.thermal strains

.temperature dependent material properties

.automatic crack growth

.fatigue crack growth models

.layered materials

.variable thickness finite elements

.pre and post processors (integrated with PATRAN)

.crack tip plasticity models

•documentation available

SAFE-3D

.internal and surface cracks

.irregular shaped crack growth

.planar cracks

.isothermal isotropy

.layered materials

.pre and post processors (integrated with PATRAN)

.automatic crack growth

.fatigue crack growth models

.documentation available

Material Library

At present, the SAFE codes have several fatigue crack growth models built in the code to be used

with the user supplied material data.

Computers

SAFE-2D and SAFE-3D codes work on UNIX based workstations/computers, VAX systems, CONVEX and CRAY computers.

21
RESULTS

2D Analysis- Layered Materials

Many structural materials have layered construction or are used with coatings. Examples are thermal barrier coatings (TBC) on turbine airfoils, hardface coatings for wear resistance, ceramic

coatings on seals etc. There are also materials with graded material properties. An example of

graded materials is an abradable rub strip material used in gas turbine engines. It is important

to know how cracks initiate and evolve in such materials so that the structural durability can be

assessed. In this paper, results for a crack near a bi-material interface are presented. The geometry

of an angled crack near the bi-material interface is shown in Fig. 2. The aim of the modelling is to obtain stress intensity factors at the two crack tips. The influence function used in determining the

coefficient matrix C is that of an angled dislocation near a bi-material interface [14]. The results of

the analysis are presented in Fig. 2. The length of the crack is fixed but the angle with respect to the interface is varied. As can be seen, the results of the SAFE-2D code are in excellent agreement

with the results of Erdogan-Aksogan and Narendran [15].

The SAFE-2D code has been developed for modelling propagation of through cracks in structural

components. There are several theories available for modelling crack growth such as the maximum

circumferential stress theory, the minimum strain energy density theory, and the maximum energy

release rate theory . In the maximum circumferential stress theory, the crack propagates perpendicular to the maximum circumferential tensile stress direction in the neighborhood of the crack

tip. In the SAFE-2D code two criteria are available for propagating cracks. One is the maximum

circumferential stress theory based on the stress at the existing crack tip ; the second is based on Kii = 0 at the end of the incremented crack tip. The first option is suitable for Mode I and Mode II cracks and involves less computations as the direction of crack growth is determined from the state of stress at the current crack tip according to equations given in [3]. The second option is suitable

especially for modelling propagation of shear bands (non-opening cracks) and is computationally

more expensive as a boundary value problem has to be solved at every iteration for determining

the direction at which Kn = 0. Many crack growth problems have been modelled and the results

have been published earlier [2,4,5,7,9,16].

Results for a 3D geometry.

Fatigue growth of a surface crack in a thick plate under tension (see Fig. 3). has been obtained

using the SAFE-3D code. The SAFE-3D results are in good agreement with experimental data(Fig.

4). A Paris type fatigue law was used using a two parameter growth model; details are provided in

[16]. Several results for three dimensional problems are provided in [7,8,9,10,16]; good correlation

with analytical and experimental results has been shown.

CONCLUSIONS

This paper summarizes the formulation and the modelling capabilities of the SAFE-2D and the

22
SAFE-3D codes for fracture mechanics. The important advantage of the method is that cracks

and their associated evolution can be modelled very elegantly and efficiently by the use of integral

equations. A significant advantage is that only remeshing along the crack surface is required as the crack grows thus eliminating remeshing of the finite element mesh for the uncracked body.

The finite element mesh has to be able to model the applied loads and any boundary correction

but does not have to model any singular fields. This method has been successfully implemented

in the SAFE-2D and the SAFE-3D codes and integrated with the PATRAN code for pre- and

post-processing.

ACKNOWLEDGEMENTS

The research performed in this paper was supported by a corporate sponsored research program

by the United Technologies Research Center.

23
REFERENCES

1. Tada, H.; Paris P. and Irwin G.: The Stress Analysis of Cracks Handbook, Del Research

Corporation, Hellertown Pennsylvania, 1973.

2. Annigeri, B.S.: Surface Integral Finite Element Hybrid Method For Localized Problems in

Continuum Mechanics, Sc.D. Thesis, Department of Mechanical Engineering, M.I.T., 1984.

3. Annigeri, B.S.; Cleary M.P.: Surface integral finite element hybrid (SIFEH) method for

fracture mechanics. International Journal for Numerical Methods in Engineering, Vol.20,

869-885, 1984.

4. Annigeri, B.S.; Cleary M.P.: Quasi-static fracture propagation using the surface integral

finite element hybrid method. ASME PVP Vol.85, 1984.

5. Annigeri, B.S.: Effective modelling of stationary and propagating cracks using the surface

integral and finite element hybrid method. ASME AMD Vol.72, 1985.

6. Annigeri, B.S.: Thermoelastic fracture analysis using the surface integral and finite element

hybrid method. Presented at the ICES-88 Conference, Atlanta, Georgia, 1988.

7. Keat, W.D.; Annigeri B.S.; Cleary, M.P.: Surface integral and finite element hybrid method

for two and three dimensional fracture mechanics analysis.

International Journal of Fracture,Vol.36, 35-53, 1988.

8. Keat, W.D.: Surface Integral and Finite Element Hybrid Method For the Analysis of Three

Dimensional Fractures, Ph.D. Thesis, Department of Mechanical Engineering, M.I.T., 1989.

9. Annigeri, B.S.; Keat, W.D.; Cleary, M.P.: Fracture mechanics research using the surface

integral and finite element hybrid method. Proceedings of First Joint Japan/US Symposium

on Boundary Element Methods, University of Tokyo, Pergamon Press, 191-202, 1988.

10. Keat, W.D.; Cleary, M.P.: Analysis of 3-D near-interface fractures in bounded, heteroge- neous domains using the surface integral and finite element hybrid method. Proceedings

of International Symposium on Boundary Element Methods, United Technologies Research

Center, October 1989, Springer Verlag , in press.

11. Cruse, T.A.; Meyers, G.J.; Wilson, R.B.: Fatigue growth of surface cracks. ASTM STP 631,

American Society for Testing Materials, 174-189, 1977.

12. Raju, I.S.; Newman J.C.: Stress-intensity factors for a wide range of semi-elliptical surface

cracks in finite thickness plates. Engineering Fracture Mechanics, Vol.11, 817-829, 1979.

13. Cruse, T.A.; Besuner, P.M.: Residual life prediction for surface cracks in complex structural

details. Journal of Aircraft, Vol.4, 369-375, 1975.

14. Dundurs, J., Mura T.: Interaction between an edge dislocation and a circular inclusion, J. Mech. and Phys. Solids, 13, 141-147, 1965

15. Narendran V.M.: Elastostatic interaction of multiple arbitrarily shaped cracks in plane

inhomogeneous regions, S.M. thesis, MIT, 1981.

24
16. Annigeri B.S., Keat W.D.: Two and Three Dimensional Crack Growth Using the Surface

Integral and Finite Element Hybrid Method, IABEM-90 Symposium, Rome Proceedings,

Springer Verlag, October, 1990.

25
-7^

ACTUAL

PROBLEM

1 1

FINITE PLATE

WITHOUT CRACK

DISCRETIZED BY

FINITE ELEMENT MODEL

FINITE PLATE WITH

CRACK CUT OUT FROM

AN INFINITE REGION

DISCRETIZED BY

SURFACE INTEGRAL MODEL

Fig. 1 Superposition of finite element and surface integral models

26
30 60 90

6 - degrees

Fig. 2 Stress intensity factors for an angled crack in a bi-material panel

27
TEST SPECIMEN

7_

6 Y 7

HYBRID MODEL:

4a FINITE ELEMENTS

60 CRACK ELEMENTS

/

Fig. 3 Specimen geometry and corresponding model

used in crack growth study

(1 .250 in.) ^- H

^£ (0.294 in.)

0.60

« 0.40

GCo

o 0.30

0.20

0.1

— WEIGHT-FUNCTION METHOD ~— TEST RESULTS

HYBRID METHOD

00 5 10 15 20 25 CYCLES X 10 ^ FROM PRE-CRACK

Fig. 4 Comparison of SAFE-3D and experimental results for

growth of a surface crack in a thick plate

28
Boundary Element Analysis of Bimaterials Using Anisotropic

Elastic Green's Functions

John R. Berger

Division of Engineering

Colorado School of Mines

Golden, Colorado 80401

Abstract

The boundary integral equations incorporating the Green's function for

anisotropic solids containing planar interfaces are presented. The fundamental

displacement and traction solutions are determined from the displacement Green's

function of Tewary, Wagoner, and Hirth [Journal ofMaterials Research, Vol. 4, pp. 113-

123]. The fundamental solutions are shown to numerically degenerate to the Kelvin

solution in the homogeneous, isotropic limit. The boundary integral equations are formulated with the use of constant boundary elements. The constant elements allow for

analytic evaluation of the boundary integrals. The application of the method is demonstrated by analyzing a copper-nickel system subjected to mechanical load.

1 . Introduction

The analysis of deformation near interfaces in solids has received renewed interest

due to reliability issues in electronic packaging. Here we present a boundary element

formulation for anisotropic interface problems incorporating the fundamental solution of

Tewary, Wagoner, and Hirth [1]. This fundamental solution for the general anisotropic

interface problem allows us to model the behavior of stresses and displacements near the

interface exactly. No discretization of the interface is required since this behavior is explicitly incorporated in the fundamental solution. A similar approach was taken in [10,

1 1] for the isotropic interface problem.

The work presented here follows the work of Cruse [4] and Snyder and Cruse [8,

9] who implemented a complex-variable fundamental solution for the crack problem in a homogeneous, anisotropic plate. They found such a formulation to be computationally

efficient as well as providing exact modeling of the singular fields near the crack tips.

The use of enriched elements [14-16] has received a good deal of attention in the finite

element literature for use in modeling singular fields; however, the use of special Green's

functions for singular fields has been limited in boundary element formulations. The

formulation presented here is useful for boundary element analysis of interface problems

where the interface is flaw-free. Future work will focus on the special Green's functions

associated with interface cracks in anisotropic solids [13]. In this paper the fundamental displacement and traction solutions are presented as determined from the Green's function given in [1]. We note that four solutions are

actually needed depending on the relative location of the source and field points in the

two anisotropic solids. The fundamental displacement solutions U are of the general

form

u-Xr«'iog(z«- zj (1.1)

29
where the 7^'^ matrix is a function of the elastic constants in either material A or B and

the a roots of the Stroh determinant [2], a = 1, 2, 3. The complex coordinates z, z' of the

field and source points are defined in the classical Lekhnitskii form [3]. The traction

fundamental solution is then derived from the displacement fundamental solution. We

compare the degenerate isotropic form of these fundamental displacement and traction

solutions with the Kelvin isotropic solution.

The details of implementing the new fundamental solutions into the boundary

integral equation

CijUj +

j TijUjdT =

j UijtjdT ( 1 .2)

are given where the integrals are evaluated in the complex plane. We focus on the

constant element case to allow for analytic evaluation of the integrals.

Throughout the paper matrices and vectors are indicated either by a bold quantity

or by subscript notation. Unless otherwise stated summation is implied over repeated

(dummy) indices. We shall use the rectangular Cartesian coordinates Xj, Xj, X3 and the

elastic state is assumed to be independent of X3. A superscript asterisk (*) on a complex

quantity indicates a complex conjugate, and a primed coordinate will indicate the

coordinate system associated with a point source.

2. Anisotropic Fundamental Solution

The problem under consideration is shown in Fig. 1. In [1], Tewary, Wagoner,

and Hirth determine the Green's function for a composite solid with a planar interface by

taking Fourier transforms of the elastic equations. The Fourier transform method avoids

the problems associated with solving a complex eigenvalue problem as arises in the

analyses of Eshelby [121 and Stroh [21. The presence of an interface in the composite

solid necessitates solving the elastic equilibrium equations

^m-Tr± = -f. (2.1)

in both the upper half plane (material A) and the lower half plane (material B). In eq.

(2.1), Cj^y, are the elastic constants, m, are the displacements, and_^ are the forces at (Xi,X2).

By definition, the Green's function solves the equation

where 6(x) is the Dirac delta function, 5.^. is the Kronecker delta, the primed coordinates

represent the location of the point source, and the superscripts are for either material A or B.

Eq. (2.2) is solved assuming perfect conditions at the interface between materials

A and B,

".1,.=o="'U- (2.3)

(2.4)

30
The Green's function determined through the Fourier transforms of eq. (2.2) subject to the

boundary conditions of eqs. (2.3) and (2.4) is found in [1]. Since the point force and the

field point may be located in either material A or B there are in fact four parts to the

Green's function. The displacement field is given by the real part of the Green's functions

given in [1],

u"'=- 7R4lr(/'^)iog(^^-z„'')+Sr(p^)Q^iog(z^-z^'**)j, (2.5) ^ [a aP J U'''=- ^Re{lr'(p^)Q^og(.f -z,"*)!, (2.6)

ap

U^''= -Mlr(p'a)Qflog{z^-z^')\, (2.7) ^ [aP I

7t lr'(p^Mz^'-zf]-lr'{pS)Q'^' log(ef -z/) . (2.8) a ap J The definitions of the terms appearing in eqs. (2.5)-(2.8) are given below. The

superscripts on the left hand side of eqs. (2.5)-(2.8) identify the location of the source and

field points, respectively, as being in material A or B. For example, a 12 superscript

indicates that the source point (location of the point load) is in the upper half plane

(material A) and the field point (calculation point) is in the lower half plane (material B).

The fundamental solutions for the boundary element analysis involve the transpose of the

Green's functions given in [1] due to the convention taken in the boundary element

literature where Uij represents a displacement in the j -direction due to a unit load in the i- direction. The convention taken in [1] for the displacement Green's function is the

reverse, that is ReG,= t/., . (2.9)

This becomes critical for bimaterial problems since the fundamental displacement

solution is not symmetric, U-j ^ Uj..

The terms appearing in eqs. (2.5)-(2.8) are defined through the roots of the sextic

equation

detA(q)= 0, (2.10)

where the Christoffel matrix A,^ is defined by = c„j,q,q, . (2.11)

The wave vector q has components qj and ^2 ^i^d the roots are obtained from eq.

(2.10) such that ^2 = Pa^i- elastically stable solids the roots are complex [1,2,

12] and the roots are labeled such that

lmp^> 0. (2.12)

The complex coordinates appearing in eqs.(2.5)-(2.8) are defined as Za'- (2.13)

Z'a''= K+Pa''^2 (2.14)

31
The remaining matrices in eqs. (2.5)-(2.8) are defined as follows:

^^1 [Pa -Pa )ii[Pa ' Pp j[Pa - Pp

j r= Cofactor of A (2.16)

^(Pr)= L(pr)rK') (2.17)

^4?^")= '^i2k\+ Pa^'^nn (2.18)

Q'p=M{a'{p',)-a:yrV{p;)} (2.19)

' Q^= N{<T*(p^)-a^rr'r*(p;)} (2.20)

:

: - QTM = m{ct(p«)- alVr'rfp^)} (2.21)

Q^''=N{cT(p«)-aM"VW)} (2.22)

; M=rr'{<Trrr'-<^?rr'}"' (2.23)

- N=rr'{<T7y;«-'-<T^yry' (2.24)

.

' \ rr = Sy(pr) (2.25)

^ 'yf-' = l<y{pi"] (2.26)

In eqs. (2.15)-(2.26), the index a generally takes values of 1 to 3. There are several

useful relations for checking the computation of the various matrices given in [13].

These relations are most helpful in checking the Green's function computation in the

boundary element code.

The traction fundamental solution may now be computed from the displacement

fundamental solution by substituting eqs. (2.5)-(2.8) into

t.= {cijki+ Pa^Cij,2)-p^nj , (2.27)

where nj is the j-th component of the local normal vector. Here, we consider the case of

cubic material symmetry where the stiffness tensor is defined in terms of three elastic

constants Cn, Cj2 , C44 as [6]

Cyw= C44[5i,,5ji+ 5^idj,,)+ c^2^ijS,,i- H5ij5^i5ij^ , (2.28)

(no sum on k or i) where the anisotropy ratio H is defined as H= lc,,+ c,2- c„. (2.29)

In eq. (2.27) we have used the relation

. - ^«"y-= ^"'t^- (2.30)

dx^ dz„

32
Note in eqs. (2.5)-(2.8) that each of the equations is a function of z or z* only. Therefore,

we will use the notation of eq. (2.27) where it is understood that if u = u(z*), the

derivative is taken with respect to z*.

The traction fundamental solution can be written as T=-^E\ (2.31)

where

Pa^C^ 2«1 + C44«2 ^44^1 + p^'^C^ 1^2

(2.32)

The post-multiplication of the displacement derivative matrix arises from the

aforementioned convention taken in [1] for t/^ The matrix of derivatives in eq. (2.31) can

be determined directly from eqs. (2.5)-(2.8). The derivatives are given in Appendix A.

We note again that for the interface problem we have four traction fundamental solutions

to consider depending on the relative location of the source and field points.

A natural question to consider is the equivalence of the displacement and traction

fundamental solutions given above with the Kelvin solution (see, for example, [5]) for the

degenerate isotropic, homogeneous case. The homogeneous case can be modeled simply

by specifying identical elastic constants for materials A and B. The isotropic case must

be considered through a limiting process for the anisotropic solution given here; however,

we can consider a "near" isotropic, homogeneous case and compare the Tewary

displacement and traction solutions with the Kelvin solution. We consider the case of

tungsten where the anisotropy ratio H = 0. From [6] the elastic constants are = 521.0

Gpa, = 201.0 Gpa, and C44 = 160.0 Gpa. These constants give us a Poisson ration v =

0.28 and a shear modulus [L = 160.0 Gpa. For the anisotropic fundamental solution we

use = 521.0 Gpa, c,2 = 201.1 Gpa, and C44 = 160 Gpa. These values provide an

anisotropy ratio of H = 0.001

The displacement fundamental solution for both the Kelvin and Tewary solutions

is shown in Fig. 2. The solution is for a homogeneous solid so C/21 = t^ii the figure.

Note that there is essentially no difference in the U21, U^2^ and U22 components; however,

the 1 components differ by a constant. This constant depends on the elastic constants

being used but is spatially invariant. In the boundary integral equations such a constant can be treated separately in the contour integrals and can be shown to be inconsequential

to the solution of the equations. Indeed, a constant in the displacement fundamental

solution must not have any effect since this is simply related to rigid body motion of the

solid.

The traction fundamental solution for the Kelvin and Tewary solution is shown in

Fig. 3. The agreement is excellent with no discernible difference between the two

solutions. For purposes of comparison, the traction fundamental solution for two cubic

crystals (copper-nickel) is shown in Fig. 4. Note the difference in the solutions as the

field point crosses the bimaterial interface and that t/21 ^12 riot equal as they are

for the homogeneous case.

33
3. Boundary Integral Equations

We now consider the use of the anisotropic fundamental solution for the interface

problem shown in Fig. 1. For the region Q. in the figure, let the boundary of material A

be ^^, the boundary of material B be and + r2 = dQ. The standard boundary

integral equation written for the region is

, . CijUj{P)+ \T,j(P,Q)uj(Q)dr= \u,j{P,Q)tj(Q)dr, (3.1)

where P and Q are points on the boundary at the source and field points, respectively. As

noted above, we must consider the relative locations of the field and source points in

writing the kernels of the integrals. Following [10, 11], we can write eq. (3.1) for the source point P in material A as

CijUj(PO+ \T^\P^,QOuj(QOdry+ jT^\Pi,Q2)uj(Q2)dr2 = n r2

jul/(Pi,QOtj(QOdry+ jujj^(Pi,Q2)tj{Q2)dr2, (3.2)

ri ri

where the superscripts on the fundamental solutions terms and the subscripts on P and Q

identify the location of the source and field points. We can write a similar equation to eq.

(3.2) for the case when P is in material B, or we can write in general

CijUj(Pp)+ \T^f\Pp,QOuj(QOdry + jT.f\Pp,Q2)uj(Q2)dr2 =

\U,^\Pp,Qi)tj{Qi)dri+ ju^^\Pp,Q2)tj(Q2)dr2, (3.3)

n r2 where |3 = 1, 2 and no summation on (3 is implied. The continuity of displacement and

traction fields across the interface is insured through the use of the fundamental solution

from [1]. Discretizing each boundary into elements eq. (3.3) becomes

^ij^j l\Ti\Pp,Ql )Uj (21 )^l + i 1 iPp,Q2 )Uj {Q2 )^2 =

I.\U^\Pl5^Q0tjm^i + l\U^\Pp,Q2)tj(Q2)dr2 . (3.4)

. n=lri n=ir2

At this point we assume that the boundary traction and displacement are constant over each individual boundary element to allow us to analytically evaluate the integrals.

Factoring out the constants from the integrals in eq. (3.4) we have

^ij^j(Pp)+ I.^jiQi)lT,f\Pp,Q0dr, + '^UjiQ2)\Tif\Pp,Q2)dr2 =

n=\ n

n=l r2

'ltj(Q0\ut^Pp^Q0dri + f,tj(Q2)\U^\Pp,Q2)dr2 , (3.5)

n=l Ti n=\ r2 where the Uj and tj are element midpoint values. Note that eq. (3.5) is applied at 2N

points around dQ, N points each on Fj and F2. We therefore need to analytically evaluate

the following eight integrals:

34
lTP'[Pp,Q,)dr, (3.6) riTif^iPp^Qi)^^ (3-6) rlU^'[Pp,Q,)dr, (3.6) r\u^^[Pp,Q2]dr, (3.6) r where (3=1,2. To analytically evaluate these integrals, we use the mapping of Cruse [4] to transform these contour integrals to definite integrals in the complex plane. The

differential contour element dF is mapped as dr=—dz„, (3.7)

where

Ba= Pa^n,-h2, (3.8)

and and «2 unit vector components of the local normal vector for dF. With this

mapping the displacement kernel integrals simply involve evaluating integrals of the form

jlog(z-^)^z, (3.9)

and the traction kernel integrals involve evaluating

J A • (3.10)

z-b

Zl where is a complex constant. These integrals are simple to evaluate; however, special care must be taken when integrating across the branch cut of the principal value logarithm

function on the negative real axis. The location of the real axis does not necessarily

coincide with the location of the physical interface because of the definition of the

complex variables used here (see eqs. (2.13)-(2.14)).

4. Internal Displacement and Stress Calculation

To calculate the displacement at an internal point we use the standard form of

Somigliana's identity,

h{p)= \u^i{p,Qyi{Q)dr- \T,i{p,Q)ui{Q)dr , (4.1)

dn do. and the corresponding form for stresses,

aij{p)= \D,ij{p,Q)t,iQ)dr- \ S,,jip,Q)u,iQ)dr , (4.2)

where p is the internal computation point and Q is on the boundary. The kernels

appearing in eq. (4.2) are determined from the fundamental displacement and traction

solutions and the anisotropic constitutive law,

(y:j= {<=ijkl+ Pa^^ijki)^ (4-3)

35
We then have,

(4.4)

"a

Siij = [cijki + Pa^Cijki)-^ (4.5)

Recall that for a cubic solid the elastic constants are given by eq. (2.28). The derivatives

appearing in eq. (4.4) have been computed previously in deriving the fundamental

traction solution, eq. (2.31), and are given in Appendix A. The derivatives appearing in eq. (4.5) can be written as -1—- , 2 ^ml (3.16)

dZa dz„

Therefore, we only need the second derivatives of the fundamental displacement solution.

These are summarized in Appendix A. Finally, we note again that the E" matrix

appearing in eq. (4.6) depends on material symmetry. We focus here on cubic material

symmetry but others are easily implemented as well.

5. Example Problem

To demonstrate the application of the method we analyze the field distribution in a copper-nickel multilayer. Multilayer materials are fabricated by depositing alternating

layers of thin-film materials such as Cu-Ni, Co-Cr, and Fe-GaAs. Here, we demonstrate

the usefulness of the analysis method detailed above in analyzing multilayers. For the

example problem we use a portion of a Cu-Ni multilayer subjected to mechanical loading.

Fig. 6. The specimen shown in the figure is under shear loading with traction free

surfaces along the remaining boundaries. We take the elastic constants as [6]:

^11 ^12 <^44 Cu 168.4 121.4 75.4

Ni 246.5 147.3 124.7

where all values are in GPa.

The elastic displacements along x = 0 are shown in Fig. 6. The displacements

exhibit qualitatively what we expect from the shear loading. Note the continuity in the

displacements as the interface is crossed and the lack of perfect asymmetry in and M2

due to the differing elastic constants in the two materials.

6. Summary

We have presented an efficient, accurate method for analyzing deformation near

interfaces in anisotropic solids. The use of a special Green's function for the anisotropic

interface problem allows us to discretize only the boundary of the problem. The interface

itself does not have to be discretized since its behavior is explicitly incorporated into the

Green's function. Embedding the Green's function in a boundary element approach then

provides us with a general tool for analyzing a variety of bimaterial problems.

A large portion of the computation time for a complete analysis is used by the

computation of the Green's function. The computation involves calculating roots of the

Stroh determinant for use in the matrix functions of the Green's function; however, once

36
calculated the matrices may be stored for subsequent analysis of the same material

system.

Acknowledgments

The author acknowledges the support received from the National Academy of

Sciences under the National Research Counsel Postdoctoral Associateships Program at the National Institute of Standards and Technology. The assistance and encouragement

of Dr. Vinod Tewary at NIST is gratefully acknowledged.

References

[1] Tewary, V.K., Wagoner, R.H., and Hirth, J.P., "Elastic Green's Function for a Composite Solid with a Planar Interface," Journal of Materials Research, Vol. 4, p.l 13-

123 (1989).

[2] Stroh, A.N., "Steady State Problems in Anisotropic Elasticity," Journal of

Mathematical Physics, Vol. 41, p. 77-103 (1962).

[31 Lekhnitskii, S.G., Theory of Elasticity of an Anisotropic Elastic Body, Holden-Day, Inc.,

(1963).

[4] Cruse, T. A., Boundary Element Analysis in Computational Fracture Mechanics, Kluwer

Academic Publishers (1988).

[5] Brebbia, C. A. and Dominguez, J., Boundary Elements, An Introductory Course,

Computational Mechanics Publications/McGraw-Hill Book Company (1989).

[61 Hirth, J. P. and Lothe, J., Theory of Dislocations, Second Edition, Krieger Publishing Co.

(1992).

[7] Cruse, T. A., Mathematical Foundations of the Boundary-Integral Equation Method in Solid

Mechanics, Air Force Office of Scientific Research Report, AFOSR-TR-77-1002 (1977).

[81 Snyder, M. D. and Cruse, T. A., Crack Tip Stress Intensity Factors in Finite Anisotropic

Plates, Air Force Material Laboratory Technical Report, AFML-TR-73-209 (1973).

[9] Snyder, M. D. and Cruse, T. A., "Boundary Integral Equation Analysis of Cracked

Anisotropic Plates," International Journal ofFracture, Vol. 11, No. 2, pp. 315-328 (1975).

[10] Yuuki, R., Cho, S.-B., Matsumoto, T., "Usefulness of Hetenyi's Solution for Boundary

Element Analysis of Crack Problems in Dissimilar Materials," Role of Fracture Mechanics in

Modem Technology, G. C. Sih, H. Nisitani, and T. Ishihara (Eds.), Elsevier Science Publishers,

pp. 823-834 (1987).

37
[11] Yuuki, R. and Cho, S.-B., "Efficient Boundary Element Analysis of Stress Intensity Factors

for Interface Cracks in Dissimilar Materials," Engineering Fracture Mechanics, Vol. 34, No. 1, pp. 179-188 (1989).

[12] Eshelby, E., Read, W. T. , and Shockley, W., "Anisotropic Elasticity with Applications to Dislocation Theory," Acta Metallurgica, Vol. 1, pp. 251-259 (1953).

[13] Tewary, V. K., Wagoner, R. H., and Hirth, J. P., "Elastic Green's Function for a Composite

Solid with a Planar Crack in the Interface," Journal ofMaterials Research, Vol. 4, No. 1 (1989).

[14] Benzley, S. P., "Representation of Singularities with Isoparametric Finite Elements,"

IntemationalJoumalfor Numerical Methods in Engineering, Vol. 8, pp. 537-545 (1974).

[15] Gifford, L. N., and Hilton, P. D., "Stress Intensity Factors by Enriched Finite Elements,"

Engineering Fracture Mechanics, Vol. 10, pp. 485-496 (1978).

[16] Heyliger, P. R., and Kriz, R. D., "Stress Intensity Factors by Enriched Mixed Finite

Elements", International Journal for Numerical Methods in Engineering, Vol. 28, pp. 1461-1473

(1989).

38
Appendix A - Derivatives of the Fundamental Displacement Solution

Case (i): 0, X2' > 0

dU11

dZr = Re-^

K

1^ A 'A ^

_A _'A* a ap

dz; = Re —<! 7t

T 1

(A.l)

(A.2)

Case (ii): X2< 0, X2' > 0

i5e = Rei

dza 7t

d^V^^ _ _ 1 *9— — RS —

dZa ^

ap Za -Zp

^p{zT-zff

(A.3)

(A.4)

Case (iii): X2> 0, X2' < 0

21

dZa

dZr K

{z'a-Zp'f

J

(A.5)

(A.4)

Case (iii): X2< 0, X2' < 0

dU^^ „ 1 — = Re —

dZa n

„ 1 *9 =Re—<! dza ^

a Za Z(x aP Zp

(A.5)

(A.6)

39
r

1

Material A

Material B

Fig. 1 General bimaterial problem

0.005

-0.005

-0.015

-0.025

-0.035

UI2 = U2I

Ul I (Kelvin Solution)

-100.0 -50.0 0 0 50 0

x-coordinatc (arbitrary units)

100.0

Fig. 2 Displacement fundamental solution comparison for an isotropic, homogeneous material.

40
0030

-0.010

-100 0 -50.0 0 0 50.0

x-coordinatc (arbitrary units)

100.0

3 Traction fundamental solution comparison for an isotropic, homogeneous material.

Fig. 4 Displacement fundamental solution for a Cu-Ni system.

41
Cu

Ni

60 mm

Fig. 5 Geometry of the Cu-Ni multilayer.

5.00

K-coordinatc (nun)

Fig. 6 Displacements along x = 0 for the Cu-Ni system.

42
Green's Functions in Elastic Fracture

Mechanics

T. A. Cruse^

1. Selected 2D Green's Function Defined

The current brief note is intended to provide only a summary of the major elements of a two

dimensional Green's function that has turned out over a number of years to have been of great

practical utility. The overall problem is that of two dimensional fracture mechanics analysis. The

sub-problems that have been addressed in the past include elastic-homogeneous, elastic- inhomogeneous, and elastic-plastic fracture mechanics.

The first sub-problem is the simplest problem of an elastic crack subject only to some state of

remote loading. The second concerns the interaction of the elastic crack with inhomogeneous strain

fields due to residual, thermal, or prior-plastic

strains. The last is the interaction of the crack

with the material under locally elastic-plastic

material behavior. The current note will summarize only the first two as they are of most

interest to the author in the context of advanced

material behavior.

The selected Green's function formulation

is that of the straight, stress-free crack in an

elastic, fiilly-anisotropic material, as shown in Figure 1 . The crack is of length 2a which is taken for convenience to be aligned with the xaxis. The crack surface is taken to be /"and is

stress-fi"ee. The Green's function which has

been formulated is that for point loads in the

infinite body with the crack geometry shown in Figure 1 . The point load is to be taken at an

arbitrary point c(x,y) which may also be located

anywhere on the crack.

The two dimensional elastic formulation is based on a complex variable representation ofthe

material and geometry, using the full anisotropic formulation of Lekhnitskii [1]. The details of the

original anisotropic boundary-integral equation are contained in the report by Cruse and Swedlow

[2]. The formulation of the Green's function, obtained formally as a Hilbert problem solution is given in Snyder and Cruse [3, 4].

S

Figure 1: Geometryfor elastic Green'sfunction

' H. Fort Flowers Professor of Mechanical Engineering, Vanderbilt University, Nashville, TN 37235.

43
A key term in the complex variable solution for the elastic fracture mechanics Green's function

of Snyder and Cruse was later determined by Cruse [5] to be missing an important constant of

integration. The constant of integration was required to keep all of the complex variable jump

terms on a consistent basis as the source point is moved throughout the anisotropic domain. The

key term has the form given in the following equation

V^^ - a ^ Vc^ - + cz - a"^ J (z ,c) - n - In

In this equation, the terms z, c refer to the complex variable mapping location for the

integration point and the point load location, respectively; i.e., z(x,y) = z + juy, where // is the

complex root of the material characteristic equation [3]. The Jfz.cj term contains square-root

factors which contain the analytical singularity associated with the unique and singular crack tip

terms for both crack tips {z = ±a). As. we shall see later, these square roots lead to the analytical

BIE expression for the crack tip stress intensity factors. Additionally, the form ofJ(z,c) contains

the critical mapping terms that differentiate applying the point load to the top and the bottom of the

crack. This necessary discontinuous behavior is also contained in some of the above square-root

terms.

The critical "constant" of integration is given by the term for the point load location in the denominator for Eq. (1). In the absence of

this term, the mapping function represented by

J(z, c) has the shape indicated by the computer

representation given in Figure 2 for a selected,

arbitrary position of the point load. It can be

seen that the mapping of the crack surface

results in repeated crossings of the negative real axis. Each of these crossings introduces a

spurious jump in the Green's function result. As

discussed in [5], the complex variable constant of integration which depends on the location

c(x,y) given in Eq. (1) rotates the mapping in Figure 2 such that there are no crossings of the

negative real axis by the contour of the mapped

crack.

The formulated Green's function can then be used in the numerical BIE solution of arbitrary

finite body loading problems, represented in Figure I. The resulting BIE is given as follows:

y» 0

CjMP) + I

t;{P ,QyXQ)dS + \ mP,Q)u,{Q)dS

y 0

44
The terms in the above BIE are based on reasonably standard notation. The points P,Q

represent the physical locations of the point load and field points, respectively. The symbols U, T

and u,t represent displacement and traction terms for the Green's fiinction and physical boundary

conditions, respectively.

The superscript-star notation indicates that the fundamental solutions contain the necessary

Green's function terms of the stated fracture mechanics model. Based on the stated stress-free conditions on the crack surface and the properties of the Green's function (also stress free on the

crack), the last integral in each of the above equation lines is zero, for any location of the point

load. We see, therefore, the essential property of a true Green's function in that the critical boundary feature (in this case, the crack) is accounted for by the kernel functions of the

fundamental solution, and not by explicit boundary integrals. As usual, the BIE is a formal

constraint equation for the boundary data which assures complete equilibrium is assured, only now

the constraint equation is for the uncracked boundary only.

There is a second, critical feature of the elastic Green's function formulation for the fracture

mechanics problem. Given the properties of the Green's function terms, represented by Eq. (1), the

elastic singularity of the crack tip stress and strain fields is exactly and completely contained in the

associated Somigliana identities. As outlined by Snyder and Cruse [3], the limiting form of the

interior Somigliana stress can by taken to obtain the stress intensity factors (SIF) of the elastic fracture mechanics problem. The resulting SIF results are represented by the following:

{K„K„) = 1

4' ") {Q)t, {Q)dS - 1 R\'-"^ {Q)u{Q)dS p,

The result above shows that the elastic SIF's in the fracture mechanics solutions are given

solely in terms of two vector-product integrals of the boundary conditions on the uncracked surface

only. Further, as discussed by Cruse and Wilson [6], the above integral identity is pathindependent, requiring only that tractions and displacements on the crack be known. They also

showed that the integral could be used as a post-processing routing for normal finite element

models of cracked bodies which gives much greater accuracy than the normal crack tip

approximation methods.

The resulting Green's function formulation and BIE implementation was shown by Cruse [5] to be extremely accurate, for very small boundary meshes in the context of boundary element

analysis. The method has been used by a variety of industries, agencies, and consultants for practical fracture mechanics solutions of technologically significant problems. A few applications are cited by Cruse [7]. The anisotropic elastic fracture mechanics capability was originally

developed for composite material problems, but the code was also used to model crack growth in nickel-base single crystals by Chan and Cruse [8].

2. Some Extended Applications

While the original Green's function derivation was for the elasticity problem, the formalism was easily extended to the elastic-plastic fracture mechanics problem in a series of papers

[9,10,1 1,12]. For the present purposes of looking at advanced materials, there is one key feature of

the elastic-plastic formulation that is a basic Green's function result. The path independent integral

result for the elastic stress intensity factor is modified in the presence of any domain distribution of

45
inhomogeneous strains These strains may arise from residual strains, such as in welding, or

plastic strains due to mechanical loading of the uncracked geometry. These applications are discussed in Ref [11].

The path independent integral for inhomogeneous strains is given by the following relation,

taken from [9]:

The Green's function plays the same role of eliminating the crack boundary from the equilibrium

formulation (BIE) and from the SIF calculation, shown above. Again, the inhomogeneous strains in Eq. (4) can derive from any numerical solution method. In the case of exact inhomogeneous strains,

the result above has been shown to give numerically-exact SIF solutions [9].

One of the important advantages of analytical Green's function formulations is the ability to extract analytical results of critical field variables, such as the SIF results that have been cited. The

elastic-plastic Green's function formulation also was used to confirm the limiting strain singularity

for perfectly-plastic material [9]. It appears to this writer that an eigenvalue approach to the

internal strain singularity power determination could be developed to confirm the singularity for

non-perfectly plastic material behavior. Such a theoretical development would be important to derive, as it would lay the foundation for determining singularity strengths for more complex

material models in the future.

3. Limitations to the Green's Function IVIethod

The principal "failing" of the Green's function method is its usual reliance on complex variable

mathematics, needed to obtain analytical results for complicated geometries. The behavior of the

Green's function terms alluded to in Figure 2 is intended to convey this additional analytical

complexity. Also, there are no analytical Green's function results for true, three dimensional crack

problems. The lack of solutions is partly due to the fact that 3D cracks have much more

complicated geometries: not only are they possible non-planer, but the crack front curvature

introduces significant three dimensionality to the problem.

Numerical fracture mechanics results using the BIE formalism have resulted in some analytical

results for 3D fracture mechanics. One theoretical study which is only partially complete was an

effort to confirm the order of stress singularity for the crack-free surface intersection problem [13]. In this case, the interior stress Somigliana identity was used to investigate the singularity order

through an expansion of the integral equation terms. The result confirmed the earlier analytical

solution by Bentham [14] based on a polar decomposition and series expansion of the local stress

The principal finding of this earlier study is that the three dimensional numerical BIE does

have the ability to resolve at least a 2D portion of the full 3D problem. However, it does not seem

likely that analytical results for the variation of the crack tip SIF around the crack front will ever

(4.)

field.

46
be found. The closest to an analytical solution is that of the elliptical crack under polynomial

surface loading, as developed by Nishioka and Atluri [15].

It is likely that researchers will, instead, seek to develop numerical Green's functions for 3D

problems. The recent literature strongly suggests that these formulations will be based on nonsingular forms of the traction-BIE, in one of its various forms. One example of that is the recent work which shows the lack of strong (Cauchy) singularities in the traction-BIE for fracture

mechanics modeling [16].

4. Advanced Material Modeling Opportunities

There are several advanced material modeling problems that seem likely to greatly benefit from

the further development and numerical implementation of 2D Green's functions. These will now be

briefly listed, without implying any specific priority ranking.

1. Material Inhomogeneity Effects: The Green's function approach has been

successfully applied to a reference problem of a crack approaching a material interface

[6]. It seems clear that the physical singularity in the inhomogeneous fracture problem

must be that of the homogeneous field so long as the crack is not overly small. Thus, a

first effort might be to include the material inhomogeneity of distinct material grains

through a domain element representation of the inhomogeneities.

2. Distinct Reinforcements: The BIE formulation can be modified to account for

discrete reinforcements through an internal compatibility matrix approach [17]. The

linkage to discrete reinforcements is through a force-displacement model of the

reinforcement, together with additional log-terms in the BIE formulation to represent

discrete point loads at specified points in the domain. This approach could be used

effectively to model crack bridging in metal matrix composites materials.

3. Cracking in Limited-Ductility Materials: Considerable attention is being given to attempting to use advanced engineering materials whose behavior might be thought of

as brittle, except for the addition of reinforcing or crack diffusing particles and fibers. One major characteristic of these materials is the presence of multiple small cracks

which are interacting with the reinforcements. A first step would be to develop

effective Green's function formulations for multiple cracks, and then to link these with

discrete inhomogeneity models.

4. Functionally Graded Materials: Considerable attention is being given in heat engine

applications to the use of ceramic coatings to reduce thermal fluxes and metallic

temperatures (a macro-example is the space shuttle surface which uses discrete

ceramic tiles). Cracks have been found to occur in ceramic thermal barrier coatings

which are parallel to and near material interfaces [18]. Green's functions which

capture the coupled flux and mechanical loading interactions are needed for cracks at and near material interfaces. The effect of interface roughness is also an issue.

5. Reliability-Based Material Design: The BIE formalism lends itself very effectively to

sensitivity analysis of boundary related features. Advanced reliability algorithms are based on the use of sensitivity algorithms (e.g., [19]). The design of advanced

47
materials for their minimum reliability instead of laboratory strength will require the

application of robustness concepts to the material configuration. The BIE Green's

fiinction methods for sensitivity analysis may be particularly well suited for use in

reliability simulation of advanced material design.

1 S . G. Lekhnitskii, Theory of Elasticity of an Anisotropic Body, Mir Publishers, Moscow 1981.

2 T. A. Cruse and J. L. Swedlow, Interactive Program for Analysis and Design Problems in Advanced Composites Technology, United States Air Force Technical Report AFML-TR-71-

268, 1971.

3 M. D. Snyder and T. A. Cruse, Crack Tip Stress Intensity Factors in Finite Anisotropic Plates,

Air Force Report AFML-TR-73-209, 1973.

4 M. D. Snyder and T. A. Cruse, Boundary-Integral Equation Analysis of Anisotropic Cracked

Plates, International Journal of Fracture, 11, 315-328, 1975.

5 T. A. Cruse, Two Dimensional Boundary-Integral Equation Fracture Mechanics Analysis,

Applied Mathematical Modeling, 2, 287-293, 1978.

6 T. A. Cruse and R. B. Wilson, Advanced Applications of Boundary-Integral Equation

Methods, Nuclear Engineering and Design, 46, 223-234 (1978).

7 T. A. Cruse, Design Applications of BIE Methods, SAE Transactions, 93, Section 3, Society

of Automotive Engineers, Warrendale PA, 3 .1078-3 .1084, 1985.

8 K. S. Chan and T. A. Cruse, Stress Intensity Factors for Anisotropic Compact-Tension

Specimens with Inclined Cracks, Engineering Fracture Mechanics, 23, 5, 863-874 (1986).

9 T. A. Cruse and E. Z. Polch, Elastoplastic BIE Analysis of Cracked Plates and Related

Problems, Part 1 : Formulation, International Journal for Numerical Methods in Engineering, 23, 429-437 (1986).

10 T. A. Cruse and E. Z. Polch, Elastoplastic BIE Analysis of Cracked Plates and Related

Problems, Part 2: Numerical Results, International Journal for Numerical Methods in

Engineering, 23, 439-452 (1986).

11 T. A. Cruse and E. Z. Polch, Application of an Elastoplastic Boundary-Element Method to Some Fracture-Mechanics Problems, Engineering Fracture Mechanics, 23, 6, 1085-1096

(1986).

12 T. A. Cruse and S. T. Raveendra, A Comparison of Long and Short Crack Elastoplastic

Response using the Boundary Element Method, Engineering Fracture Mechanics, 30, I, 59-

75 (1988).

13 T. A. Cruse, Three Dimensional Elastic Surface Cracks, in Fracture mechanics: Nineteenth

Symposium, ASTM STP 969, T. A. Cruse, Editor, American Society for Testing and

Materials, Philadelphia, 19-42 (1988).

14 J. P. Bentham, State of Stress at the Vertex of a Quarter-Infinite Crack in a Half Space,

International Journal of Solids and Structures, 13, 479-492 (1977).

48
15 T. Nishioka and S.N. Atluri, Analysis of Surface Flaw in Pressure Vessels by a New Three

Dimensional Alternating Method, Journal of Pressure Vessel Technology, 104, 299-307

(1982).

16 Q. Huang and T. A. Cruse, On the Non-Singular Traction-BIE in Elasticity, International

Journal for Numerical Methods in Engineering, 37, 2041-2072 (1994).

17 Fernando De Falco, Metodo degli elementi di contomo per problemi di interazione terrenostruttura in campo elastico non lineare. University of Rome, 1992.

18 T. A. Cruse, S. E. Stewart, and M. Ortiz, Thermal Barrier Coating Life Prediction Model

Development, Journal of Engineering for Gas and Turbine Power, 112, 251-260 (1990).

19 T. A. Cruse, K. R. Rajagopal, and, J. B. Dias, Probabilistic Structural Analysis Methodology

and Applications to Advanced Space Propulsion System Components, Mathematical

Computer Modeling, 15, 3-5, 103-118, 1991.

49

Mechanical Properties of Metal-Matrix Composites

L. C. Davis, Research Laboratory

Ford Motor Co.

MD 3028 SRL, Dearborn, MI 48121-2053

ABSTRACT

Metal-matrix composites (MMC) show great potential as light weight materials for a variety

of automotive applications. The high specific modulus and strength of these materials make

them suitable for disk brake rotors, connecting rods, cylinder liners, and other high

temperature components. Typically, MMC's consist of an aluminum matrix reinforced by

particles of SiC or other ceramics. The volume fraction of the reinforcing phase is in the

range 15 - 30% and the particle size is generally 1 - 10 /xm. Considerable effort has been

put into developing models that accurately describe the properties of MMC's. The inputs

to the models are the properties of the constituents (e.g., moduli, yield strength, coefficient

of thermal expansion, creep rate), the volume fraction of particles, and information

regarding the phase geometry (particle shape and distribution). Unlike fiber-reinforced

composites where fiber-matrix interface properties are often dominant, the metal-ceramic

interface can usually be assumed to be ideal. Some linear properties (for example modulus)

can be estimated with sufficient accuracy by applying rigorous bounds, either HashinShtrikman or more accurately the third-order bounds of Milton, Torquato, and others.

However, for non-linear properties such as yield strength, finite element analysis (FEA) has

proven to be extremely useful. FEA modeling at Ford Research Laboratory has included:

(1) residual stresses induced by thermal mismatch between the matrix and the reinforcing

particles, (2) stress-strain relation under uniaxial loading, and (3) creep of MMC's. In most

cases, the FEA results agree well with experiments done by Allison and coworkers in the

Materials Science Department of our laboratory. An interesting example of the verification

of the modeling has been the measurement by neutron diffraction (with Los Alamos

National Laboratory) of internal particle strains due to residual and applied stresses. On

the other hand, comparing FEA calculations with measurements of creep rates has not

shown as good agreement, possibly demonstrating the significance of changes in matrix

microstructure in the composite relative to the unreinforced material. One of the current

limitations of the FEA approach is the restriction to unit cell models. It has not been

possible to analyze realistic particle arrays (e.g., random) in three dimensions, although some

work in two dimensions has been reported.

L EraiODUCriON

Particulate-reinforced metals, or metal-matrix composites (MMC) as they often

referred to, are being considered for several automotive apphcations [1]. These applications

take advantage of the reduced weight characteristic of the matrix, which is almost always

aluminum, and the increased modulus and yield strength imparted by the ceramic (SiC, TiC)

reinforcements. For many components (e.g., connecting rods and disk brake rotors), the

volume fraction of inclusions is in the range 15-30%. Particle sizes are generally 1 -10 /xm.

Considerable effort has been put into developing models to understand the properties of the

materials and to provide a basis for design rules. Modeling of MMC's is done at two levels.

Analytical expressions for bounds on the linear properties (e.g., moduli) have been

computed. Since the upper and lower bounds are close together, these are sufficient for

estimating properties in many cases [2]. For more accurate results and for all non-linear

properties, such as yield stress and creep rates [3], finite element analysis (FEA) of the

micromechanics of the material is required.

51
The purpose of this paper is to present several examples of FEA modeling of the

mechanical properties of MMC's. Most of the FEA was done with the commercial code

ABAQUS, most recently running on a Cray Y-MP2E. Examples include the bulk modulus

as a function of the volume fraction of inclusions, residual stresses due to thermal expansion

mismatch, stress-strain relationships with metal plasticity, and creep of composites. In all 3D cases, the modehng is restricted to unit-cell models where the particles are of uniform

size and shape (e.g., sphere, cylinder). In 2D, random structures can be studied [4,5,6].

The inputs to these models are the properties of the constituents (e.g., moduh, yield

strength, coefficient of thermal expansion, creep rate, the volume fraction of inclusions, and

phase geometry information such as particle shape and arrangement). Unhke fiberreinforced composites where the fiber-matrix interface characteristics can be dominant, the

metal-ceramic interface can be treated as ideal (perfectly bonded) in most cases.

n. ELASTIC MODULI

The simplest composite mechanical property to calculate is the bulk modulus K. In

Table I we give the properties of the matrix (Al) and of the inclusions (SiC,TiC). The

Young's modulus E for the ceramics is significantly larger (by a factor of 6) than that of Al.

Thus we expect a substantial increase in K as the amount of ceramic included increases. In

Fig. 1, we show the geometric model used for this calculation. It consists of spheres on a simple cubic lattice. Since the structure is periodic, we can restrict the analysis to the unit

cell. Furthermore, by applying symmetry, we reduce the computational volume to the

wedge-shaped region in Fig. 1. The boundary conditions on the wedge are:

' uSi)=U^ (y=0\ u(i)=Q (y=L\ u^(i)=U^ (x=L)

(1)

w.(0=0 (z=0), u(i)=u(i) (diagonal surface)

A load (external force) can be applied to the node represented by Uy giving a uniform

(across the front surface) apphed stress Oyy^. No force is applied to the node represented

by U^, which is left as a free variable to be determined by FEA . Poisson's ratio is determined according to

(2)

This calculation gives Young's modulus for the composite:

• • yy

where L is the length of wedge (2L is the lattice constant). K can be determined from the

expression

52
~3(l-v)

(4)

The calculated bulk modulus K is displayed in Fig. 2 as a function of volume fraction, 0.

The bulk modulus of the composite increases with 0 up to a critical composition at which

percolation occurs (0^, = 0.52). The rate of increase is, however, much less than a simple

rule of mixtures (hnear interpolation between the end points) estimate. The solid and

dashed lines in Fig. 2 are the rigorous bounds of Milton [7] calculated from the work of

Torquato et al. [8] It is reassuring that the FEA results fell between the upper and lower

bounds as required. These bounds are third order bounds that are specific to the simple

cubic lattice. Expressions for other lattices and random arrays, which are more appropriate

for comparing to experiment, also exist. (See Ref. 2.)

m. RESIDUAL STRESSES

Residual stresses can develop in MMC's due to the difference in coefficients of

thermal expansion (CTE) of the matrix and the inclusion [9]. The CTE's of ceramics are

1/6 to 1/3 of the Al CTE . If, during fabrication, the MMC is stress free at an elevated

temperature (aging temperature, for example), cooling to room temperature produces

significant compressive stresses in the particles. In linear elastic materials, it can be shown

that the pressure developed in an isolated particle is given by

m p

where (oi^) is the CTE of the matrix (particle), AT is the temperature drop, Gj^ is the

shear modulus of the matrix and Kp is the bulk modulus of the particle.

For a typical temperature change (AT = 256 °C), the stress in the matrix near the

particle becomes large enough for the metal matrix to flow plastically. Using a von Mises

stress criterion, we take the matrix to yield when the stress invariant

(where is the deviatoric stress) exceeds the yield stress oq (taken to be 245 MPa). Here

we are assuming the matrix to be elastic/perfectly plastic. Arbitrary strain hardening can be

accommodated by ABAQUS and measured curves on unreinforced material have been used

in other calculations that were compared to experiments.

A large plastic zone develops around the particle as shown in Fig. 3. The pressure

p in the particle is now less ( 284 MPa compared to 357 MPa) than that given by Eq. (3)

since Al has yielded in the region next to the particle . The pressure in the particle also

decreases with increasing volume fraction of inclusion. A discussion of this effect has been

given by Bullough and Davis [10].

53
The plastic zone influences any subsequent responses to external loading of the

composite. For example, if one wishes to determine Young's modulus for this material, it is difficult since there is no linear elastic region, no matter how small the applied stress. A

calculation of the stress-strain relationship for 30%TiC/Al is shown in Fig. 4. For tension,

there is a reasonable quasi-linear portion from which an approximate E can be found, but

for compression the curve never appears linear. In any event, the inferred Young's modulus

(slope of the curve at zero strain) is less than the true value determined in the absence of

residual stresses (top curve). Various aspects of this approach have been verified by neutron

diffraction studies [9].

IV. CREEP

Another property of composites that is of interest to materials scientists and designers

is creep [3]. Metals, particularly at high temperatures, slowly deform (over times of hours

or days) under stresses comparable to the yield stress. This slow deformation, which is in

addition to the initial elastic and possibly plastic response to apphed stress, is known as creep. The creep rate for metals can be written as

where the exponent n can be large (5-15). This is the expression for steady state creep. An

additional factor, typically f"^, can be included to describe primary creep, which occurs just

after loading. The steady state creep rate of the MMC can also be expressed in the same

form, Eq. (7), with an identical n, but with a different coefficient A^. (if only stress concentration due to the high modulus inclusions is important). Measurements of the creep

rates of MMC's sometimes give different exponents n, probably because of changes in matrix

microstructure. If the microstructure does not change, however, it can be shown that A^,/A

depends only on the phase volume and geometry (e.g.., volume fraction of inclusions, particle

shape and arrangement, etc.), not material properties such as elastic moduh. An example

is shown in Fig. 5 in which the composite strain rate vs time is shown for two different values

of the particle Young's modulus. [These FEA calculations were done using an axisymmetric

unit cell method.] In the steady state or asymptotic region at large times, the two curves come together. Other tests with various parameter changes have verified this lack of

dependence of A^./A on material properties. In the initial region of time, the stress redistributes from its original pattern following loading to the final stress state that is characteristic of the steady state. Comparison of the calculated creep rates to experimental

rates measured on 15%TiC/Al composites shows only modest agreement, indicating that

changes in microstructure of the matrix relative to the unreinforced material (and thus

changes in the matrix creep properties) have occurred. It may be necessary to consider the

dislocation network on scale finer than those employed in the continuum mechanics model.

V. 2D RANDOM COMPOSITES

A limitation of the resuhs presented so far is that they pertain to an ordered structure or unit-cell method. In real materials, particles are generally not in ordered arrays, although

sometimes in aligned-fiber composites, the arrangement can be fairly regular. Some

simulations of random materials in two dimensions have been reported. Calculations in

three dimensions are too difficult computationally at present.

Thorpe and co-workers [4] have performed simulations of the elastic properties of

two dimensional composites. Their models consist of ordered and random arrangements of

54
circular inclusions (equivalent to aligned fibers). They have used a digital-image-based

method or spring-grid scheme that is equivalent to a finite element method. A liquid

algorithm is used to generate the random, non-overlapping configurations. Other FEA

simulations on 2D random composites have appeared in the literature recently.

The area modulus k for composites of rigid fibers embedded in an elastic matrix are shown as a function of matrix area fraction 0^ in Fig. 6. Three ordered arrays as well as random arrays are displayed. The modulus for an ordered array diverges as k « (0^ - <t>ic)'^^ near the percolation concentration 0^^,. The behavior for random

composites (either non-overlapping or overlapping) differs. A mean field theory given by

Davis et al. [11] provides an adequate description (See Fig. 7).

VI. CONCLUSIONS

The modeling work reviewed in this paper has pertained to the following properties

of metal-matrix composites:

o Elastic moduli

o Residual stresses and their effects on stress-strain relationships

o Creep

o Plasticity and yielding phenomena in the matrix

Most of these properties are fairly well described by approaches based on rigorous bounds

and finite element analysis. Yet, limits on the accuracy of the results exist; the primary

limitations of the FEA modeling to date are:

o Restriction to unit-cell models for 3D

o Inability to treat properly microstructural and mesoscopic effects

In addition to the properties listed, a fundamental particulate-composite property that has

not been addressed is fatigue behavior. It is important to be able to model the effects on

fatigue of these parameters:

o Volume fraction of the ceramic phase

o Particle shape

o Particle arrangement (e. g., clustering) and size distribution

o Apphed stress 0 Microstructure of the matrix in the reinforced material

Some of these may be amenable to analysis by Green's functions and boundary elements,

while others may not. As a start, two-dimensional analyses on linear materials might provide

insight, but ultimately non-linear material properties and three dimensions must be treated

for a complete understanding. Finally, the role of interfaces has not been discussed. It is generally thought that in particulate systems interfaces are nearly ideal, although this

assumption probably needs to be re-examined.

ACKNOWLEDGMENTS

1 would hke to acknowledge my collaborators in this work: John Allsion, Ron

BuUough, Jidong Chen, Ken Hass, and Mike Thorpe.

55
Table I. Parameters used in finite element calculations.

Young's Modulus

h (CjFa)

72.0 (22 °C)

70.8 (150 C)

431.0

Poisson's ratio v 0.34 0.19

CTE 23.6 8.6

Creep rate amplitude A 2.79x10"^^

(150 °C) —

Creep rate stress exponent n

10.5

Temperature decrease

from stress-free state

AT(°C)

256 250

REFERENCES

1. J. E. Allison and G. S. Cole: JOM, 1993, vol. 45, pp. 19-24.

2. L. C. Davis, Metall. Tf-ans. A, 1991, vol. 22A, pp. 3065-3067

3. L. C. Davis and J. E. Allison, Metall. TYans. A, 1994, vol. 24A, (in press).

4. J. Chen, M. E Thorpe, and L. C. Davis, to be published; A. R. Day, K. A. Snyder,

E. J. Garboczi, and M. R Thorpe, /. Mech. Phys. Solids, 1992, vol. 40, 1031-1051.

5. S. Ghosh and S. N. Mukopadhyay, Computers and Structures, 1991. vol. 41, 245-256.

6. Y.-L. Shen, M. Finot, A. Needleman, S. Suresh, Acta Metall Mater., 1994, vol. 42, 77-

97.

7. G. W. Mihon, Phys. Rev. Lett. 1981, vol. 46, 542-545.

8. S. Torquato, Appl. Mech. Rev., 1991, vol. 44, 37-76.

9. L. C. Davis and J. E. Allison, Metall. Ti-ans. A, 1993, vol. 24A, pp. 2487-2495.

10. R. Bullough and L. C. Davis, Acta Metall. Mater., 1994.

11. L. C. Davis, K. C. Hass, J. Chen, and M. R Thorpe, Appl Mech. Rev., 1994, vol. 47,

S5-S9.

56
Simple Cubic Array

Finite Element Volume

o

Q-O

m

o2

13 m

150

125

100

75

50

25

1 1

A 1 /C? ^ Al/olU

—1

1 1

• FEA

- — Uppar Bound

Lower Bound

1

1 1 1 1 _. 0.0 0.1 0.2 0.3 0.4 0.5 Volume Fraction

0.6

Fig. 1 Unit cell for the simple cubic array. The volume

used for finite element analysis is 1/16 of the unit cell (the wedge-shaped region indicated).

Fig. 2 Calculated bulk modulus vs volume fraction of spherical SIC particles (E = 431 GPa, = 0.19) in an Al matrix (E = 72 GPa and y = 0.34) . The filled circles are from a finite element analysis, and the solid and dashed

lines are third-order bounds.

o

Q.

0)

10

200

50

150 -

100

r 30^ TiC/AI

Simple Cubic

-No residual - Tension

0 • Compression

0.000 0.001

Stroin

0.002

Fig. 3 Region of plastic flow in the Al matrix surrounding

the TiC particle due to cooling 256 K. The matrix is assumed to be elastic/perfectly plastic with tensile yield stress CTq = 245 MPa.

Fig. 4 Absolute value of the macroscopic stress vs strain

for tensile (filled circles) and compressive applied loads

(open circles) including the effects of residual stress due

to coefficient of thermal expansion mismatch. The solid

line is for no residual stress and is the same for tensile and compressive loading.

57
CO

(D

4->O

"o

CO

ISOMPa

200 400 600 800 1000

t(hr)

1.0

0.8

0.6

0.4

0.2

4^ 1.0

0.8

0.6

0.4

0.2

0.00

(a) Triangular ^^^^

o

(b) Honeycomb ^^^^

e

o e

(c) Kagome

( ,

^^^^

o

e

o

(d) Random

^

^ ft/ 0 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0

Fig. 5 Strain rate vs time after loading 15%TiC/2219 A! to 150 IVlPa with an initially elastic matrix for two different values of the particle Young's modulus.

o

0.1 0.2 0.3 0.4 0.5 0.6

Area Fraction p

Fig. 6 Reciprocal of the relative area modulus for rigid

circular inclusions vs matrix area fraction for ordered

arrays (a, b, and c) and overlapping and non-overlapping

random arrays (d) [from Ref. 4].

Fig. 7 Reciprocal of the relative area modulus vs are

fraction of matrix for random, non-overlapping array of

rigid inclusions. Open circles are from simulations and

the solid line is a mean field theory [from Ref. 11].

58
Green's Function BEM Research at Rutgers CMAS and its

Applications to the Modeling of Mechanical and Piezoelectrical

Behaviors of Advanced Materials

Mitsunori Denda

Rutgers University

Mechanical and Aerospace Engineering Department

RO. Box 909, Piscataway, New Jersey 08855-0909, U.SA.

1. Introduction

At Rutgers Center for Computational Modeling of Aircraft Structures (CMAS) sponsored

by FAA we are developing the boundary element method for elasto-plastic multiple crack

problems using micromechanics concepts and complex variables. The method is used for the

damage tolerance analyses of aging aircraft structures. An amalgamation of the method of

singular integral equations for cracks, the boundary element method for the finite boundary,

and the plastic source method for crack-tip plastic deformation is achieved with the help of

micromechanics tools such as dislocations, point forces, and their dipoles. The formulation is carried out in terms of complex variables to facilitate analytical integration of the boundary

and the crack face integrals and the area integral of the plastic source. The method consists

of three independent modules: complex variable boundary element method (CVBEM), crack source method (CSM), and the plastic source method (PSM). The CVBEM, formulated with

the help of a physical interpretation of Somigliana's identity, uses the continuous distributions of point forces and dislocation dipoles. The CSM models a crack by the continuous

distribution of dislocation dipoles along a line. The PSM is the Green's function representation of the plastic deformation by the continuous planar distribution of dislocations.

In the first half of the paper the BEM research at Rutgers CMAS is introduced. The

power of complex variable techniques coupled with micromechanics is demonstrated in the

derivation of the fundamental solutions and the Green's functions for the BEM. The extension

of the methods to anisotropic elasticity and piezoelectricity is suggested in the second half

of the paper.

59
2. Micromechanics Tools in Complex Variables

Muskhelishvili's [1] complex variable formalism for plane isotropic elasticity uses two analytic

functions or complex potential functions, (l>{z) and ii^{z), of a complex variable z = x -\- iy to

express the displacement, stress, and strain components according to

2fiu{z) = K(f){z) - z(j)'{z) - ^P{z), (2.1)

and

<7xx + cryy 2

a

= <l>'{z)-\-<f>'{z),

+ ia^y = z<f)"{z) + ip'{z), (2.2)

where fi is the shear modulus, k is given by /c = 3 — 4i/ in plane strain and k = [3 — v) /{I -\- u)

in plane stress in terms of Poisson's ratio u. A prime indicates the differentiation with respect

to z and a bar the complex conjugate.

Consider a point force with the magnitude f = fx + ify (per unit thickness) and an edge

dislocation with the Burgers vector h = bx iby, independently located at ^ in the infinite

isotropic medium. Their complex potential functions are given by [2] = -llog{z-0,

^C- V = -^7log(^-0+7-^, (2.3)

where

{k — —K, 7 = //27r(AC + 1) for point force, ro a\

k = 1, 7 = ii.ib/7r{K + 1) for dislocation. ^

' ^

The corresponding dipole solutions (i.e., the force and the dislocation dipoles) are given by

<i>^'\z;0 = -ld{\og{z-a].

.

• V, . i^^'Hz;0 = -kjd{log{z-C)} + ^dl^-^Y (2.5)

where ,

di- ) = -^i- )d^ + )dl (2.6)

is the total differentiation operator.

The continuous distribution of point forces over an arc L (with arc parameter s) is given,

from (2.3) and (2.4), by

IL

i^^'\z) = K jW)\og{z-i)ds^- jj{s)-^ds, (2.7)

60
with T{s) = t/2'n-{K + 1), where t = -\- ity is the traction. The continuous distribution of

dislocation dipoles is given, from (2.5) and (2.4), by

cf>('Hz) = -

I 7(5)^{log(^-0},

iP^'^z) = _|^^,/{log(^-0}+/^7(5)^|7^}, (2.8)

with 7(s) = ifib/irlK + 1), where b = -\- iby is the dislocation.

3. Boundary Element Method in Complex Variables

3.1. Physical Interpretation of Somigliana's Identity

Consider a body R with its boundary OR subject to the traction T = -\- iTy and the

displacement U = Ux + iUy. According to a physical interpretation of Somigliana's identity

[3], the displacement field in this body is obtained by assuming that the region R is embedded

in an infinite medium and OR, which is simply a line marked out in the infinite domain, is covered by a continuous distribution of point forces with density T and by a continuous

distribution of dislocation dipoles with the Burgers vector U.

We discretize the boundary and interpolate the the boundary displacement and traction

following the standard BEM procedure. For the straight boundary element the boundary

integrals are evaluated analytically providing the explicit BEM formulae in complex variables

[4, 5]. However, it is necessary to separate the real and imaginary parts of the displacement

and the traction before establishing the displacement and the traction boundary equations.

3.2. The Green's Function BEM

The availability of Muskhelishvili's formalism is the major motivation of the complex variable

approach to the BEM which also provides access to the analytic function theory. The

derivation of the 2D Green's functions, if not the fundamental solutions, almost exclusively

depends on the complex variable formalism. (The fundamental solutions that satisfy a set of

special boundary conditions are called the Green's functions.) The analytic Green's function

solutions are limited to simple configuration such as the half-plane, bonded two half-planes

(bimaterial plane), an elliptical hole or inclusion in the infinite body and a crack in the

infinite body, etc. Nevertheless, these Green's functions play important roles for problems

involving stress singularities and concentrations. The combination of these Green's functions

with the method such as the crack source method (introduced later) will enlarge the area of

applicability even further. For example, using the bimaterial domain Green's function in the

crack source method one can analyze the problem of near interface cracks very accurately.

Further extension of the crack source method itself is possible to accommodate interfacial

cracks, which enables us to analyze multiple interfacial cracks. The systematic derivation and

cataloging of these Green's functions is a much waited action in the BEM and mechanics

of materials community. Once the required Green's functions are obtained the boundary

element method is formulated with the help of the physical interpretation of Somigliana's

61
identity. Note that the role played by the fundamental solutions is now simply taken over by

the Green's functions. An example is given below for the Green's functions for the infinite

domain with an elliptical hole, which covers the cases of a circular hole and a straight crack.

Consider an elliptical hole with the semi-axes a and 6 ( a > 6 ) in the x and y directions,

respectively. The Green's functions for the point force F and the edge dislocation 6, located

individually at , are given in the form

<f>\^\z;a = <f>^^Hz;0 + = 4^^'\z;0 + i^\'\^;0. (3.1)

where 4>^'\z;^) and ^p^'Hz;0 are the fundamental solution given by (2.3) and (l)\"\z;^) and

ip\^\z;^) are the image terms. Introduce the conformal mapping function

••^

z = M{w) = R{iu -{-—); i?=i(a + 6), m = (3.2)

w 2 a + 6

that maps the points z and ^ in the z-plane to points w and p in the lo-plane according to

z w = —; ana p = ^ 2R ^ 2R

Then the image terms are given [6] by

i;\'\w;p) = 4l\ra;p)^ + ^^^{ro;pn-^^^4'\w-p), (3.3)

ivj.[wy ow

where

(t)'ii{w] p) = L{w, —) + kL{w, i),

i^n{w;p)

P P

p{l + m.p'^) - p{m -\- p'^) 1 pp{m-p^) w-y

p{7n^ p'^) — mp{m -\-p'^) 1 m ' pp{m — p'^) ^ ~

i^ni^'^P) - kL{w,—) + L{w,]:),

P P

and the constants k and 7 are given by (2.4). The function L{w, r) is defined by

.

. . , L{w, t) - log{w — r) — log w,

where r = m/p or 1/p and the function M^{w) by

MM = Mi-) = R (mw + -)

w \ wj

The BEM which uses the Green's functions in place of the fundamental solutions will

be called the Green's function BEM. For the elliptical hole problem it is formulated in the

same way as the ordinary BEM on the bases of the physical interpretation of Somigliana's

identity. The advantage is that the boundary of the hole need not be modeled since the

traction-free boundary condition is automatically satisfied. This is especially useful for the

crack problems (6 = 0) for which the crack-tip stress singularity is build in the Green's

function and the explicit formula for the stress intensity factor is available.

62
4. Crack Source Method

4.1. Single Crack in the Infinite Body

Since the dislocation dipole over an infinitesimal segment gives rise to a displacement

discontinuity, it is considered as a source of crack opening displacement and is called the

crack source. In addition, the continuous distribution of the crack sources over an arc is called the crack element. Consider a straight center crack of length 2a subject to a self- equilibriating traction

±t=±{;.}. (4.1)

in the. global coordinate system xy, where ± indicate the upper and the lower surfaces of the

crack. In the non-dimensional local coordinate system XY (normalized by a), with the origin

at the crack center and the X-axis along the crack, the crack extends in the range — 1 <

X < +1 and the density function of the crack element is given by 7(A^) = ij.i6{X)/7r{K -f 1)

in terms of the crack opening displacement 6 — Sx -\- iSy- In order to build in the known

crack-tip behavior into the solution, we interpolate the density function by

in

where Um-iiX) is Chebyshev polynomial of the second kind. When this is substituted into

the complex potential functions (2.8), the resulting integrals are of Cauchy-type which can be evaluated analytically [7] with the result

+ 1 m-l

where

T("^)(Z) = [Z-VZ^-l)" (m>0),

, , (Z -y/Z^-lV

U(m-i)^Z) = yz2_i {^^>0)^ (4.2)

and Z = X -i- iY is the non-dimensional complex variable. The traction on the upper and

the lower faces of the crack is given by

{tx+itYf{X)^±-^^f:7n6^^^Um-i{X) {\X\<1), (4.3)

which is used to determine the crack opening displacement coefficients ^^""^ by collocation

method. Once the crack opening displacement is determined the stress intensity factor is given by

A'(±l) = Kii±l) + iKn{±\) = -^#E(±ir-'-m « + 1 V a

63
4.2. Multiple Cracks in the Infinite Body

Consider the problem of A'^ multiple straight cracks, ( j = 1, . .

. , TV ), in the infinite

body when the individual crack surface is loaded according to (4.1). First, the traction

contribution from the individual crack is obtained in the associated local coordinate system

as described above. This is transformed to the global coordinate system so that contributions

from all the cracks are assembled giving rise to the total traction t^^^"*" = {i^'^^"'', ^y'^^"^}'^ on

the upper surface of the crack Lk,

- t^''* = z{j:"tx]' (4-4)

where is a coefficient matrix derived in [7]. The unknown crack opening displacement

coefficients = ^^(j)^x^^ij)l}^ i j ~ l?---!-^; = ^i---^P(j)) ^re determined, from

Equations (4.4) and (4.1), by collocation.

4.3. Effect of the Finite Boundary

We now consider multiple center cracks in a finite body R whose boundary OR is subject

to the traction T = {Tj;^Ty}^ and displacement U {U^, Uy}^ while each crack surface is loaded according to (4.1). The total traction on the upper surface of crack Lk is given in the

form

N ( P{j) ) M

tW+ = E E«S'C +E{Gf'T„-^Hf'U„}, (4.5)

j=l (^m=l J n=l

where the first term in the right hand side comes from (4.4) and the second term comes from

the traction BEM. The quantities U„ and T^, are the boundary displacement and the traction

vectors. The coefficient matrices G*J^^\ H*^''*^ are derived in [7]. The total displacement on

the non-crack boundary dR is given in the form

N (

P(j) 1 ^ f 1 1 2U. = E E ^w'*!? + E 7G-.T„ - H„U„ , (4.6)

j=l [m=l J n=l I }

where Gn, and H„ are coefficient matrices derived in [7]. The solution is obtained by

setting up traction equations on the upper surface of each crack from equations (4.5) and

(4.1) and displacement boundary equations on dR from equation (4.6) and the boundary

condition on dR.

4.4. Numerical Results

A single center crack in a plate in uniaxial tension was analyzed b)^ the present method using one Chebyshev polynomial. The same problem, with identical mesh, was analyzed by the

crack Green's function BEM, which uses the Green's function that satisfy the traction free

crack surface boundary condition automatically. The stress intensity factor results agreed up

to seven significant digits. Figures 1, 2, and 3 show two collinear cracks, two parallel cracks,

64
and three parallel cracks, respectively, in the infinite body. The numerical results have been

obtained for a large plate, compared to the cracks, using seven Chebyshev polynomials for

each crack in each case. Comparison of the numerical results and the results from the stress intensity handbook [8] is listed in Tables 1, 2, and 3. The crack source method is as accurate

as the crack Green's function BEM for single crack problems. Furthermore, it is the natural

extension of the crack Green's function BEM for multiple crack problems.

5. Plastic Source Method

In order to deal plane plasticity problems correctly it is necessary to introduce the fictitious

in-plane plastic strain components (or the plane plastic strain) [2]

^a(3 "t" ^^33^a0 (plane strain)

^a0 = < ( c^, 1,2 ), (plane stress)

instead of the 3D plastic strain components ef^. Modified elastic constants

A* = i

A (plane strain)

. (plane stress)

u (plane strain)

1^ (plane stress)

are used instead of Lame constants, A and fi, and Poisson's ratio v [2]. The plastic deformation in an infinitesimal area dA at ^ is imagined to be a source of

residual stress and is called the plastic source. The fundamental solution for the plastic source is given by [2]

d(l>*{z) T.^l0g(2 - <^)C?A,

1 dA, (5.1)

where

1

(7 =

2tt{k + 1)

^77 =

1 2-k{k + 1)

7r(«; + l)

7r(« + 1)

(^11 ^22 + 22£i2) • (5.2)

Consider a region D of plastic deformation where the plane plastic strain components e*^

are prescribed. The complex potential functions reflecting the effect of the plastic deformation in D are given by integrating (5.2) over D to give

r(^) = /|(.-lK|log(.-0 + .-|(j^ dA, (5.3)

65
when z is outside the plastic region. Additional treatment required when z is in the plastic

region was discussed by Denda [2]

Application of the divergence theorem to (5.3) will result in an alternative representation

of the plastic deformation in D which consists of the effective body force and traction distributions in D and on its boundary dD^ respectively. The potential functions for the effective

body force distribution is given by

and those for the effective traction on the boundary by

... <i>\[z) ^ [ {a*d^-T.d()log{z-0,

ll JdD ^

^

/ • rA^) = f - ^*dA log{z -0-^^ f U*d^ - r.d^) (5.5)

Zl JdD ^ ^ ll JdD ^ ^ Z — C, where z is assumed to be outside the plastic region.

In the numerical implementation the plastic region is discretized into a collection of

triangular or quadrilateral plastic elements in each of which the plastic strain distribution is approximated by an interpolation function. The analytical integration of (5.3) was performed

in [2] for the constant interpolation function and the case of linear interpolation coupled with

the triangular plastic element is treated in [9] using the effective force method. The linear

interpolation of the plastic strain used in the mesh of triangular plastic elements guarantees

the continuity of the plastic strain across the boundary of two adjacent plastic elements.

This results in the cancellation of the effective traction on the boundary with the exception

where the plastic strain is inherently discontinuous such as on the crack surface. For crack

problems the effective body force in the plastic elements and the effective traction on the

crack surface represents the influence of the crack-tip plastic deformation completely.

The solution procedure whereby the unknown plastic strain distribution is determined was

given by Denda and Lua [10] for standard elasto-plastic problems. In order to use the crack

-source method it is convenient to break the problem into the elastic and the plastic solutions.

The former is the elastic solution under the applied load and the latter the solution of the

plastic elements. The crack source method is used for the determination of each solution, once for the elastic solution and several times, iteratively, for the plastic solution. Note

that each solution gives rise to a 1 / x/r stress singularity. A singularity cancellation scheme,

whereby the final solution is obtained by elimination of the total stress intensity factor, is used as a part of the convergence criterion of the procedure as reported in [11].

6. The BEM in Piezoelectricity and Anisotropic Elasticity

The micromechanics and complex variable approach to the Green's function BEM for plane

isotropic elasticity (and plasticity) can be generalized to plane anisotropic elasticity using

66
the Stroh formalism [12]. Further generalization to piezoelectricity is described below. The

key results for the anisotropic elasticity is obtained from those of the piezoelectricity by

setting the piezoelectric constant to zero, which has the effect of decoupling the elasticity

and electricity solutions. Basic equations in piezoelectricity and the notation used below are summarized in Appendix.

6.1. Fundamental Solutions in Piezoelectricity

The fundamental displacement solution for a point force fi [i = 1,2,3) and point charge p

located at = Ci + K2 is given by

u(^)(^;0^| |^'^ = 23?{^Alog(Z-S)A^f}, (6.1)

where f = {fi,<i>}'^ is a 4 x 1 column vector and log(Z — S) = Diag {log(ZQ — H^)} is a 4x4 diagonal matrix with = Ci +PoC2- The solution for a dislocation bi {i = 1,2,3) and

a potential jump A</) at ^ is given by

u(^)(^;e) =

| }^^ = 23?{^AIog(Z-S)B^b}, (6.2)

where b = {6,-, A^}'^ represents displacement and potential jumps for an arbitrary counterclockwise circuit surrounding ^. These solutions are obtained following the complex variable

formalism of piezoelectricity outlined in Appendix.

The dipole solutions are obtained by applying the total differential operator (2.6). to

the fundamental solutions as in the case of isotropic elasticity. A force dipole and a charge

dipole located at ^ has the solution

uf{z; 0 = 23?

{

^A dlog(Z - S) A^f} , (6.3)

while the solution for a dislocation dipole and a potential jump dipole at ^ is given by

ui'\z-, 0 = 2^[ ^A dlog(Z - S) B^b} , (6.4)

where

dlog(Z — S) = Diag {c?log(Za — Ha)} = Diag — g^ni on a).

The continuous distribution of the point force and charge over a contour L has the solution

u^;Hz) = |^23?|^Alog(Z-S) A^f}tf5, (6.5)

where 5 is the arc length. The solution of the continuous distribution of the dislocation

dipoles and the potential jump dipoles over L is given by

ul'\z) = |^25R|^Adlog(Z-S)B^b}. (6.6)

67
6.2. Physical Interpretation of Piezoelectric Somigliana's Identity

In place of the displacement {ui}^ and the traction {ti}^ in elasticity, vectors {ui,(l)}^ and

{ti, —Lo}^- are used in piezoelectricity, which may be called the piezoelectric displacement

and traction, respectively. The generalization of Somigliana's identity in elasticity to piezoelectricity along with its physical interpretation can be performed in the same fashion. Thus

the piezoelectric displacement in the domain R is given in terms of the layers of piezoelectric

traction and displacement over the boundary OR which are embedded in the infinite body.

The formulation of the piezoelectric BEM can be performed in the same way as in elasticity.

7. Extension to 3D

In 3D BEM the micromechanics approach introduced earlier coupled with the physical interpretation of Somigliana's identity still is promising. However, the lack of 3D formalisms

comparable to those in 2D (Muskhelishvili [1], Stroh [12], Barnett-Lothe [13]) is the major

handicap in the derivation of the Green's functions. We hope that a 3D formalism as effective

and manageable as those in 2D will be established in the future. Those who are working on

such a venture are Piltner [14], Giirlebeck and Sprossig [15], and Zhdanov [16] among others.

Acknowledgment

This work was supported in part by the FAA Center for Computational Modeling of Aircraft Structures and by the Department of Aerospace and Mechanical Engineering, Rutgers

University.

References

[1] N.I. Muskhelishvili. Some Basic Problems of the Mathematical Theory of Elasticity.

Noordhoff, Groningen, 1958.

[2] M. Denda. Complex variable Green's function representation of plane inelastic deformation in isotropic solids. Acta Mechanica., 72:205-221, 1988.

[3] M. Denda. A complex variable approach to inelastic boundary value problems. In

H. Chung et al., editor. Advances in Design and Analysis in Pressure Vessel Technology,

ASME PVP - Vol. 130, NE - Vol. 2, pages 23-32, New York, 1987.

[4] M. Denda and Y. F. Dong. A boundary element method for multiple crack problems

based on micromechanics and complex variables. In Proceedings for BETECH 94 (9th

International Conference on Boundary Element Technology)., pages 227-234, 1994.

[5] M. Denda and Y. F. Dong. Explicit BEM formulae in complex variables, (submitted

for publication)., 1994.

68
[6] Denda M. Two-dimensional singularity solutions in the presence of an elliptical hole.

J. AppL Mech., 56:231-233, 1989.

[7] M. Denda and Y. F. Dong. Micromechanics and complex variable approach to the

boundary element method for multiple crack problems, (submitted for publication),

1994.

[8] Y. Murakami et al. Stress Intensity Factor Handbook. Pergamon Press, Oxford, 1987.

[9] M. Denda and Y. F. Dong. Explicit complex variable formula for nonuniform plane

plastic deformation, (in preparation), 1994.

[10] M. Denda and Y. J. Lua. Formulation of the plastic source method for plane inelastic

problems, part 2: Numerical implementation for elastoplastic problems. Acta Mechanica, 75:111-132, 1988.

[11] M. Denda and Y. F. Dong. Principle of superposition for elastic-plastic crack analysis.

(in preparation), 1994.

[12] A. N. Stroh. Dislocations and cracks in anisotropic elasticity. Phil. Mag., 3:625-646,

1958.

[13] D. M. Barnett and J. Lothe. Synthesis of the sectic and the integral formalism for

dislocations. Green's function and surface waves in anisotropic solids. Phys. Norv.,

7:13-19, 1973.

[14] R. Piltner. On the representation of three-dimensional elasticity solutions with the aid

of complex valued functions. J. of Elasticity, 22:45-55, 1989.

[15] K. Giirlebeck and W. Sprossig. Quaternionic Analysis and Elliptic Boundary Value

Problems. Birkhauser Verlag, Basel, 1990.

[16] M. S. Zhdanov. Integral Transforms in Geophysics. Springer-Verlag, Berlin, 1988.

[17] Z. Suo, C. M. Kuo, D. M. Barnett, and J. R. Willis. Fracture mechanics for piezoelectric

ceramics. J. Mech. Phys. Solids, 40:739-765, 1992.

69
Table 1. Two Collinear Cracks {Fja = KiA/cry/Tra and Fjb — Kib I

2a/d FlA FlA Fib Fib

(Handbook) (Numerical) (Handbook) (Numerical)

0.05 1.00031 1.0018 1.00032 1.0018

0.1 1.0012 1.0027 1.0013 1.0028

0.2 1.0046 1.0061 1.0057 1.0071

0.3 1.0102 1.0117 1.0138 1.0153

0.4 1.0179 1.0194 1.0272 1.0287

0.5 1.0280 1.0295 1.0480 1.0495

0.6 1.0410 1.0426 1.0804 1.0821

0.7 1.0579 1.0596 1.1333 1.1351

0.8 1.0811 1.0827 1.2289 1.2314

0.9 1.1174 1.1187 1.4539 1.4639

Figure 1: Two collinear cracks in the infinite body under uniaxial tension.

70
Table 2. Two Parallel Cracks [Fj = Ki

(T\/t^)

2a/d Fi Fi

(Handbook) (Numerical)

0.0 1.0000 1.0011

0.2 0.9855 0.9870

0.4 0.9508 0.9517

0.8 0.8727 0.8732

1.0 0.8319 0.8440

2.0 0.7569 0.7746

5.0 0.6962 0.7129

d

< 2a

u

2: Two parallel cracks in the infinite body under uniaxial tension.

71
Table 3. Three Parallel Cracks {Fia = Kia/(^\/t^ and Fjb = KiB/cry/Tta)

(Handbook digital values for Fjb are not available.)

2a/d Fia

(Handbook)

Fia

(Numerical)

Fib

(Handbook)

Fib

(Numerical)

0.1 0.99500 0.99687 — 0.99410

0.2 0.98198 0.98379 0.97306

0.3 0.96299 0.96430 0.94156

0.4 0.94010 0.94100 0.90361

0.5 0.91535 0.91650 0.86789

0.6 0.89080 0.89254 0.82333

0.7 0.86851 0.87041 0.78603

0.8 0.85052 0.85062 0.75234

t'

d

d

< 2a

I'

Figure 3: Three parallel cracks in the infinite body under uniaxial tension.

72
Appendix

A. Piezoelectricity

A.l. Basic Equations

For a linear piezoelectric material the stress cr and the electric induction D are given by the

strain e and the electric field E by

<7ij = Cijki€ki — ekijEk, (A.l)

Di = CikiEki + UkEk^ (A. 2)

where Cijki-, Giki, and Cik are the elastic stiffness, piezoelectric, and dielectric permittivity

constants, respectively. Summation convention over the repeated indices is implied. The

strain is given by the displacement u and the electric field E by the electric potential (f) according to the relations

where a comma followed by an index i indicates differentiation by the coordinate Xi. The

stress and the electric induction satisfy the equilibrium equations

+ fj = 0, A,i - p = 0, (A.4)

where f and p are the body force and extrinsic bulk charge densities, respectively. Across an

interface of the piezoelectric domain + and the adjacent domain — , the jumps in the stress and the electric induction are given by

t% n, \dI - D-] = -co', (A.5)

where n is the unit normal to the interface pointing from the + side; t° and co' are the force

and the charge per unit area of the interface. Obviously, the interface jump conditions can be rewritten as the boundary conditions of the domain in the form

maji = tj, UiDi = —u. (A. 6)

A.2. Complex Variable Formalism

The equilibrium equations (A.4) in 2D, with no body force and bulk charge, are given by

Da,a = {eaipui - tap(f>) ^ Q ( f , j = 1, 2, 3; cv, ^ = 1, 2), (A. 7)

using (A.l), (A. 2), and (A. 3). The general solution of (A. 7) is obtained in the form

{W„ (t>f = Rf{Z), Z = mi.Ti + 1712X2, (A.8)

73
where a = {o.i}^ and {iii, (j)}'^ are column matrices of dimension four and f{Z) is an analytic

function of Z. Without loss of generality we can set

•

' " mi — I, 7112 — PThe stress and the electric induction are given by

'', . • (Taj — {Cajipai -\- epjaa4) mp f'{Z),

Da = {eatpai - ea0a4)mp f{Z).

The quantities p and a are determined by substituting (A. 8) into (A. 7), which yields

(A.9)

or in a matrix form

(Cajipai + epjaa4)mpma - 0,

{eaipai - eapa4)mpma = 0,

{q+p(R + R^) +/T}a-0,

(A.IO)

(A.ll)

where

Qik = Ciikl, Rik = Ciik2i Tik = C2ik2i

Qi4 = Q4i = Ciii, Ri4 = e2il, i?4i = 61^2, Ti4 = T^i = 6212 (hk = 1,2,3),

Q44 -- —cll, R44 — —612, 744 = —€22-

A nontrivial a exists if p is a root of the determinant polynomial

. |q + p(r + r^) +/t|| -0.

The eigenvalue a is obtained from (A.ll). Write (A.ll) as (Q+pR)a= -p(R^+pT)a

and introduce a new vector

b = (r^ + pT) a = -i (Q + pR) a.

V / p

Then, (A.9) can be written as

{au,Dif =-phf'{Z), {a2^^D2f = hf'{Z).

Introduce the stress functions

(A.12)

(A.13)

(A.14)

Wi,<i^4y

(A.15)

(A.16)

(A.17)

such that

d

0x2

_d_

dxi

74
Then, equilibrium equation (A. 7) are automatically satisfied by the stress functions (A. 17)

and, from (A. 16), we find = hf{Z). (A.18)

It is shown that [17] the polynomial (A. 14) cannot have real roots and that the roots

consist of four pairs of complex conjugates denoted by

Pa+4=Pa = 1,2,3,4).

The corresponding eigenvectors have the similar property

aa+4 (o; = 1,2,3,4),

each of which introduces a new vector

ba+4 = ba {q - 1,2,3,4).

If we define the 4x4 matrices A and B by

A = [ai,a2,a3,a4] , B = [bi, b2, as,

then the general solution of (A. 7) is given in the form

{ui,4>f = Af(Z) + A f(Z) = 2Sft{Af(Z)} (A.19)

and the stress functions in the form

Wi,<P4f = Bf(Z) + B f(Z) = 2SR{Bf(Z)} (A.20)

where

i{Z) = {faiZ,)f, Z^ = x,+paX2 (a -1,2,3,4),

and f (Z) = f(Z).

75

APPROXIMATE OPERATORS FOR

BOUNDARY INTEGRAL EQUATIONS

IN TRANSIENT ELASTODYNAMICS

Thomas L. Geers and Brett A. Lewis

Center for Acoustics, Mechanics and Materials

Department of Mechanical Engineering

University of Colorado, Boulder, CO

Abstract

The integral equation of three-dimensional transient dynamics for an isotropic elastic

medium has been known for about a century, and numerical methods for solving the equation have

been under development for about a quarter-century. However, applications have been limited to

idealized problems because of the complexity of the integral operators and the intense consumption

of computational resources. The state of development for anisotropic media is much further

behind. The present situation has motivated searches for accurate approximate operators that

facilitate straightforward implementation and rapid computation. One such search is described

herein, which involves the formulation and evaluation of singly and doubly asymptotic operators

for unbounded domains. The former approach exactness at either early time or late time; the latter approach exactness at both early and late time. Singly asymptotic operators yield satisfactory

results only in restricted circumstances, whereas doubly asymptotic operators have proven to be

quite robust.

Introduction

Materials characterization, flaw detection, medical diagnosis, earthquake-resistant construction, oil exploration and defense technology are some of the areas in which dynamic boundaryelement analysis has been productively used. Considerable progress has been made in timeharmonic applications, because the pertinent integral operators facilitate accurate and efficient

numerical solution. However, computational obstacles have hindered the development of transient

boundary-element analysis, which has motivated the development of approximate operator

representations.

77
This paper addresses the development of approximate integral operators in threedimensional transient elastodynamics, following similar work in transient acoustics (see, e.g., Geers and Zhang, 1994). The focus here is on a body embedded in an unbounded, isotropic,

linear-elastic medium. The body itself, and/or some medium in the immediate vicinity of the body,

may exhibit nonhnear behavior; in the latter case, the boundary on which the integral operators act encloses both the body and the nonlinearly behaving medium. In the discussion following, the

domain inside the integral-operator boundary will be called the body, whether or not it includes any

medium.

Outside the integral-operator boundary, the total elastodynamic field is separated into the

(known) incident and (unknown) scattered fields. The incident field is that which would exist if

the body were replaced by linear-elastic medium. The scattered field, then, is merely the

difference, at any point and any time, between the total field and the incident field, i.e., it is the

field caused by the presence of the body. Treatment of the scattered field on the integral-operator

boundary is then greatly facilitated by the formulation of a temporal impedance relation (TIR) that

provides the body's time-domain view of the surrounding medium. When the TIR is combined

with the equations of motion for the body, the interface compatibility conditions at the integraloperator boundary, and the pertinent initial conditions, a complete mathematical formulation is

obtained that lends itself to efficient numerical solution.

" Exact TIR's are compUcated and costly to employ in computations, as they are nonlocal in

both space and time, i.e., they require full computational matrices and long-memory response data.

Hence, accurate approximate TIR's are needed for efficient computation. This paper is about

approximate TIR's that approach exactness at early time, at late time, or both. Hence, the

developments below pertain exclusively to scattered fields.

Three-Dimensional Transient Elastodynamics

In the absence of body forces, the scattered displacement field in an infinite, isotropic,

elastic medium may be expressed in terms of a scalar potential (()(x,t) and a vector potential

Y(x,t) as the Helmholtz decomposition (see, e.g, Eringen and Suhubi, 1975)

ii(x,t) = V(l)(x,t) + Wx\if(x,t), V-\iif(x,t) = 0. (1)

With this decomposition, the displacement equation of elastodynamics separates into the uncoupled

wave equations

clV^(p(x,t)='^(x,t), clV^\i/(x,t)=^(x,t), (2)

78
where an overdot denotes a time derivative, and Cq and Cg are the dilatational and shear wave

speeds, respectively, given in terms of the Lame parameters X, \i and the mass density p by

ci=(X + 2^L)lp 4 = ^ I p. (3)

The associated stress tensor may be expressed

G(xj) = XV^<l>(x,t)I +2ix^^(p(xj)^- ix(V x\lr(x,t)] + (4)

where / is the identity tensor. An exact TIR for both (pixj) and \i/(x,t) is Kirchhoffs retarded potential formula

(Eringen and Suhubi, 1975)

27i(p(F,t)= \{R-^^(P,tj^) + R-2-^[(p(p,tj^) + ^(p(P,tji)]}dSp, (5)

J orip dnp c

S

where (p is either 0ot ifr and c is either cd or cs, F is a field point and P is an integration point

(both on the smooth boundary S), R is the distance from P to F, np is the surface normal at P, and

the retarded time tj^ =t-R/ c. An exact TIR that directly links the boundary displacement vector

u(F,t) and the boundary traction vector t(F,t) is Love's integral identity, which may be written

in Laplace-transform space as (Cruse and Rizzo, 1968)

ja(F,s) + j

U(P,s)f(F,P,s) dSp = Jt(P,s)U(F,P,s) dSp , (6)

5 s

where T(F,P,s) and U(F,P,s) arc tensor Green's functions. As mentioned above, computational

implementation of these TIR's is currently impractical.

Early-Time Approximations ETAi and ETA2

An analysis of (5) at early time (Felippa, 1980a; Nicolas-Vullierme, 1991) for a continuous surface field on a smooth boundary yields the curved-wave approximation

(p(F,t)+ K(F)c(p(F,t) « -c^(F,t), (7) on

where k is the mean curvature of the boundary. Now (7) is also produced by ray theory (NicolasVullierme, 1991), which yields in addition (Lewis, 1994)

79
^ ^ dn^ dn RlRr ^ dn^ dn r} R,

Nonzero Component of VV0: fVV0j„„ « —j, (8)

dn

Nonzero Components of VfVx y/j+fV x y/jV:

- - d^Vr 1 ^Wf ) [V(VxY)+(^xii/)W]„^ = mvxv/;+rvxvA;V7^„- '"dT"^^ d^'^^^^^'

/VfVxv>;+fVxi^;V7„^ = /vrvxv>;+fVx^;V7^„== ^^+R-^^-Rfyf^,

where R^ and R^ are the principal radii of curvature at the field point, ^ and ^ are the

corresponding surface coordinates, and n is the surface normal

With boundary tractions /„, and on the body being the negatives of (T„„, (T„^ and <J„^,

respectively, the preceding equations may be manipulated (Lewis, 1994) to yield ETA2 in local

coordinates

V(F,t)+K(F)CV(F,t) = pCU'(F,t) + 2fXK(F)b'ti'(F,t), (9)

where t' = {t„ t^} and u' = {u„ u^} , and C and D' are diagonal matrices with nonzero

elements Cd,Cs,Cs and 2,1,1 respectively. ETA2 may be expressed in global Cartesian coordinates

as (9) without primes, in which t(F,t) = Q^ (F)t'(F,t), u(F,t) = ff(F)u'(F,t),

C(F) = Q^ (F)CQ(F) and D(F) = Q^(F) D'Q(F), where Q(F) is the appropriate coordinate

rotation matrix.

A lower-order ETA may be obtained from ETA2 by noting that, at very early time, the first terms on either side of (9) overshadow their lower-derivative counterparts. Thus, neglecting the

second terms on either side of (9) and integrating in time, one obtains ETAi, expressed in global

coordinates as (Underwood and Geers, 1981)

t(F,t)=pC(F)U(F,t). (10)

80
Late-Time Approximations LTAi and LTA2

This approximation is obtained by expanding the Green's functions of (6) in Taylor series

to get

f(F,P,s) = f^(F,P) + 0(s^),

(11)

U(F,P,s) = U^(F,P) + sUUf,P) + 0(s^).

Note that the 0(s^ ) term in the f - expansion is zero. Introducing these expansions into (6) and

inverse transforming, one obtains LTA2 as

A(F,P)f(P,t) + B(F,P)t(P,t) = r(F,P)u(P,t)^

(12)

where the spatial operators A(F,P), B(F,P) and r(F,P) are defined by

A(F,P)f(P,t) = jf(P,t) U^(F,P)dSp, B(F,P)t(P,t) = jt(P,t)U^(F,P)dSp,

s s r(F,P)u(P,t) = j

u(P,t) [S(F -P)+f^(F,P)] dSp. (13)

s A lower-order LTA may be obtained from LTA2 by noting that, for slow motions at very late

time, the second term on the left overshadows its higher-derivative companion. Thus, neglecting

the first term on the left side of (12), one obtains LTAi

B(F,P)t(P,t) = r(F,P)u(P,t). (14)

Doubly Asymptotic Approximations DAAi, DAA2 and DAA1.2

These approximations are obtained by a procedure called operator matching, which is akin

to solution matching in the method of matched asymptotic approximations (Van Dyke, 1964). A

scalar form of the method was presented in Felippa, 1980b, extended to operator form in NicolasVullierme, 1991 and refined in Geers and Zhang, 1994. The procedure is outlined here for DAAi;

the more complicated extension to DAA2 and DAA1.2 is described in Lewis, 1994.

To begin, (10) and (14) are Laplace-transformed to yield

t(F,s) = pCiF)sn(F,s),

(15)

t(F,s) = B-^(F,G)r(G,P)U(P,s).

Next, a DAA trial equation is chosen as

81
t(F,s) = [sUj(F,P) + Uo(F,P)]u(P,s), (16)

where the operators Uo(F,P) and U](F,P) Bcce unknown; (16) is then expressed in the asymptotic

forms

, . f(F,s) = [Uj(F,P) + 0(s~^)]sn(P,s), s^oo,

(17)

f(F,s) = [Uo(F,P) + 0(s)]U(P,s), 5=>0.

The first of these is now matched to the first of (15) to obtain Uj(F,P) = pC(F)5(F - Pj and the

second is matched to the second of (15) to obtain Uq(F,P) = B~^(F,G) r(G,P). Finally, these

results for Uq(F,P) and Ui(F,P) 3(rt introduced into (16) and the resulting equation is inversetransformed to obtam DAAi as

t(F,t) = pC(F)U(F,t) + B~^(F,G)r(G,P)u(P,t). (18)

Through a similar procedure, DAA2 is obtained by assuming the trial equation

sf(F,s) + fo(F,P)f(P,s) = [s^U2(F,P) + sUj(F,P) + Uo(F,P)]a(P,s), (19)

matching the two leading terms for 5 => <» to their counterparts in the Laplace-transformed, global

form of (9), and matching the two leading terms for s=>0 to their counterparts in the Laplacetransformed (12). The result is

t + Qt = pen + (B~^r- QB~^AB~^r)U + QB'^rn, (20)

in which Q = [B~^f + K(pC^ - D)J(pC + B~^A B'^tf^; here, in the interest of compactness, spatial and temporal dependencies have been suppressed. A useful simplification of DAA2 is

effected by taking A=0, which corresponds to matching the two leading terms for j' => to their

counterparts in the Laplace-transformed, global form of (9), but only the leading term for s =>0 to

its counterpart in the Laplace-transformed (12). This doubly asymptotic approximation is DAA1.2.

Boundary Element Analysis

All of the seven approximate temporal impedance approximations given above lend

themselves to numerical computation through boundary element discretization. The discretized

forms appear as (10), (9), (14), (12), (18), and (20) with A=0 or A^O, except that the spatial

vectors and tensors are replaced my computational vectors and matrices, respectively. Thus, each

discrete TIR reduces to an ensemble of ordinary differential equations in time.

82
Example Problem

Perhaps the simplest transient elastodynamics problem for an unbounded domain is the

suddenly pressurized spherical cavity (Timoshenko and Goodier, 1970). Analytical DAA1.2 and

DAA2 treatments of this problem turn out to be equivalent, and produce, in fact, the exact solution.

The figure below shows the coincident DAA1.2 and exact analytical displacement histories, as well

as corresponding DAAi and DAA1.2 boundary-element results. The boundary-element DAAi

history, which has been verified as virtually identical to its analytical counterpart, is clearly

inadequate, whereas the boundary-element DAA1.2 history lies very close to the coincident

analytical histories. These results are representative of behavior observed in studies of more

complicated idealized problems (Lewis, 1994; Oberai, 1994).

Conclusion

Four singly asymptotic and three doubly asymptotic approximations have been presented as computationally practical substitutes for the exact temporal impedance relation pertaining to a body

embedded in an infinite, isotropic, elastic medium. Among these approximations, only DAAi.

and DAA2 provide levels of accuracy sufficient for treating problems involving broad-band

excitations.

2.0 1

' ^ :

'

\ ^

' 1

f5

0 5 10 15

Time Cpt/a

Radial Displacement of a Step-Pressurized Spherical Cavity in an

Infinite, Isotropic, Elastic Medium: Exact, DAAi and DAA1.2 Solutions

83
Acknowledgments

The study described herein was sponsored under Contract DNA001-92-C-0038 with the Defense

Nuclear Agency, Alexandria, VA. Computing resources were provided by Aptek, Inc., Colorado

Springs, CO.

References

Cruse, T.A., and Rizzo, F.J., 1968, "A Direct Formulationa dn Numerical Solution of the General

Transient Elastodynamics Problem", J. Math. Anal. App., Vol. 22, pp. 244-259.

Eringen, A.C., and Suhubi, E.S., 1975, Elastodynamics - Linear Theory Vol.2, Academic Press,

New York.

Felippa, C.A., 1980a, "A Family of Early Time Approximations for Fluid Structure Interaction",

/. Appl. Mech., Vol. 47, pp. 703-708.

Felippa, C.A., 1980b, "Top-Down Derivation of Doubly Asymptotic Approximations for

Structure-Fluid Interaction Analysis", Innovative Numerical Analysis for the Engineering

Sciences, Shaw, R.P., etai, eds.. University Press of Virginia, Charlottesville, pp. 79-88.

Geers, T.L., and Zhang, P., 1994, "Doubly Asymptotic Approximations for Submerged

Structures with Internal Fluid Volumes: Formulation" (with P. Zhang), to appear in J. Appl.

Mech., Vol. 61, No. 4.

Lewis, B.A., 1994, Doubly Asymptotic Approximations for Non-Reflecting Boundaries in

Transient Elastodynamics, Ph.D. Dissertation, University of Colorado, Boulder.

Nicolas-Vullierme, 1991, "A Contribution to Doubly Asymptotic Approximations: An Operator

Top-Down Derivation", J. Vib. Acoust., Vol. 113, pp. 409-415.

Oberai, A.A., Transient Response Analysis of Cylindrical and Spherical Shells Embedded in an

Elastic Medium: Residual Potential and Doubly Asymptotic Solutions, M.S. Thesis, University of

Colorado, Boulder.

Timoshenko, S.P., and Goodier, J.N., 1970, Theory of Elasticity, McGraw-Hill, New York.

Underwood, P.G., and Geers, T.L., 1981, "Doubly Asymptotic Boundary-Element Analysis of

Dynamic Soil-Structure Interaction", Int. J. Solids Struc, Vol. 17, pp. 687-697.

Van Dyke, M.D., 1964, Perturbation Methods in Fluid Mechanics, Academic Press, New York.

84
Simulation of the Electrochemical Machining Process

Using a 2D Fundamental Singular Solution

H.A. Nied

Aug. 1994

GE CORPORATE RESEARCH & DEVELOPMENT

Schenctady, NY 12301

Abstract

Since the Electrochemical Machining (ECM) process is essentially a surface phenomenon,

a solution is readily obtained by employing classical potential theory. By using a 2D Green's

function or fundamental singular solution, the boundary value problem for ECM process can be

reduced to the solution of a boundary integral equation. When the surfaces are discretized, a numerical solution for the potential and the derivatives at the surface being machined can be

derived. The boundary element method (BEM) developed by Brebbia was used to solve the field

equations for this non-conventional machining process. The electrochemical anodic reaction was

furnished by Faraday's Law, which provided the relationship for the rate of dissolution of the

workpiece. By this approach, a 2D computer process model was developed for the purpose of

simulating the ECM process. As an example, application of this process model for determining

the machined shape of an airfoil compressor blade is described.

85
Introduction

Electrochemical machining is a process for the removal of an electrically conductive metal

by anodic dissolution using large surface current densities when the anode and cathode are separated by a narrow gap containing a rapidly flowing electrolyte. ECM offers many advantages over traditional machining processes like milling and turning. ECM is a non-contact machining

process that can quickly shape any electrically conductive material regardless of the material's

hardness or toughness. Since the workpiece is shaped through anodic dissolution, the process

does not produce residual stresses in the parent material, resulting in a stress free component.

Additionally, due to the nature of the process, ECM produces a smooth, burr free surface that

usually does not require post machining operations.

The major characteristics of the ECM process are illustrated in Figure (1). In this figure,

which shows two flat electrodes for simplification, metal is dissolved from the anodic workpiece

due to the electric current flow between the cathode and anode by ion transport. The rate of

material dissolution depends entirely on the total charge that is passed to the workpiece. Different

alloys have different rates of material removal based on the dissolution valency of their

constituents. It is worth mentioning that the theoretical material removal rate usually does not

correspond to the experimental rate. This difference is called the current efficiency and is usually

expressed as a percentage of the theoretical amount predicted by Faraday's law.

The primary functions of the electrolyte are to provide a conductive media to carry the

current, remove heat generated from the tool and workpiece, and flush away the dissolution

products and gases generated in the thin boundary layer downstream. As the primary by-product

of the ECM process, hydrogen gas is generated which tends to collect at the cathode. Presence of

hydrogen in the interelectrode gap has the effect of reducing the electrolytes' bulk conductivity.

While the current density controls the amount of hydrogen that is generated, maintaining high back

pressure in the gap can minimize the effect by reducing the size of hydrogen bubbles that are generated. The basics of electrochemistry applied to the ECM process, including equilibrium

electrode potentials, irreversible electrode reactions, polarization and overpotential are described in

detail by both De Barr[l] and McGeough[2].

The ECM process has been used to fabricate difficult to machine mechanical components.

In particular, it is very difficult to machine airfoils made of super and titanium alloys using

conventional machining methods. Figure (2) shows schematically the typical tooling arrangement

used for machining airfoils and blades in the aircraft engine industry using cross flowing

electrolyte. Electrolyte under high pressure flows in the gap between the airfoil and the electrodes

from the leading to trailing edge of the airfoil. ECM process parameters such as temperature,

voltage, feedrate, and pressure are tightly controlled during the ECM process. Other parameters

such as electrolyte conductivity, pH, and turbidity are monitored and adjusted as necessary.

Cathodes for the convex and concave shapes are usually feed from opposite sides at equal

feedrates during the process to maintain a prescribed rate of material removal.

Since the design of ECM tools is currently a costly iteration process, a 2D ECM process

model was developed for the purposes of improving general process understanding, predicting

machined shapes, aiding in the design of cathode tools, and process optimization.

Modeling Description

Process ModeUng Background:

Early ECM simulation was focused on the solution of the Laplace equation for single

electrode geometries and mathematical methods for cathode design. It was recognized by these

researchers, Ref. [3]-[71, that the accuracy for the prediction of anode material dissolution was

highly dependent on the current density distribution. The approach used by these authors was

86
analytic mathematics. For example, Nilson and Tsuei[41 used complex variables for obtaining

solutions to simple single tool geometries. None of these highly mathematical methods could be

practically used for simulating the entire complicated ECM process. A numerical method was

introduced by Jain and Pandey Ref.[81-[12] using the finite element technique which would allow

handling more complicated geometries. This approach, although superior to previous efforts,

required discretizing the interelectrode region which in most ECM processes is extremely small, in

the order of 0.010 inch. Narayanan et al Ref.[131-[14] modeled the ECM single tool process by

using the boundary element method. The BEM provides a more efficient numerical method for

updating only the workpiece surface rather than remeshing the entire model, as would be required

in finite elements.

All of the afore mentioned models of the ECM process were primarily directed to the

solution of the Laplace equation for the determination of the current density distribution without

detailed consideration of the effects of electrolyte flow on the machining process.

BEM Process Model Development:

Many problems of engineering can be classified as mixed boundary value problems when

the function and its derivative are specified on sections of the boundary surrounding the interior

domain. Such is the case for the ECM process when voltage or currents are specified on surfaces

for single or multiply connected domains. The numerical method usually used for solving such

mixed boundary value problem is by applying the Finite Element Method(FEM) Ref.[91. Using

the FEM method, one discretizes the material interior with a mesh that includes the boundaries. If gradients of the potential near the surface are sought, the interior mesh must be accordingly

refined. The FEM formulation can then be reduced to a set of linear equations for the solution of

the potential at each interior node when boundary condition nodal values are specified on the

boundary. This method provides a solution at each interior point whether it is needed or not. The

resulting linear matrix equations are banded, however, and can be easily solved on a high speed

digital computer. The major difficulty with this method is the complexity required during the

remeshing of the model with each time step after the material removal from the anode.

On the other hand, the BEM formulation does not require an internal mesh Ref.[13]. It requires only a mesh on the boundary which is an advantage when applied to the ECM process

since we are chiefly interested in the solutions on the surface. In addition, since the ECM process

model requires a moving cathode(moving boundary), the application of the BEM method provides a better numerical model since remeshing is only required on surfaces to simulate gap closure and

workpiece dissolution on a surface. In the case of using FEM, remeshing at each time step would

require updating the mesh in the interior as well as on the boundary.

Solution of the integral equations for the BEM formulation is accomplished by using a fundamental singular solution(Green's function) and applying Gauss's theorems of potential theory

that allow converting the problem from the solution of any interior domain to the distributed

contribution of the potential and its normal derivative on the bounding surface. The application of

the singular solution and the derivation of the boundary integral equations used in this formulation

are contained in the appendix. The numerical solution of these boundary integral equations was

developed using the method developed by Brebbia. For a more thorough exposition on this topic,

the reader is referred to chapter 2 of Ref.[151. His formulation and approach was used as the

building block for developing the electrochemical machining process model for simulating the

ECM'ing of gas turbine airfoils and blades. The most important variable that must be determined

in the ECM process is the current density distribution on the workpiece with the subsequent

application of Faraday's Laws of electrochemistry. Consider the airfoil ECM application shown in

Figure (2). If a cutting plane is passed through the tooling and the blade normal to the stacking

axis, a 2D model at that section is obtained. Figure (3) shows at the left a CAD/CAM drawing of

87
an airfoil blade with the dual opposed cathode electrodes and the equivalent BEM model on the

right.

The electrolyte conducts the electric current between the cathode and anode by an ionic

transport phenomenon which provides the mechanism for surface dissolution of the airfoil(anode).

This dissolution of the anode surface produces a slowly moving boundary condition in the

mathematical formulation. Crank [16] rigorously derives the moving boundary conditions which

must be incorporated into a precise formulation of the ECM process. In the case of a 2D problem,

the current density distribution produced on the anode will not be uniform in this process. It is the

purpose of this model to accurately determine the current density distribution by the BEM method.

The electrolyte flow in the gaps between the cathode and anode is a high velocity turbulent

flow. The electrolyte flow will split at the inlet for airfoil machining and merge at the outlet. In the

present development, the universally accepted velocity profile for turbulent channel flow is assumed for the electrolyte fluid mechanics. The basic equations used to model the electrolyte

flow characteristics are contained in the appendix. The most important consideration is the

temperature rise that occurs in the machining gaps during the machining process from inlet to

outiet. The temperature rise in the electrolyte gap is due to Joule heating and is proportional to the

local current density. The mrbulence in the flow produces essentially a uniform temperature across

the gap which varies inversely as the gap size but increases in the flow direction. The wall shear

stress resistance due to the viscosity of the electrolyte follows the widely accepted Darcy friction

factor method which is dependent on the Reynolds number of the flow. The electrolyte properties

such as electrical conductivity, viscosity, density and specific heat are a function of temperature and

must be accounted for in order to have a realistic simulation.

Using the Boundary Integral Equation method, many complicated geometries can be

conveniently represented by a simple 2D contour of the ECM process model consisting of the

boundaries enclosing the tool, electrolyte gap, insulation and workpiece. The advantage of this

formulation is that only the boundary composed of the tool, the workpiece and insulators with the

interelectrode electrolyte have to be discretized to construct a 2D model without the necessity of

discretizing the interior of the electrolyte between the tool and workpiece as would be the case if the

finite element method were employed. This analytical modeling method is most appealing since in

the electrochemical machining process we are only interested in the reactions and dissolution of the

anode at the surface. Another reason for using the boundary element method is that greater

numeral accuracy can be achieved on the surfaces in the model since the functions and their

derivatives are calculated directiy at the surfaces. In the finite element formulation, the functions

and their derivatives are calculated at interior Gauss points and the values at surfaces are extrapolated from these interior points. The kinematics of the model is also easier to implement

for both the velocity of the surface mesh of the tools and the surface removal of the workpiece. It turns out that an initial coarse surface mesh on the airfoil becomes a finer mesh as material is removed since the distance between element nodes at the surface is reduced. At each time step,

the current density distribution and the nodal displacement due to machining the workpiece, the

total accumulated weight removed and tool displacement are calculated in sequence until the final

process time is reached.

Electrolyte Properties

For an effective process model, it is necessary to incorporate electrolyte properties into the

ECM computer model to be able to simulate the electrochemical machining behavior at the

surface. The electrochemical surface dissolution of the anode is a complicated phenomenon and is different for various cell combinations of metal and electrolyte. The properties of the electrolyte are chiefly dependent on the concentration of the solute and temperature.

88
The ECM2D computer code has a materials properties library for both electrolyte and

workpiece material properties which were determined in the laboratory by experiments. The

material library requires electrolyte data as a function of temperature, concentration, and current

density to determine the electrolyte's electrical and thermal conductivity, viscosity, density,

oveipotential and cutting efficiency. It is also essential to include the overvoltage in the model for

the precise determination of the local current density distribution on the surface of the workpiece.

Applications

The 2D process model computer program had been developed for simulating the ECM

process to predict the workpiece shape change including the total accumulated weight removed by

using boundary element method formulation. Provisions for cross-flowing electrolyte have also

been incorporated into this model. When this option is exercised, electrolytic heating is accounted

for based on the pressure differential across the electrolyte inlet and outlet. The model predicts the

surface current density distributions on both the tool and workpiece, final machined shape of the

workpiece, total current consumed, and total weight machined from the workpiece.

EC Airfoil Machining Simulation: Dual Electrode

As an example, the results of modeling a typical fan blade in the aircraft engine industry

is described. Actual machine and process parameter settings were used in the computer program

to analyze the machined shapes, weight of material removed, current density distribution around

the perimeter of the airfoil, interelectrode temperature increase in the electrolyte, and the current

history used in the process. Since this code is a 2D computer program, the full height of the blade

was modeled by several cutting planes normal to the stacking axis.

Figure (4) is an illustration of a BEM model that was constructed to simulate the 2D ECM

process for machining a typical cross section of fan blade near the root. This figure shows the

nodal points for the cathode tooling surrounding a forged airfoil preform at t= 0.0 sec. Note that

only a surface mesh is required to define the model. This model was oriented such that the trailing

edge was at the electrolyte inlet on the left and the electrolyte exit at the leading edge on the right.

The contour of tiie model is composed of series of elements on the surface.

Boundary conditions are specified on the entire model. For example, voltage potentials for

each element are specified on both the cathode and anode. Elements used to represent insulation or

the inlet and outiet channels of the electrolyte are specified by setting their normal derivatives to zero, tiius indicating no flow of current. During the machining process, the cathode tools have a machine prescribed velocity and voltage as shown in Fig.(5). The velocity and voltage are set

initially high, but are reduced as the machining gap becomes smaller.

As the machining process progresses in time, only the surface mesh for the the tool and

insulation remain the same. The relative position of the tooling mesh changes only if a feedrate is applied to the tools. Remeshing occurs over the entire surface of the workpiece at each time step

due to anodic dissolution. Workpiece remeshing is required during the simulation to maintain

optimal nodal spacing over the surface to prevent nodal crossover and improve the leading and

trailing edge definition. Figure (6) shows the location of the cathodes and the final machined

shape after 215 seconds of simulated process time. Note that the interelectrode gap at the end of

the process is very small and is approximately 0.010 inch.

Fig.(7) also shows the final position of the tooling and the machined airfoil near the tip at t= 215. seconds. With the zoom capability of the computer program, the analyst is able to investigate tiie geometry of the airfoil at the trailing and leading edges. As an illustration, Fig.(8)

shows the shape of the airfoil at the trailing edge. The definition and accuracy of die simulated

profile depend on the number of elements in the model, with a greater number yielding superior

89
results. The rate of dissolution at the airfoils leading and trailing edges limits the number of nodes

that can applied in practice. Higher nodal density requires smaller time steps to prevent nodal

cross-over on the workpiece surface. The analyst must make the final determination of ultimate

model resolution verses the cost of running a larger model with smaller time steps.

An accurate determination of the current density distribution around the outer contour of

the airfoil is extremely important to precisely simulate the ECM process. Figure (9) is a plot of the

current density distribution around the airfoil at three different times during the process. At

t=212.0 seconds in the process, the cathodes are close to the airfoil and produce the final shape.

The current densities are largest at both the leading and trailing edges, which is due solely to geometric effects. The current density on the convex and concave sides of the airfoil are essentially

uniform.

One of the key results of the analysis is the difference in the total temperamre rise in the

gap. While the maximum temperature rise near the root was 12.86 ^F, near the tip the

temperature rise was slightly less, being equal to 1 1.88 ^F. This temperature rise is directly

proportional to the square of the current density and the flow rate of the electrolyte in the gap. The

reduced flow rate at the root can be attributed to a narrower interelectrode gap. This temperature

rise at the outlet is small due to turbulent mixing when the flows merge. There is also a slight

difference in the electrolyte flow rates on the convex and concave sides due to the split flow that occurs at the inlet and the initial position of the preform. If there is a significant difference in the

balance of the flow through the two gaps, there can be a significant overall effect on the workpiece

shape.

Conclusions

The boundary element method(BEM) is an effective analytic tool that can be used as an ECM process simulator for predicting the final machined shape of 2D workpiece geometries

based on specified materials and process parameters. It has also been shown the BEM technique can be used to successfully predict the shape of an aircraft engine airfoil or bucket, considered to be one of the most difficult surfaces to model due the complex surface geometry and sharp radii at both leading and trailing edges.

ECM modeling offers unique insights into the process that can only be effectively obtained

from modehng. This includes such items as the interelectrode gap size, current density

distribution, and the electrolyte pressure, temperature and velocity distribution. Models can be

used to conduct parametric smdies of future tool designs and to aid with the interpretation of

experimental results. In addition, the ECM simulator furnishes an alternative rational way to establish the initial process parameters during the process development stage. In practice during electrochemical machining only external parameters such as inlet and

oudet pressures' applied voltage and total electric current can be readily measured and monitored.

A process model provides a tool to determine the state of machining in the interelectrode gap

which is of major interest. This allows the user to investigate the electrolyte temperamre, velocity,

current density on the workpiece and cathode tooling, pressure distributions and gap size along the

flow path. These variables are virtually impossible to measure except by laboratory methods. For

example, consider the current density distribution on the anode workpiece. This single variable has

a significant effect on the local machining mechanism by virtue of Faraday's laws of electrolysis.

Being able to predict the current density distribution would allow the user to make the necessary

changes in tool design and to adjust the tool velocities for better performance and surface finish.

90
References

[I] A.E. De Ban- and D.A. Oliver, Electrochemical Machining, American Elsevier Publishing

Co. Inc., New York, 1968.

[2] J.A. McGeough, Principles of Electrochemical Machining, William Clowes & Sons,

London, 1974.

[3] R.H. Nilson and Y.G. Tsuei, Inverted Cauchy problem for the Laplace equation in

engineering design ", Jour of Eng Math, Vol.8, No.4, October 1974, pp. 329-337.

[4] R.H. Nilson and Y.G. Tsuei, "^Tree Boundary Problem for the Laplace Equation with

Application to ECM Tool Design", Transactions of the ASME, March 1976, pp. 54-58.

[5] Y.G. Tsuei, C.H. Yen and R.H. Nilson, 'Theoretical and Experimental Study of Workpiece

Geometry in Electrochemical Maching", ASME publication, 76-WA/Prod-6, 1976, pp. 1-5.

[6] R.H. Nilson and Y.G. Tsuei, 'Tree Boundary Problem of ECM by Alternating-Field

Technique on Invened Plane", Computer Methods in Applied Mechanics and Eng, 6(1976),

pp. 265-282.

[7] P. Forsyth, Jr. and H. Rasmussen, "A Kantorvich Method of Solution of Time Dependetn

Electrochemical Machining Problems Computer Methods in Applied Mechanics and Eng.

23 1980, 129-141.

[8] V.K. Jain and P.C. Pandey, ''Design and analysis of ECM TooHngs", Precision Engineering,

Vol 1, No.4,Oct. 1979, pp. 199-206.

[9] V.K. Jain and P.C. Pandey, "Finite element approach to the two dimensional analysis of

electrochemical machining". Precision Engineering, Vol 2, No. 1 , Jan. 1980, pp. 23-28.

[10] V.K. Jain and P.C. Pandey, "Tooling design for ecm", Precision Engineering, Vol 2,

No.4,Oct. 1980, pp. 195-206.

[II] V.K. Jain and P.C. Pandey, "Tooling Design for ECM- A Finite Element Approach",

ASME Jour, of Eng. for Industry Vol 103, May 1981, pp. 183-191.

[12] V.K. Jain and P.C. Pandey, "Tooling design for Electrochemical Machining of Complex

Shaped Workpieces", J. Inst. Eng. India, Vol 62, Nov. 1981, pp. 95-99.

[13] O.H. Narayanan, S. Hindua and C.F.Noble , "The Prediction of Workpiece Shape During

Electrochemical Machining by the Boundary Element Method", Int J. Mach. Tool Des.

Res., Vol.26 No.3, 1986, pp. 323-338.

[14] O.H. Narayanan, S. Hindua and C.F.Noble , "Design of tools for electrochemical machining

by the boundary element method", Proc. Instn Mech Engrs, Vol.200 No. C3, 1986, pp. 195-

204.

[15] C.A. Brebbia & S. Walker, Boundary Element Techniques in Engineering NewnesButterworths, London, 1980.

[16] John Crank, Free and moving boundary problems. Clarendon Press, Oxford, 1984.

[17] G.F. Hewitt and N.S. Hall-Taylor, Annular Two-Phase Flow, Pergamon Press, Oxford,

1970.

Appendices

(1) Integral Equation Formulation of Boundary Value Problem : Since the ECM process is essentially a surface phenomenon which requires a precise

determination of the 2D current density distribution on the workpiece, the boundary element

method(BEM) was selected as the preferred numerical method for solving the governing

differential equations. The boundary element method was found to be a powerful computational

tool for determining 2D shape changes during the machining process. As previously stated, one

of the major features of the BEM technique is that only a mesh on the outer contour or on an enclosed surface is required for constructing a model. Since one is only interested in surface

91
dissolution, it would not be necessary to discretize the interior of the model which is the case if finite elements were to be used. By virtue of the surface formulation, ie, current density and

potential, the model will provide a more accurate representation than if the finite element method

was employed for this particular manufacturing process.

The basic concept of boundary integral equations has its roots in potential theory.

Referring to Fig.(16) consider a material bounded by a surface. On this surface either the potential or its normal derivative must be specified. The solution objectives for such boundaiy value

problems is to seek the value of the potential and its first derivative everywhere in the interior and on the bounding surface. On the surface where a potential has been specified as a boundary

condition, the first derivative is sought. On the boundary where the derivative or flux has been

specified, the value of the potential is sought.

The direct method for the formulation of the mixed boundary value problem is based on

Green's identity. Referring to Fig.(16), consider the potential problem written in terms of the

Laplace equation:

V^<D=0 in domain n (1.1)

where 0 is the potential. The mixed boundary conditions for the specification of the potential and

the flux can be written as:

<I) = 0 on the surface

q = —= (7 on the surface r, (1.2)

dn

where the total boundary is given by: F = Fj + Fj. The notation of the bar denotes the known

boundary condition specified. Consider a weighting function <£>* which has continuous first derivatives. This function will be required to satisfy the govenering equation. Using the weighted

residual method, the integral equation based on Greens 's theorems can be written as:

J(V2<I))0*dn = J(q

- q)0'dF - J(cD - 0)q*dF (1.3)

' n r,

where q = — and q = — . The first integral is a volume interal, whereas, the two integrals on

on on

the right hand side are surface integrals. When the integral on the left side is integrated by parts, an

integral equation can be obtained which will be the starting point for the boundary element method.

Accordingly, (1.3) can be written as: Jo(W)dn = -JqO*dF- Jq<I)*dF +

J<I>q'dF + Joq*dF (1.4)

a r, r, r,

(2) Boundarv Element Formulation of Boundarv Value Problem : The integral equation in appendix (1.4) can be recast into a numerical form by introducing a fundamental solution which satisfies the governing equation and the associated boundary

conditions. In the case of an electrostatic problem, if we assume a concentrated charge is acting at point "i", the original governing equation can be written as:

' V'O +A^=0 (2.1)

92
where A' is a Dirac delta function. The solution of this equation is the fundamental solution. If

eqn.(2.1) is satisfied by a fundamental solutions, then

J

<D(V'a)*)dQ = -0' (2.2)

a and represents the unknown function at point "i". When this relationship is applied to eqn

(1.4), the following equation will be obtained.

& + J<Dq*dr +

J

<Dq*dr = JqO'dF + JqO'dT (2.3)

where the fluxes are given by q = — and q = —

dn dn

The integral equation solution can be obtained by using the fundamental singular solution for the

2D domain.

(3) Integral Equation Formulation of Boundary Value Problem : Using Brebbia's method Ref.(4), the integral eqation can be rewritten in the following form

by applying the singular solution and taking limits. Accordingly, the integral equation for the

potential at any point on the surface can be written as: io' + Jo^ dr + Jo^ dr = Jq<D*dr

+

JqO*dr(3.1)

Tj r, fj r, Equation (3.1) can be generallized in a compact form for a later numerical format.

io' + Jo q*dr = JqO'dr (3.2)

r r—

where r = Fi + and <I> = 0 on Fjand q = —= ^ on Fj. on The final desired numerical form for application in the ECM code is obtained as:

The integrals in (3.3) are evaluated using the Gauss quadrature method.

(4) Basic Equations for Electrolvte Flow Formulation and Temperature Rise: The derivation of the integral form for tiie momentum equation is derived in Ref.[17]

which is essentially a statement of the balance of die pressure, wall friction and momentum forces

93
along a control volume on a streamline assuming that the gravitational effects are negligible from

inlet to outlet. This equation is written as:—dP f C d

^sdA= I x^^sdS+ —(mv)AsdA (4.1)

ds JJ ds

s where s is the distance along the streamhne, is the wall shear stress, m is the mass flux rate, v

is the velocity along the streamline, A is the cross sectional area and S is the surface area of the

stream tube.

This relationship can be reduced to a differential form that can be employed in the

computer code:

-A/' = -^Aj + pvAv (4.2)

where //^^ is the hydraulic radius and p is the mass density of the electrolyte. In this form the

units on each side of the equation are gm

The other relationships necessary for the flow formulation are the wall shear stress , the dimensionless wall friction factor for mrbulent flow X and the Reynolds number

respectively. Accordingly, they can be written as: K=\W.. (4.3)

where X is the widely known Darcy friction factor. These relationships have not been derived

from first principles, but have been validated by extensive experimentation. For a highly turbulent

flow this friction factor becomes:

X = 0A95{Log,,R„f' (4.4)

which is dependent on the Reynold's number of the flow.

K = ^

(4.5)

/J. with /J. being the electrolyte viscosity. In this formulation all the electrolyte properties are temperature dependent.

The temperature rise in the electrolyte can be determined from the basic energy equation.

In Ref.[15] , it was shown that the heat conduction term for the high Reynolds number application

plays an insignificant role compared to the convective term. The following equation results: P.c,v^ = yV, (4.6)

where the subscripts refer to the electrolyte property. The density, specific heat and specific

resistivity are given by the symbols , and cr^ , respectively. The temperature rise in the

electrolyte is essentially due to the local Joule heating where / is the cuirent density. This

equation can be rewritten in a numerical form for programming as: AT =

PeCe)

•As (4.7)

94
Figures

Tool

Electric

current

^

Cathode ~~

flow

boundary layer region_o__o_o _C) 2_ Q-Q—O-lJ

Hydrogen gas generated

boundary layer region ^7"jiP — *^ T "

dissolution

Workpiece

Anode +

Fig.(l) Schematic of the Electrochemical Machining (ECM) Process

TOOL

FEED

BRASS

ELECTRODE

ELECTROLYTE

FLOW

INLET

OUTLET

BRASS

ELECTRODE

Fig.(2) Airfoil ECM process tooling configuration

95
Fig.(3) CAD/CAM drawing of airfoil and tooling.

Tooling/Workpiece initial conditions

6 j——1 • 1 • 1 • 1 • 1 5

4

• 3

2

-24 •

-5_g I. I -6 -4 =2 0 2 4 6

X(cm)

Fig.(4) BEM representation of tooling and airfoil preform at t=0.0 sec.

96
Velocity/Voltage Vs Machining Time

0.14

0.12

O 0.10 w

c

0.08 u

O 0.06 >

0.04

0.02

I I I

I

I

I

f

It--^

I

I I

I I

'

V«locity (ir/sec) Voltagt

-r

I

I

I

I 3=^

-15

14

-13

Itage

-12 Vol

-11

-10

0 25 50 75 100 125 150 175 200 225 250

Machining Time (sec)

Fig. (5) Cathode velocity and voltage as a function of time

Final Tooling/Airfoil Position

E

u

2.0

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

-2.0

1 1 I 1 1

Tool Position

4 Airfoil Position

n

X(cm)

Fig.(6) BEM model of Tooling and airfoil near the root at t= 215.0 sec.

97
Final position of tooling/airfoil

1.00 0.75

0.50

0.25 ? 0.00

-0.25 • >•

-0.50 •

-0.75 •

-1.00

-1.25 •

-1.50

Fig.(7) BEM model of Tooling and airfoil near the tip at t= 215.0 sec.

1.00

Final position ot tool Ing/alrtol

I

Fig.(8) Zoom view of BEM model of Tooling and airfoil at t= 215.0 sec.

98
Current Density Distribution Comparison

Section EE (near root)

t- too 3CC3. B t- 1 89 3ec3. S t-212jecs

-80-

Numbers

Fig. (9) Current density distribution over airfoil at different times

q =

3n

X

X,

Fig(lO) Boundary Element Method (BEM) concept

99

Boundary Elements and a Green's Function Library

F.J. Rizzol, P. A. Martin^, and R.A. Roberts^

1 AEEM Dep. Iowa State University, Ames lA 5001

^ Dept. of Mathematics, University of Manchester, Manchester Ml 3 9PL, U.K.

3 Center for NDE, Iowa State University, Ames lA 5001

1 Introduction

The boundary element method (BEM) is now an established procedure

for obtaining numerical solutions for a variety of problems in engineering

and applied mathematics. The formulation of the BEM relies heavily on the

existence of suitable Green's functions. Indeed, Green's functions are the

BEM's main analytical ingredient.

With the conventional BEM, a significant analytical step is taken at the

outset, and this involves only the simplest Green's functions. The result is a representation integral for desired fields in terms of boundary values of the

fields. This result is obtained before elements and approximations of any kind

are introduced. Unfortunately, when only the simplest Green's functions are

used, about half of the boundary values of the fields are unknown in the

representation integral at this stage. Elements and approximations are needed

afterwards, in essence, to numerically solve a boundary integral equation.

This is done to obtain the mentioned unknown boundary data. Then, with all boundary data known, the representation integral provides the desired field

solution throughout the region of interest. The process, in effect, reduces a three-dimensional problem to a two-dimensional one. This is one of the great

features of the BEM.

Now if more sophisticated region-specific Green's functions were to be

used at the outset, less unknown boundary data would appear in the

representation integral, fewer elements and associated approximations would

be needed with the BEM, and accuracy could be increased while computing

demands would be reduced. Carrying such reasoning to its end, one would

need no elements at all if exact, region-specific Green's functions could be

found. No unknown boundary data would appear in the representation

integral in such cases, and they would provide the desired fields everywhere

after the first step. Errors for real problems would be limited to numerical

quadrature errors on integrals with known integrands.

The simplest Green's functions (G) which give fields due to point

disturbances in all of space, are well-known in analytically-convenient,

closed form for broad classes of problems. The more advantageous, more

sophisticated Green's functions (G*), which give such fields in the presence

of bounding surfaces and other problem-specific features, have not been

known, in any form, except for relatively few simple geometries and boundary

101
conditions. Indeed, this is perhaps the main reason why the conventional BEM

is based on the simplest functions, despite considerable computational effort in

using it for industrial problems.

However, there has been some recent success in constructing more

sophisticated Green's functions G*, in analytical or approximate-analytical

form, for specific shapes S, for a variety of classes of problems e.g.

[ 1 ]. Therefore, one might try to preserve many of the mentioned advantages

with the BEM, when G* is known analytically, by first assembling a 'library' of

such existing G*. Then simplified computer codes, which could be nothing

more than quadrature routines involving prescribed boundary data and the G*

for specific S, could be written by users to solve their own real problems. In

fact, many such codes already exist as "the last step" of existing BEM codes.

With a library of G*, the last step is all that would be required for many

problems. Successful experience with such a library could be motivation to

construct more and more sophisticated G* entries over time.

Despite this scenario, it is clear that many classes of problems, for many

geometrical shapes, will probably defy, indefinitely, the construction of G* in

even an approximate-analytical form. For such cases it is possible to construct

G* in discretized or numerical form, for a variety of difficult but commonly

ocurring surfaces S. These G* could be placed in another type of library.

Such library entries would take time and computational effort to create.

However, like the creation of analytical G*, these would be one-time tasks.

With modern technology for storage and quick retrieval of massive amounts of

data, e.g., on compact disk, and transmission of such data on computer

networks, it is possible to take a fresh look at the advantages of using even

numerical G*, versus conventional BEM which uses the simpler analytical G.

How to construct a library of discretized G*, which would consist of

matrices of numbers as its main ingredient, and how to use the library, are

less clear than with analytical G*; but these matters are the subject of this

paper. In fact, the conventional BEM is a prime vehicle for constructing

discretized G*. We intend to show that use of discretized G*, with some

attention to standardized protocol, could be almost as convenient and accurate

for problem solving as having their analytical or approximate-analytical

counterparts. Speed might be improved too. To obtain discretized G* for the

two-surface problem, as described below, requires a number of matrix

multiplications. This can take some time. However, existing G* in analytical

form often require extensive integrations or series summations, and these are notorious in the amount of time and effort to get numerical values from them.

The focus of the rest of this paper is to look in some detail at the

consequences of getting more sophisticated Green's functions G* into the BEM

picture. To fix ideas, we outline the essential aspects of creating a library of

discretized or numerical Green's functions for problems of scattering of timeharmonic acoustic fields by one or more three-dimensional obstacles, each

bounded by a smooth surface S. Ideas for one and two scatterers are treated

explicitly. How to deal with more than two scatterers will be obvious, and

simplifications which occur whenever G* is known in at least approximateanalytical form will be apparent. How to proceed for problems other than

ones in acoustics, such as elastostatics, time-harmonic elastodynamics, steady

state heat conduction, or any problems governed by linear elliptic partial

differential equations, should be apparent as well.

102
Nevertheless, the NIST workshop, which gave rise to these proceedings,

was concerned with Green's functions and boundary elements with

applications for modeling the mechanical behavior of advanced

materials. Thus we close this introduction with some comments on the kind

and character of G* useful for this purpose.

Specifically, for materials-related problems, we wish to create for the

library as many Greens functions G* as possible, which contain the

geometrical and constitutive features of the most important advanced

materials. Fields of greatest interest are likely to be elastostatic or elastodynamic, when modeling mechanical behavior; but again, acoustic,

thermal, or even electromagnetic responses to certain inputs may be relevant

to such behavior.

A typical strategy would be to model generic or model problems for

"cells" of composite materials, first for materials in their perfectly-bonded,

undamaged state. Then, the G*, in either analytical or discretized form, for

these problems could be formed and stored in the library, ready for nominal,

representative, static loads to be applied by materials analysts. From these

loads and library of G*, a variety of responses could be easily calculated from

which, in turn, a "cell" stiffness or modulus could be constructed.

Macrostructural behavior of bodies made of such composites could thereby be

assessed.

Subsequent models could include damaged materials, i.e., those with

common disbonds in advanced composites or cracks which grow in

characteristic patterns. Then, library entries G*, with the essential, difficult,

geometrical features of these models, could be created, ready for scientists to assess the influence of the damage on local stiffness and thus, subsequently,

on the macrostructural behavior of bodies containing such defects. In any

case, with a library of proper G*, the desired fields may be generated by

representation integrals similar to the ones discussed above. The value of

"what-if' experiments, that could be quickly and easily run with a good

library of G* for complex materials, damaged or undamaged, would evidently

be considerable.

2 The one-surface problem

Consider time-harmonic scattering of acoustic waves in 3-D by a bounded body B' with surface S. The representation integral for the acoustic

field u at a point P in B exterior to B' is

where G(P,Q) = G(aP) = - eikR/(2:7r R), R is | Q - P|, with point aalso in B and

points q (and p et. seq.) on S, k is a frequency parameter, and u satisfies a

radiation condition at infinity. Representation (1) is obtained by applymg

(1)

103
Green's theorem to u and G in in the usual fashion [ 2 ]. Now if du(q)/dn = f(q)

is prescribed on S, expression (1) does not give the solution for u(P) since u(q)

is unknown.

However, suppose G* = G + w, where w is a regular function (satisfying

the same governing differential equation as u), can be found, and G* used in

place of G, such that

G*(^,F) = 0 forqeS, (2)

then (1) reduces to the integral

2u{P) =

( f{q)G*{q,P)dS^ ^ S'J (3)

which gives, in fact, the solution u(P) or u(p). With G* instead of G,

everything under the integral sign is known.

Thus with known G* for given S and k, it is clear that solutions for

arbitrary / are obtainable with a simple quadrature. Having a library of G*'s,

for as many shapes S as possible, would therefore have obvious advantages.

The result is

(4)

(5)

Therefore in light of (3) and (5) it is true that on S

' A'% = ^s;, (6)

such that solving the boundary integral equation (BIE) (4), is formally

equivalent to finding the Green's function G* (cf. [3 ], and [ 4 ] eqs. (28), (29)).

Pursuing this line of reasoning a bit further, it can be shown [ 2 ], by

applying Green's theorem to G and G*, that G* satisfies the BIE (4), with 2G in

place of the integral on the right hand side of (4). Since this is true, it follows

that

Alternatively, consider the limit as P goes to p in (1).

u(p)+-Cu{q)-^(q,p)dS^ I u{q)—G{q,p)dS^ = Irjf{q)G{q,p)dS^

s s

or in operator form

Au^SJ.

104
G''{q,P)^2A''G{q,P), (7)

where it is important to note that relation (7) between G and G* holds so long

as at least one point in the argument of each function is on S. Formally

inserting expression (7) for G* into (3) we obtain

«(^)= Cf(q)A-'G{q,F)dS{q). (8)

If A Ms assumed known, (8) like (3), represents the solution u(P). Equation

(8) is hardly a new result. Nevertheless, its use as considered below does

appear to be new.

With existing boundary element methodology, it is possible to regard (8) as having many of the advantages of (3). Specifically, it is possible to first

form a nodal-value-approximation to A'\ say in the form of a (NxN,

(N=number of nodes)) matrix, for as many shapes S and frequencies k as desired, using a good robust BEM. Next, multiply A^' by G evaluated at nodes q

for chosen P. We now have a nodal approximation to G*/2. Now improve this

approximation over the boundary elements using appropriate shape functions

in the variable q. Similarly represent / , so that it is now possible to integrate

the only remaining variables under the integral in (8), namely, products of

shape functions. The result of this integration is another square (NxN) matrix

K, similar in character to the stiffness matrix in finite elements. With this

process, an approximate form of (8) may be written

(sum on i,j = 1... N )

u(P) = f(l^)A^a,q^)G(q^,P) (9)

where /(/,) is a (IxN) row matrix of nodal values of /, G{qj,P) is the (Nxl)

column matrix of nodal values of G for desired (parameter) P, and A^\l^,qj) is the product of A^^ with K.

The ingredients in (9) and the strategy surrounding their formation

and use deserve more discussion.

Suppose a library of A^^ were available for a sequence of (say oblate)

spheroidal-shaped rigid scatterers in an acoustic medium. Let the entries be

specified by an eccentricity parameter e and frequency parameter k. Now a

library user, with known /'s in hand (i.e. known input waves at a certain k),

would like to know the scattered field at desired P, from a rigid spheroid of

certain e . This user could proceed as follows.

Locate the proper A^^ for chosen e and k, specify / in a standard (easy)

format, and specify a list of specific P locations for the desired fields. The

software to pick the necessary A^^ from the library, multiply by K to get A^'

(i.e. do the integration in (8)), and finally do the multipHcations in (9) with

105
the row / and column of G for chosen P, could all be part of a black box. The

box itself could be part of the library. The point is, with a good library,

accurate reliable u(P) values could be obtained rather quickly with little or no

knowledge of the underlying process required from the user.

In fact, if a quadrature scheme (order of shape functions, etc.) could be

decided upon in advance, it would be possible to store A~' rather than A~\ and

save some execution time. This and other such issues should be transparent to

the user. In any case, interested knowledgeable users, who might wish to

write some of their own library-access software, are faced primarily with

tasks involving formation and multiplication of matrices. The main, complex,

time-consuming task of getting A^^ would already have been done.

3 The two-surface problem

In this section we consider the same scattering problem as in section 2

but with two bounded scatterers. For now, consider their surfaces. Si and S2 to

be disjoint; other configurations will be treated later. Here, we assume the

desired scattered fields satisfy a radiation condition as before and also satisfy

the boundary conditions

^=fj 7 = 1,2 (10)

dn

where /, and are given functions on Sj.

Proceeding as in section 2, the counterpart of equation (1) for the two

surface problem is

Again, (11) contains unknown data on the surfaces, namely, m(^i) and uiq^) on

Siand S2, respectively.

However if G* for S2 were known and used in place of G in (11), i.e.,

where

^ G*{q„P)=0, (12)

then (11) reduces to

106
2«(F)=

J

[fM.)GHq^^P)-u(q.)j^G^{q„P)]dS^ +Jf,{q,)G*{q„P)dS^ (13)

which is the two-surface counterpart of (3) with G* for S2. Note the

dependence on the unknown function M(^2)is missing, such that the integral

over S2 is known. Of course unknown on Si, is still present.

To simplify the subsequent discussion, assume further that is zero.

(It is a simple matter conceptually to add this integral back for nonzero , and

dealing with this term is no more difficult than with (3) via (9)). Without the

last integral, (13) is formally the same as (1), with G* in place of G.

One might be tempted now to try to replace G* in (13) with a function

G**, say, where the normal derivative of G** vanishes, not only on S2, as does

G*, but also on Si. With such a G**, an equation like (3) could be written.

Then, in principle, the strategy described above following (3) would pertain.

However, there is difficulty enough in trying to find analytical G* or numerical G* as described for the one-surface problem, such that a comparable strategy for a G** is best postponed, perhaps indefinitely.

Nevertheless, it is worthwhile to view (13) in a fashion similar to (1),

with (4), and (5). That is, imagine solving a BIE like (5) for u on Si where only

Si needs be discretized with boundary elements. This is possible, without

discretizing S2, if G* is used in the process. Information about S2 is contained

in G*. Thus, with G*, the two surface problem, via the BIE/BEM, is formally no

more difficult than the one surface problem.

If one must use a numerical G* in (13), many of the issues already

addressed in connection with (7), (8), and (9) are still applicable. However,

there are some new issues as well.

Specifically, note first that the arguments of G* and its normal

derivative as appear in (13) involve qi as well as P. Thus any BIE/BEM

methodology involving Si, and subsequent use of (13) as a solution for u(P),

requires G* as a two-point function, where neither point is on S2. Thus (7) is insufficient for this purpose.

In light of this, consider another consequence of applying Green's

theorem to G and G*, namely (cf. [ 5 ], [ 6 ]),

2G*((2,P)=2G(P,0- CG*(l,P)^G{l,Q)dS, (14)

where S in (14) should be regarded as S2. Now since / in the argument of G* in

(14) is on S2, we may use (7) to write (14) as

107
G*{Q,P) = G{Q,P)-r-^G(l,Q)A-'G(l,P)dS,, (15)

where Q_ and P may be interchanged in any of these expressions, since both G

and G* are symmetric in these variables. Similarly, (cf. [ 2 ]) it is true that

^ G*(q„P) = -^G(q„P)-r-^G{P,l)A-'^G(l,q,)dS^ (16)

where again integrals over S mean over S2.

Now with (15) and (16), and existing boundary element methodology,

the reasoning leading to a library of information about S2 is similar to that

surrounding (8). Specifically, it is possible to write (15) as

G^{Q,P) = GiQ^P) -G"'{Q,l,)A],\l,,q.)G{q.,P) (17)

where the n superscript means normal derivative at / , and l^ and q. are the

arguments for row and column matrices as before, andA^^(/,,^p, is exactly the

same matrix encountered earlier. A similar expression exists for the normal

derivatives in (16), namely,

? G'^"' {Q,P) = G"HQ.P)-G"'{P,l^)A-^\l^,q^)G"'{q^,Q) (18)

Note that a row-times-square-times-column multiplication is required

for each P, Q_ choice in expressions (17) and (18). The results of such

operations are the G* function evaluations as needed in a BIE/BEM treatment

of the surface Si. Other aspects of the treatment are the same as if the freespace G were usable, i.e., as if S2 were not present.

Observe that expressions (15) and (17) have the classic form for a region-dependent Green's function, i.e., G* = G + w, where (minus) the integral

term in (15) is w, and (minus) the triple-matrix product in (17) is (an

approximate) w. Moreover, each expression for G* can be interpreted as "the

field at Q.due to a point disturbance at P in the presence of a surface (S=S2) on

which the normal derivative of the field vanishes". This is precisely the

interpretation of (one of a class of) region dependent Green's functions.

4 Partitioning

Consider equation (11) again (with f2 zero), and suppose limits are taken as P goes to pi on Si and p2 on S2, respectively. The result is

\,u, + \^u^ = S,J, (19)

108
where the operator notation of (5) is invoked with the following additional

considerations: the first subscript on the A and S operators refers to pi

locations on Si, whereas the second subscript on those operators refers to qj

locations on Sj; the single subscript on u and f refers to qj locations on Sj; for

the Aij operators with i not equal to j, the free u(p) term (cf. (4)) is zero. Now if U2 is formally eliminated from (19) and (20), the result is

(Al ~ Al-A^lAl)^! ("^11 ~ Az-Az'^Zl)/! • (21 )

Next, if the group of terms in parentheses on the left side of (21) is called A*,

and the group of terms in parentheses on the right side is called S*, (21) is formally the same as (5). Thus the presence of S2 is manifest in the operatortriple-products with A^l in their centers.

In light of the observations about (21), made possible by the

partitioning process, and the form of equations (1) and (13) [with zero], and

the form of equations (15) and (16) [or (17) and (18)], it probably occurs to the

reader to question whether the partitioning process and the process of

defining and using G* are equivalent. Indeed, this is the case as is rigorously

shown in [ 2 ],

Specifically, if G* and its normal derivative, as given by (15) and (16), are used in place of G and its normal derivative in (4), the "n" terms in (21)

come from the "G" part of the Green's function, and the triple-product terms

come from the "w" (integral (over S2)) part.

Note finally that the surfaces Si and S2 above need not be disjoint. The

partitioning of a single surface into two parts is arbitrary and really a matter

of convenience. For example, S2 may be a common-shaped appendage on a variety of problem-specific shapes Si. Specifically, S2 may be a thin antenna

attached to the surfaces Si of different vehicles. Most of the discussion in

sections 3 and 4 applies in such cases.

5 Discussion

Consider again the work involved in the creation and use of the Green's

function library. For clarity, consider these matters for the one- and twosurface pi-oblems separately.

For the one-surface problem, the suggested strategy is to form a

library of for single surfaces S, dependent on shape parameters, and for

material parameters characterizing common acoustic media. The desired field

for a particular S and k would then be given by (9). Of course it would take

time and effort using existing BEM to form A'^\ However, we believe that ease and speed of subsequently obtaining scattered fields, using A^' from the

library, would more than justify the creation of individual entries. This has

certainly been the case in our experience with the small library we have

109
already created for our own use. Like any library, the larger the better, but

creation costs would suggest some advance determination of parameters

characterizing those A^^ which would be used repeatedly. The number of

entries for acoustic scattering would probably be more a function of the

variety of shapes rather than the number of different acoustic media. The size

of the individual entries would be a function of the degrees of freedom, or number of nodes N in the BEM used to get the In any case, all of the work

to get A^^ would be done in advance of their use by the library formers.

For the users, the time to get u(P) via (9) would depend primarily on N.

The user specifies /, and picks the particular A^^ needed, then /(/,) would be

formed by library software. Next, the software would multiply f(l^) by

A^\l^,qj) to get a row dependent on qj. Subsequent operations by the same

software to get u would involve "(IxN) row times (Nxl) column" operations for

each (perhaps many) P from a list of P specified by the user. A new

multiplication involving the (NxN) A^\/.,^p would be required only for new /.

For the two-surface problem the situation is a bit more complicated

and more computationally intensive. Here, the strategy is to set up and solve

the BIE equivalent of (4) or (5) wherein only Si needs be discretized. For this, as already noted, G* replaces G. In principle then, solving the two-surface

problem using G* is formally identical to solving the one-surface problem via

the conventional BEM using G. However, to actually implement this formality,

both G*(/7j,(7j) and its normal derivative at qi are required, where neither pi

nor qi are on S2. Then, to finally get u(P) via (13), one needs G*(P,q^) and its normal derivative at qi. Again, neither point is on S2.

Getting the values of G and its normal derivative at the mentioned points

is easy, but getting the corresponding values of G* and its normal derivative at

the same points requires (17) and (18). (Expression (7) is insufficient since

q=qi is required). Thus, just to obtain values of the required Green's functions

at points off of S2, more matrix multiplications via (17) and (18) are needed

than to get u(P) via (9), for the one-surface problem. Then, after obtaining

those values, there still remains the job of using them to set up the BIE (5),

solving (5), and finally using (13) to get u(P) for the two-surface problem.

Nevertheless, despite the number of steps, there is considerable

conceptual and strategic advantage in utilizing a library-collection of A^^ for

the two-surface problem, as well as for the one-surface problem. Matrix

multiplications are the main computational burden in (17) and (18). However,

such are less formidable than discretizing a commonly-occuring but complex

shape S2 each time - not to mention that the operator (matrix) A* which needs

to be inverted is smaller than the A (for the union of Si and S2) when only G is used.

With the library in place, this, in a nutshell, is the tradeoff: less user

expertise and modeling effort plus certain CPU time required to solve a problem using G*, versus less of the first two items, and probably less CPU too,

than when using G. We believe in the long run this is a good trade.

110
In any case, having a library of G*, in analytical or discretized form,

the BEM should become a more powerful, more accurate and even a faster tool,

which would be usable on smaller computers, by less-expert users.

Acknowledgement

Author PAM is grateful to Iowa State University, through the Glenn Murphy

Professorship in Engineering, for funding a summer 1994 visit to Ames, lA,

during which some of the research for this paper was done. FJR and RAR are

grateful to David Budreck for valuable conversations on discretized Green's

functions. FJR wishes to acknowledge several associates and students,

especially colleagues Tom Rudolphi and Tom Rogge, for some clarifying

discussions on the structure of the equations in this paper, and Abhijit

Chandra [ 7] for some ideas on Green's functions for modeUng of materials

behavior.

References

1. V. Tewary, Elastic Green's Functions for Anisotropic Solids, Proc. NIST

Workshop on Green 's Functions, BEM, and Advanced Materials

Modeling, Boulder, Colorado, August, 1994.

2. P.A. Martin and F.J. Rizzo, Partitioning and Boundary Integral

Equations, (preprint).

3. D.E. Budreck, A Discretized Green's Function for an Elastodynamic Crack,

{preprin t).

4. R.A. Roberts, The Asymptotic Computational Ansatz: Application to Critical

Angle Beam Transmission Boundary Integral Equation Solution, J. Acoust.

Soc. Am. 95 (4), 1994.

5. B.A. Boley, A Method for the Construction of Green's Functions, Q, Appl.

Math. 14(3), 1956.

6. S. Bergman and M. Schiffer, Kernel Functions and Elliptic Differential

Equations in Mathematical Physics, Academic Press, New York, 1953.

7. A. Chandra, Y. Huang, X. Wei, and K. X. Hu, A Hybrid Micro-Macro

BEM Formulation for Micro-Crack Clusters in Elastic Components, Int J. Num.

Meth. Eng. ( to appear).

Ill

Fracture and Crack Growth of Polymeric Composites for Use in Dynamic

Applications

Ronald A. L. Rorrer

Senior Materials Project Engineer

Advanced Materials Research - Power Transmission Division

Will Fehringer

Research Engineer

Advanced Technology & Engineering Services

The Gates Rubber Company

Denver, Colorado

Abstract

Examination of the technical concerns and issues in failure of polymeric, specifically elastomeric,

composites used in dynamic applications. One of the dominant failure mechanisms is related to the initiation and propagation of cracks in the composite, which often lead to subsequent

catastrophic failure. With multi-component composites, the design issues center around material

selection, engineering of the component materials and interfaces, and geometry changes.

Engineering of the component materials and interfaces modifies the composite strength, while

geometry changes modify the applied loading (stress or strain). The overall goal is to extend the

useable life of the composite by elimination or reduction of crack propagation as a dominant

failure mode. Examples of typical failure paths will be presented and discussed for an

elastomeric composite.

Introduction

While there are certainly many uses of polymeric composites in static structural applications, a

great deal of polymers and elastomers are used in dynamic load and strain applications due to their light weight and flexibility. Due to non-linear viscoelastic nature of these materials, the

resulting responses to the dynamic loading become very complex. For example, even if the

materials are considered linear, the viscoelastic moduli, or specifically the losses, result in

internal heat generation which can impact the final response. Additionally, the material

properties have a time and temperature dependence (1), that is usually not a consideration with

metals.

Polymeric Composites

Advanced polymeric composites can be as simple as single component polymer matrices with

reinforcement fibers. Or they can be as complex as an elastomeric composite with. multiple

113
layers of different elastomers, tensile members made of steel and polymers (such as aramid,

polyester), and other polymeric components such as fabrics. Typical examples of the latter are

belts, tires and hoses. For example, a belt may contain up to 5 different layers of elastomer.

Many of the failure modes of these dynamic composites are from propagating cracks that reach a point where the operation of the product is compromised or catastrophic failure ensues. When

selecting materials for a multi-component composite, it is desirable to optimize the properties of

the various components relative to one another. Generally, any enhancement of material

properties incurs a cost penalty. Therefore, it is rare in elastomeric composites that there will

be orders of magnitude difference in material properties between layers of similar materials. The

exception to this generality is the specific case of tensile members, which will possess an order

of magnitude difference in modulus and strength.

The interfaces between different materials are another concern of a polymeric composite. The

interface should posses sufficient strength or adhesion to avoid failure at stress or strain levels

much below the properties of the adherends. Additionally, the interfacial bond between fibers

and matrix are often neglected (2) in fiber reinforced composites. This becomes important, since

for many polymers and virtually all elastomers, fillers are used for reinforcement.

Where does this lead? Optimizing a construction typically means ensuring that a "weak" material or interface can meet the demands of the rest of the composite. In fact, the ideal composite may

be one, where all materials and interfaces possess the same factor of safety or lifetime. This

situation yields the desired life as well as minimal construction costs.

Due to the time and financial cost of testing final products, it is desired that the individual

materials and their composite constructions response be predictable. Given that each material

and interface in a composite possesses characteristic inherent flaws and propensities of crack

propagation, modeling of the composite should identify a priori the location of potential problems

due to the distribution and magnitude of stress and strain.

Theoretical Considerations

Fracture properties are usually determined from either the "pure shear" test or the tensile test as shown in Fig. 1. The nomenclature of the testing mode refers to the strain state of energy

storage in the sample. Adhesion between layers of dissimilar materials is also determined by a modifed tensile test or 180 ° peel test as shown in Fig. 1(c). The deformation energy partitions

between stored energy that is available for recovery and fracture events and the energy dissipated

as heat (3). Dissipation of tear energy in polymeric materials result in an overestimate of the

energy available for the fracture or crack propagation event as predicted from stress-strain

data (4).

Figure 1. Fracture/Fatigue Testing Geometries.

114
Simple analytical modeling is applicable to simple geometries and single material system, such

as the test specimens shown in Fig. 1. However, when there is more than one material in a system and the geometry is no longer simple, FEA techniques must be relied on to predict the

final behavior. For example, in Linear Elastic Fracture Mechanics (LEFM) the total stored

strain energy can be estimated for a multi-component composite for a given loading and also the

direction and propensity of crack propagation with FEA. In elastomeric fracture there has been no universally accepted failure theory, such as Tresca or von Mises under combined states of

stress (5). Measurements (4) have shown that fracture energies in Mode III (shear) are 50%

higher than Mode I (tensile). Historically, this difference has been attributed to friction in the

shear failure mode. However, one view is that rubber fails not along shear planes, but normal

to the principal stresses (5). An alternative view of this, is that it is virtually impossible to maintain a state of shear at the crack tip due to the high extensibility.

The general case of of crack propagation with 3-dimensional loading and anisotropic material

properties is of interest in analyzing polymeric composites. While there are many cases of 1 or 2 dimensional loading and orthotropic material properties, the more complex state of stress can be present. Crack initiation and initial propagation directions occur in areas of obvious stress concentration. However, the propagation path is not always obvious, due to the different layers

of material and interfaces and their relative properties.

Example Cases

Examples of elastomeric composites experiencing crack growth and subsequent failure under

dynamic loading are V and V-ribbed belts. A cross-sectional view of the belts are shown in

Figs. 2(a) and 2(b) respectively. As in any composite, there are n multiple phases or layers of

material that in this specific case result in n-1 interfaces. In the V-belt case there can be 7 or more individual constituents with 6 interfaces. While the V-Ribbed Belt can have 5 or more

individual constituents with 4 interfaces. Obviously, the number of layers are not fixed, but vary

with application demands.

Figure 2. Belt Configurations.

While quite commonplace the loading and operating demands on a V-belt can be quite complex.

Two general configurations of V-belts are shown in Fig. 3(a) and 3(b). The first configuration

3(a) is known as fixed center distance on the drive setup. The second configuration 3(b) is known as a tensioned drive due to the backside idler. Both systems have many similarities in the

loading. The tension of the drive varies around the pulleys and thus the drive system. The

variation of tension puts a cyclic tension on the elements in the belt. The wedging of the V-belt

in the pulley creates a three dimensional state of stress and strain while the belt is in the pulley.

Additionally, the belts experience flexural stress and strain, when wrapped around the pulleys.

The major difference from a stress or strain point of view between fixed center and tensioned

drives is the presence of the backside idler which creates a reversal of stress and strain on the

undercord.

115
Figure 3. General Drive Configurations of a V-belt.

The following figures (Figs. 4,5, and 6) show three separate crack locations on a single V-belt

in the undercord region. All three cracks are initiated at the apex of a the cog. Initial evidence

of cracking is usually not considered failure of the belt. However, as the crack progresses

undercord sections will eventually "chunk out". Figure 4 shows a single crack starting at the

apex whcih eventually truns as it approaches another material layer. There is a secondary crack

which has started off of the primary crack and is progressing to the left of the photomicrograph.

Figures 5 and 6 show double apex cracks. The primary cracks progress in the direction of the

single apex crack shown in Fig. 4.. In both cases there is the presence of a secondary apex

crack. The path and propagation of the secondary crack can lead to "chunk out", which is on

the verge of occurring in Fig. 5.

Figure 4. Single Apex Crack in Cog of a V-Belt

Figure 5. Double Apex Crack with Incipient Material Loss in Cog of a V-Belt

Figure 6. Double Apex Crack in Cog of a V-Belt

Figure 7 shows a sidewall crack in a V-belt which occurs in a transverse and longitudinal

direction to the belt, as compared to the apex cracks which originate in a radial (or through the

thickness direction).

Figure 7. Sidewall Crack in a V-Belt

V-Ribbed Belt

Figure 8 shows an example of a flex crack in a V-Ribbed belt. The flex cracks occur

preferentially in drives with backside idlers due to the tensile stress/strain imposed by the

backside idler.

Figure 8. Flex Crack in Undercord of V-Ribbed Belt.

Summary

In summary, there are many design parameters that can be varied to minimize the crack initiation

and propagation in an elastomeric composite. Reiterating, the parameters are material storage

and loss moduli, the interfaces between dissimilar materials, and overall construction. These

parameters in turn determine the level of stress and strain as well as the fundamental fracture and

fatigue properties of the composite. As can be seen from the apex cracks in the V-belt examples,

direction of the crack will determine the failure mode.

This raises the following question. What is needed from fracture mechanics, fatigue and

116
modeling? Given the fundamental material properties (storage and loss moduli, strain energy

release rate) and geometry, the propensity of a composite to propagate a crack should be

predictable. This predictive ability would enable the product designer to evaluate the relative

impact of changes without having to build prototypes of each proposed modification. This is especially significant when the product development cycle is on the order of 1/2 to 2 years.

References

(1) Abou-Hamda, M., Mai, Y., and Wu, S., "Analysis of Fatigue Crack Growth in A

Rubber-Toughened Epoxy Resin: Effect of Temperature and Stress Ratio," Polymer, 34,

pp. 4221-4229, (1993).

(2) Latour, R.A.J., Black, J., and Miller, B., "Fracture Mechanisms of the Fiber/Matrix

Interfacial Bond in Fiber-reinforced Polymer Composites," Surface and Interface

Analysis, 17, pp. 477-484, (1991).

(3) Borowczak, M., and Causa, A.G., "Fatigue Behavior of Cord-Reinforced Rubber

Composites," Polymers & Polymer Composites, 1, pp. 211A-217A, (1993).

(4) Ahagon, A., Gent, A.N., Kim, H.J., and Kumagai, Y., "Fracture Energy of Elastomers

in Mode I (Cleavage) and Mode III (Lateral Shear)," Rubber Chemistry and Technology,

48, pp. 896-901, (1975).

(5) Kinloch, A.J., and Young, R.J.

; "Fracture Behavior of Polymers," Elsevier Science

Publishers LTD, Essex England, 1990, pp. 1-496.

117
(a) Pure Shear (b) Tensile

Figure 1. Fracture/Fatigue Testing Geometries.

(c) Tensile or Peel

(a) Fixed Center Drive (b) Tensioned Drive

Figure 3. General Drive Configurations of a V-belt.

118
Figure 4. Single Apex Crack in Cog of a V-Belt
120
Figure 8. Flex Crack in Undercord of V-Ribbed Belt.

121

Elastic Green's Functions for Anisotropic Solids

V.K. Tewary

Materials Reliability Division

National Institute of Standards and Technology

Boulder, Colorado 80303

1. Introduction

We briefly review the Green's function method for solution of the Christofifel equation and

a computationally efficient method for calculating the Green's function for anisotropic

solids. We describe the virtual-force method for satisfying the boundary conditions to

account for discontinuities in the solids. The virtual-force method consists of applying a

distribution of virtual forces just outside the domain of solution. The solution obtained by

using the virtual-force distribution and the Green's function gives a solution of the

homogeneous equation. The virtual-force distribution is then determined by imposing the

prescribed boundary conditions. This method is similar to the image-charge method in electrostatics [1] and is the basis for the boundary-element method [2] for solving

elastodynamic problems.

The Fourier representation of the Green's function is quite general and, subject to certain

well-known conditions of integrability and convergence [1], can be used for most physical

problems. In the case of elastodynamic Green's functions, the Fourier representation is CPU intensive and is not computationally efficient for anisotropic solids. We have

developed a delta-function representation [3] that is particularly suitable for anisotropic

solids. In this paper, we describe the deha-fianction representation for elastodynamic as well as elastostatic Green's functions for infinite solids and its application to bounded

solids by using the virtual- force method.

2. Green's function method

We represent the space and time variables by x and t, respectively. The Cartesian

components of a vector will be denoted by indices i, j, k, etc., which assume the values

1,2, or 3, corresponding to x, y, and z coordinates. Summation over repeated Roman

indices will be implied unless stated otherwise.

The Christoffel equation for elastic equilibrium is written in the operator form as

Ly Uj (x,t) = Fi (x,t). (1)

where

Lij =Cikji5Vaxkaxi - SijpaVat^ (2)

123
c is the fourth-rank elastic-constant tensor, p is the density of the solid, and F(x,t) is the

applied force. The Green's function G(x,x'; t,t') is the solution of the equation

*: LijGjk(x,x';t,t')= 5u,5(x-x')5(t-t'X (3)

where x' and t' are variables in the same space as x and t, respectively; 5ik is the

Kronecker's delta tensor that is 1 for i=j, and 0 otherwise; and 5 (x) is the Dirac delta

fiinction defined by the relation

w(a) = j 5 (p-a) w(p) dp. (4)

In eq (4), a is a constant on the real axis and w(p) is any arbitrary integrable function of p, and the integration is over the entire real axis. We have used the same symbol for the

Dirac delta function and the Kronecker delta tensor since they can be identified by their

arguments.

The particular solution of eq (1), which can be verified by applying the operator L and

using eqs (3) and (4), is given by

Up(x,t) = J G(x,x';t,t') F(x',t') dx' dt'. (5)

The integration in eq (5) is over the entire space of x and t. For an infinite solid with no

boundary conditions prescribed over space, eq (5) gives the final solution. For infinite

solids, G depends on x and x' and t and t' only through their differences x-x' and t-t', respectively. In such cases G can be denoted by a single space and time variable as G(x-x', t-t') or G(x,t) that implies x'=t'=0.

For solids with spatial discontinuities, such as bounded solids or those containing holes,

crack, interfaces, etc., we need to satisfy some prescribed boundary conditions. Let S

specify the space of the solid over which we need to solve eq (1) with boundary conditions

prescribed over the surface of S. Equation (5) still gives the particular solution of eq (1)

with the integration over x' restricted to S. To obtain the homogeneous solution, we apply

a distribution of virtual forces f(xs,t) over the surface of S. The function f(x,t) is 0

everywhere except for x = xs, where xs lies over the surface of S. The homogeneous

solution, as in eq (5), is given by

UH(x,t) = I G(x-x'; t-t') f(x',t') 6 (x'-xg) dx' dt'. (6)

Since eq (6) gives the homogeneous solution for arbitrary f(xs ,t), we determine this

function by imposing the boundary conditions. This is the essence of the Green's function

method. For solids with simple geometrical surfaces, we can determine f(xs ,t) analytically.

For solids with complicated geometrical shapes, f(xs,t) has to be determined numerically

as is done in the boundary-element method. The solution of eq (1), that is the

124
displacement field u(x,t) for an applied unit force, subject to all the prescribed boundary

conditions, is the Green's function for the solid.

The above derivation shows how Green's function for any geometrical structure can be

built up in stages or modules. For example if G(x,t) is the free space Green's function for an infinite solid, then to obtain Green's fianction for a semi-infinite solid with one free

surface, we apply a virtual force just outside the free surface and add a homogeneous

solution to the solution for the infinite solid. We then determine the virtual force by

applying the boundary condition at the free surface. The Green's flinction thus obtained

will be the Green's function for the semi-infinite solid. If we want to add another free

surface or any other discontinuity like a hole in the solid, we apply another virtual force at the new surface or the discontinuity. We add this homogeneous solution to that obtained

by using the Green's function for the semi-infinite sold. We determine the new virtual

force by imposing the additional boundary conditions at the discontinuity. Only the

additional boundary condition needs to be satisfied since the semi-infinite Green's function

will automatically satisfy the boundary condition at the first free surface.

In the linear case, f(xs, t) will be proportional to F(x,t). Hence the boundary conditions

will be satisfied for all F(x,t). The Green's function thus obtained will therefore be

independent of F(x,t). The Green's function is a characteristic of the solid including

discontinuities, if any, and does not need to be recalculated for a different applied load. It should be usefial, therefore, that the Green's functions for typical geometrical shapes can be stored in a central computer and made available to other users. The idea of such a

library of Green's fiinctions has been recently suggested by Rizzo [4].

3. Integral representation for Green's functions

In general it is possible to use the three-dimensional (3D) Fourier integral representation

for the Green's function as given below:

G(x,t) = (27i)-^ I Gf(]K, CO) exp i(K.x - cot) dK dco (7)

where = -1, Gf(K, cd) is the Fourier transform of the Green's fianction, K is the wave

vector, and co is the frequency; and the integration in eq (7) is over the entire space. From

eqs (2) and (3)

Gf(K, q) = [A (K) - 1 pco^]-' , (8)

where

A ij (K) = Cikj, Kk K| . (9)

The 3x3 matrix A(K) is the Christoffel matrix in Fourier space. It is the long-wavelength

(low-K) limit of the Born-von Karman dynamical matrix [5]. Its eigenvalues co^(K) are the

125
K

squares of the phonon frequencies and its eigenvectors are the polarization vectors of the

corresponding phonons. Equation (7) along with eqs (8) and (9) can be used for

calculating the Green's function. The Green's function thus calculated will not be causal.

To ensure causality, we can introduce a small imaginary part in co [5] and take the limit as the imaginary part approaches 0. Alternatively, we can take the Laplace transform over

time and use the Laplace inversion integral in eq (7).

For an isotropic solid, the matrix inversion in eq (8) and the integration in eq (7) can be

done analytically. The Green's fianction can also be obtained analytically [6] from eq (7)

for a line force, or a 2D approximation, since the component of K in the direction of the

line force is 0. For a general 3D anisotropic solid, eq (7) requires a 4D numerical

integration- over three components of K, and one frequency variable. In general Gt{K, co) has singularities (resonances) on the real axis. The integral in eq (7) is defined in the

Cauchy sense. It involves evaluation of principal values that creates problems of numerical

convergence. A numerical evaluation of the general 3D anisotropic Green's function using

the Fourier representation is CPU intensive. Some shortcuts for evaluating the Green's

function are available in the literature (see refs. [3] and [7] and other references quoted

therein).

We have developed a delta-function representation of the Green's function [3] that is computationally convenient even for 3D anisotropic solids. In this representation we write

G(x,t) = (471^)-' j Gq(q) 5<"(t-q.x ) H(t) dq, (10)

where

Gq(q) = Lim,^.oIm[A(q)-(l-i8)I]^ (11)

A ij (q) = c,kji qk qi

,

H(t) is the Heaviside step function, being unity for t >0 and 0 for t<0 , 6^^^ is the first derivative of the delta function with respect to its argument, and q, that has the

dimensions of inverse velocity, is a vector in slowness space . We identify A(q) as the

Christoffel matrix and Gq(q) as the Green's function in the slowness space. The delta

function in eq (10) is a statement of the physical fact that a phonon of slowness value q

will reach the distance x in time t, then its velocity 1/q must be x/t. The function Gq (q)

represents physically the weight or number of phonons of slowness vector q that the solid can provide. If, instead of the imaginary part on the RHS of eq (1 1), we take the real part, use the delta function of t-q.x instead of its derivative, and remove H(t), then eq (10) is the Radon transform of the elastodynamic Green's function. The Radon representation of

the Green's function has been developed by Wang and Achenbach [7] and applied to

several interesting cases.

(12)

126
Using the representation of the delta function, we obtain from eq (1 1)

G, (q) = -% Is Cs, (q) Csj (q) 6 [E\ (q) - 1], (13)

where Cs (q) (s =1,2,3) and E\(q) are, respectively, the eigenvectors and eigenvalues of

A(q). The right hand side of eq (10) requires integration over three variables - the three

components of q. However, the integrand is a product of two delta functions. Hence the

integration over any two of the three variables is done analytically simply by substituting

for their values determined by the delta functions. Numerical integration is required over

only one variable. Moreover, the integrand is not singular and does not contain oscillatory

functions. Consequently the delta-function representation is computationally much more

efficient relative to the Fourier representation in eq (7).

Finally, the solution of the homogeneous equation, which can be verified by direct

substitution for arbitrary f(q), is given by

UH(x,t) = (471^)-' I G,(q) f(q) 5('>(t-q.x) dq. (14)

The virtual force f(q) has to be determined by imposing the boundary conditions.

Application of eq (14) to calculate the elastic waveforms in anisotropic solids has been

given in [3]. Equations (10) and (14) reduce to corresponding elastostatic Green's

function in the limit t=+0.

References

[1] Morse P.M. and Feshbach H., Methods of mathematical physics. New York: McGrawHill; 1953, Part 1.

[2] Berger J.R., Boundary element analysis of anisotropic bimaterials with special Green's

functions. Eng. Anal. Boundary Elem. 14: 123; 1994.

[3] Tewary V.K., Computationally efficient representation for elastodynamic and

elastostatic Green's functions for anisotropic solids. Phys. Rev. B 121: 15695; 1995.

[4] Rizzo F.J., Boundary element modeling for nondestructive evaluation using uhrasonic

and elctromagnetic fields and the idea of discretized Green's functions. These procedings.

[5] Maradudin A.A., Montroll E., Weiss G., and Ipatova I., The theory of lattice

dynamics in the harmonic approximation. Solid State Physics, Supp. 3: New York:

Academic Press; 1971.

[6] Tewary V.K., Wagoner R.H., and Hirth J.P., Elastic Green's function for a composite

solid with a planar interface. J. Mater. Res. 4: 113; 1989.

[7] Wang C.Y. and Achenbach J.D., Elastodynamic fijndamental solutions for anisotropic

solids. Geophys. J. Int. 118: 384; 1994.

127

Appendix A. Organizing Committee

Bradley K. Alpert

Computing and Applied Mathematics Laboratory

National Institute of Standards and Technology

Boulder, Colorado

John R. Berger, Conference Secretary

Division of Engineering

Colorado School of Mines

Golden, Colorado

Thomas L. Geers

Department of Mechanical Engineering

University of Colorado

Boulder, Colorado

Sharon Glotzer, ex officio member

Polymer Division

National Institute of Standards and Technology

Gaithersburg, Maryland

Ronald A. L. Rorrer

Advanced Materials Research, Power Transmission Division

The Gates Rubber Company

Denver, Colorado

Vinod K. Tewary, Chair

Materials Reliability Division

National Institute of Standards and Technology

Boulder, Colorado

129

Appendix B. Program

NIST Workshop on

Green's Functions and Boundary Element Analysis

for Modeling of Mechanical Behavior of Advanced Materials

Sunday, August 14, 1994

7:00-9:00 p.m. Registration and reception for workshop attendees

Monday, August 15, 1994

7:30-8:00 a.m.

8:00-8:15 a.m.

Session 1

Chair

8:15-8:45 a.m.

8:45-9:15 a.m.

9:15-10:00 a.m.

10:00-10:15 a.m.

Session 2

Chair

10:15-11:00 a.m.

10:00-11:45 a.m.

Registration

V. K. Tewary, NIST

Welcome and Opening of the Workshop

R. J. Rizzo, Iowa State University

H. I. McHenry, NIST

Overview of NIST

D. Hall, NIST

The Center for Theoretical and Computational Materials Science

at NIST

B. Eichinger, Biosym Technologies

Materials Research and Software Development at Biosym

Break

B. Eichinger, Biosym Technologies

G. Krishnasamy, PDA Engineering

Bouldary Element Method for Computer Aided Design

L. C. Davis, Ford Research Laboratory

Mechanical Properties of Metal Matrix Composites

11:45 a.m.-l:15 p.m. Lunch Break

131
Monday, August 15, 1994, continued

Session 3

Chair

1:15-2:00 p.m.

2:00-2:45 p.m.

S. Mukherjee, Cornell University

T. A. Cruse, Vanderbilt University

Green 's Functions in Elastic Fracture Mechanics Analysis

B. S. Annigeri, United Technologies Research Center

Crack Growth Modeling Using the Surface Integral and Finite

Element Hybrid Method

2:45-3:00 p.m. Break

Session 4

Chair

3:00-4:30 p.m.

B. K. Alpert, NIST

Poster Session

M. Denda, Rutgers University

Green 's Function BEM Research at Rutgers CMAS and its Applications to the Modeling of Mechanical and Piezoelectrical

Behaviors of Advanced Materials

M. Dunn, University of Colorado

Electroelastic Green 's Functions and their Use in Piezoelectric

Material Modeling

L. Gray, Oak Ridge National Laboratory

Stress in Epitaxial Thin Films

Z. Q. Ziang and A. Chandra, University of Arizona

A Hybrid Micro-Macro BEMfor Ceramic Composites

I. Kaljevic, Ohio Aerospace Institute Stochastic Boundary Elements for Two-Dimensional Potential Flow

in Nonhomogeneous Media

S. Mukherjee, Cornell University

Overview of BEM Research in Mukherjee 's Group at Cornell

E. Pan, University of Colorado

3-D Boundary Element Formulation of Anisotropic Elasticity

with Gravity

S. Saigal, Carnegie Mellon University

Inverse Elasticity Solutions Using Boundary Elements

132
Tuesday, August 16, 1994

Chair

8:15-8:45 a.m.

8:45-9:30 a.m.

9:30-9:45 a.m.

Session 6

Chair

9:45-10:30 a.m.

10:30-11:00 a.m.

11:00-11:30 a.m.

11:30 a.m. -1:00 p.m.

Session 7

Chair

1:00-1:45 p.m.

1:45-2:15 p.m.

2:15-2:45 p.m.

2:45 p.m. -3:00 p.m.

R. Thomson, NIST

R. A. L. Rorrer, Gates Rubber Company

Fracture and Crack Growth of Polymeric Composites for Use

in Dynamic Applications

A. Neid, General Electric Company

Simulation of the Electrochemical Machining Process

Using a 2D Fundamental Singular Solution

Break

T. A. Cruse, Vanderbilt University

R. Thomson, NIST

Lattice Green 's Functions for the Study of Defects in Solids

V. K. Tewary, NIST

Elastic Green '5 Functions for Anisotropic Solids

J. R. Berger, Colorado School of Mines

Boundary Element Analysis of Bimaterials Using Anisotropic Elastic

Green's Functions

Lunch Break

C. Fortunko, NIST

F. J. Rizzo, Iowa State University

Boundary Element Modeling for Nondestructive Evaluation Using

Ultrasonic and Electromagnetic Fields and the Idea of Discretized

Green 's Functions

T. L. Geers, University of Colorado

Approximate Operators for Dynamic Boundary Element Analysis

B. K. Alpert, NIST

Accurate Discretization of Integral Operators

Break

133
Tuesday, August 16, 1994, continued

Session 8

Chair

3:00-4:30 p.m.

S. Glotzer, NIST

Group Planning Session

Fracture (Defects, Interfaces, Composites)

Quantitative Nondestructive Evaluation (Ultrasonics, Inverse Problems)

Computer Aided Engineering (Stress Analysis, Design,

Optimization)

Wednesday, August 17, 1994

Session 9

Chair:

8:15-9:45 a.m.

9:45-10:00 a.m.

10:00-11:30 a.m.

11:30 a.m.

T. Geers, University of Colorado

Presentation and discussion of summaries from the

group planning sessions.

Break

Continuation of sessions

Workshop adjourns

134
Appendix C. Participants

NIST Workshop on

Green's Functions and Boundary Element Analysis

for Modeling of Mechanical Behavior of Advanced Materials

Bradley K. Alpert

Computing and Applied Mathematics

Laboratory

National Institute of Standards and

Technology

325 Broadway

Boulder, CO 80303

(303) 497-5920

alpert@boulder.nist.gov

Balkrishna S. Annigeri

Mail Stop 73

United Technology Research Center

41 1 Silver Lane

E. Hartford, CT 06108

(203) 727-7129

ba@utrc.utc.com

John R. Berger

Division of Engineering

Colorado School of Mines

Golden, CO 80401

(303) 273-3682

jberger@mines.colorado.edu

Abhijit Chandra

AME Department

University of Arizona

Tucson, AZ 85721

(602) 621-2080

chandra@idea.ame.arizona.edu

Li Chang

NASA Lewis Research Center

21000 Brook Park Road

Cleveland, OH 44135

(216) 433-8119

Thomas A. Cruse

Department of Mechanical Engineering

Vanderbilt University

Box 1592, Station B

Nashville, TN 37235

(615) 343-8727

tcruse @vuse.vanderbilt.edu

L. Craig Davis

Physics Department

Ford Research Laboratory

Mail Drop 3028

P. O. Box 2053

Dearborn, MI 48121

(313)322-7006

Idavis @smail.srl.ford.com

Mitsunori Denda

College of Engineering

Rutgers Unviersity

P.O. Box 909

Piscataway, NJ 08855

(908) 932-4391

denda @j ove .rutgers.edu

Marty Dunn

Department of Mechanical Engineering

University of Colorado

Boulder, CO 80309

(303) 492-6542

dunnm@spot.colorado.edu

Vincente Dyreyes

Grumman Aerospace Corporation

B43-035 Stewart Avenue

Bethpage, NY 11714

(516) 346-8475

135
Bruce Eichinger

Biosym Technologies

9685 Stranton Road

San Diego, CA 92121

(619) 546-5540

bruce@biosym.com

Chris Fortunko

Materials Reliability Division

National Institute of Standards and

Technology

Boulder, CO 80303

(303) 497-3062

fortunko@boulder.nist.gov

Thomas L. Geers

Department of Mechanical Engineering

University of Colorado

Boulder, CO 80308

(303) 492-6355

geers @spot.colorado.edu

James Given

National Institute of Standards and

Technology

Gaithersburg, MD 20899

(301) 975-2447

given@tiber.nist.gov

Sharon Glotzer

Polymer Division

National Institute of Standards and

Technology

Gaithersburg, MD 20899

(301) 975-6775

glotzer@polygon.nist.gov

Leonard Gray

Oak Ridge National Laboratory

P.O. Box 2008

Building 6012

Oak Ridge, TN 37830

(615) 574-8189

ljg@ornl.gov

Dale Hall

Office of Intelligent Processing of

Materials

National Institute of Standards and

Technology

Gaithersburg, MD 20899

(301) 975-5727

hall@micf.nist.gov

Igor Kaljevic

Ohio Aerospace Institute

22800 Cedar Point Road

Brook Park, OH 44142

(216) 962-3000

Guna Krishnasamy

PDA Engineering

2975 Red Hill Ave.

Costa Mesa, CA 92626

(714) 540-8900

krishna@pda.com

Harry I. McHenry

Materials Reliability Division

National Institute of Standards and

Technology

Boulder, CO 80303

(303) 497-3268

mchenry@micf.nist.gov

Subrata Mukherjee

Theoretical and Applied Mechanics

Kimball Hall

Cornell University

Ithaca, NY 14853

(607) 255-7143

subrata@khan.tam.cornell.edu

Graham Mustoe

Division of Engineering

Colorado School of Mines

Golden, CO 80401

(303) 273-3661

gmustoe@mines.colorado.edu

136
Arthur Neid

Process Physics Laboratory

General Electric Company

Building KW, Room D276

P.O. Box 8

Schenectady, NY 12301

(518)387-5529

Emian Pan

Department of Civil Engineering

University of Colorado

Boulder, CO 80309

(303) 492-7636

pane@spot.colorado.edu

Frank J. Rizzo

Department of Aerospace Engineering

and Engineering Mechanics

2019 Black Engineering Building

Iowa State University

Ames,IA 50011

(515) 294-7399

frizzo@cnde.iastate.edu

Ronald A. L. Rorrer

Advanced Materials Research

Gates Rubber Co.

P. O. Box 5887

Denver, CO 80217

(303) 744-4289

rrorrer@mines.colorado.edu

Sunil Saigal

Department of Civil and Environmental

Engineering

Carnegie Mellon University

Pittsburg, PA 15213

(412) 268-2077

ss6o@andrew.cmu.edu

Vinod K. Tewary

Materials Rehability Division

National Institute of Standards and

Technology

325 Broadway

Boulder, CO 80303

(303) 497-5753

tewary@boulder.nist.gov

Robb Thomson

Materials Science and Engineering

Laboratory

Building 223

National Institute of Standards and

Technology

Gaithersburg, MD 20899

(301)975-5665

robb@phlogiston.nist.gov

137

Appendix D. Abstracts

139

Accurate Discretization of Integral Operators

B. K. Alpert

Computing and Applied Mathematics Laboratory

National Institute of Standards and Technology

Boulder, Colorado

The numerical solution of integral equations involves reduction of equations on a function space to a finite system of equations, typically over a chosen set of basis

functions (the Galerkin method). The accuracy of the computed solution is determined by

the accuracy with which the true solution can be represented by the finite basis, as well as the accuracy of projections onto the basis. An alternative discretization technique, the

Nystrom method, which reduces the integral equations to a finite system of equations by

replacing each integral with a quadrature, can often be designed to obtain better

convergence. We discuss some quadrature techniques for handling kernel and solution

singularities arising in integral equations and give numerical examples of their accuracy.

141
Efficient Modeling of Two and Three Dimensional Crack Growth

Using the Surface Integral and Finite Element Hybrid Method

B. S. Annigeri

United Technologies Research Center

East Hartford, Connecticut

WilHamD. Keat

. ;^ / Clarkson University

, , . Potsdam, New York

The need for modeling crack growth in materials is important for assessing the

fatigue life of structural components. In a variety of industries: aerospace, nuclear,

automobile and others; fatigue crack growth can be a significant source of failure of

structural components. This paper describes recent advances in modeling crack growth

using the Surface-Integral and Finite Element Hybrid method. Cracks in an infinite or

semi-infinite domain are modeled using a continuous distribution of dislocations in 2D or

force multipoles in 3D resulting in a singular integral formulation. The uncracked body is modeled using finite elements. The coupling of the surface-integral and finite element

models is obtained by traction and displacement matching on the external and internal

boundaries of the finite plate with crack(s). A significant advantage is that the finite

element mesh remains fixed as the crack propagates. The remeshing then is required only on the surface integral discretization which is accomplished automatically without user

intervention. This paper describes the recent effort of modeling through cracks in layered

materials and modeling of embedded and surface cracks in 3D. Numerical predictions

have shown good correlation with experimental data and other analytical solutions. The

research efforts have resulted in the development of the SAFE-2D and SAFE-3D codes

for effective modeling of crack growth.

142
Boundary Element Analysis

of Bimaterials Using Anisotropic Elastic Green's Functions

J. R. Berger

Division of Engineering

Colorado School of Mines

Golden, Colorado

The boundary integral equations incorporating the Green's function for anisotropic

solids containing planar interfaces are presented. The fundamental displacement and

traction solutions are determined from the displacement Green's function of Tewary,

Wagoner, and Hirth [Journal of Materials Research, Vol. 4, p. 113-123]. The

fundamental solutions are shown to numerically degenerate to the Kelvin solution in the

homogeneous, isotropic limit. The boundary integral equations are formulated with the use of constant boundary elements. The constant boundary elements allow for analytic

evaluation of the boundary integrals. The application of the method is demonstrated by

analyzing a copper-solder system subjected to mechanical loading.

143
Green's Functions in Elastic Fracture Mechanics Analysis

T. A. Cruse

Dept. of Mechanical Engineering

Vanderbilt University

Nashville, TN 37253

The focus of my discussion of Green's functions will be on the treatment of the

special boundary conditions associated with elastic fracture mechanics analysis, rather

than on developing new Green's functions for generalized differential operators. The talk

will briefly review the difficulty of BIE analysis for elastic fracture mechanics problems

which derives from the degeneration of the BDE formulation. The approach to eliminating

the difficulty which has been successfully used in two dimensions is based on the

appropriate Green's function for the problem. The analytical benefits deriving from the use of the Green's function for 2D elastic fracture mechanics analysis will also be

reviewed. The principal finding concerns path independent integrals in cracked bodies.

Extension of the use of this GF formulation to elastoplasticity will also be highlighted.

The computational burdens for the use of the Green's function will be delineated.

The 3D elastic fracture mechanics problem formulation does not provide a

suitable Green's function, in analytical form. Recent work on the traction-BIE provides a possible numerical approach to constructing general 3D Green's functions for complex

crack geometries. The basic algorithm for the proposed approach will be given. Other

problem areas needing advances in formulation include Green's functions for cracks at and near material interfaces in both 2D and 3D. I will attempt to outline some of these

problem needs.

144
Mechanical Properties of Metal-Matrix Composites

L. C. Davis

Research Laboratory

Ford Motor Co.

Dearborn, Michigan

Metal-matrix composites (MMC) show great potential as light weight materials

for a variety of automotive applications. The high specific modulus and strength of these

materials make them suitable for disk brake rotors, connecting rods, cylinder liners, and

other high temperature components. Typically, MMC's consist of an aluminum matrix

reinforced by particles of SiC or other ceramics. The volume fraction of the reinforcing

phase is in the range 15 - 30% and the particle size is generally 1-10 |j.m. Considerable

effort has been put into developing models that accurately describe the properties of

MMC's. The inputs to the models are the properties of the constituents (e.g., moduli,

yield strength, coefficient of thermal expansion, creep rate), the volume fraction of

particles, and information regarding the phase geometry (particle shape and distribution).

Unlike fiber-reinforced composites where fiber-matrix interface properties are often

dominant, the metal-ceramic interface can usually be assumed to be ideal. Some linear

properties (for example modulus) can be estimated with sufficient accuracy by applying

rigorous bounds, either Hashin-Shtrikman or more accurately the third-order bounds of

Milton, Torquato, and others. However, for non-linear properties such as yield strength,

finite element analysis (FEA) has proven to be extremely useful.

FEA modeling at Ford Research Laboratory has included: (1) residual stresses induced by thermal mismatch between the matrix and the reinforcing particles, (2) stress- strain relation under uniaxial loading, and (3) creep of MMC's. In most cases, the FEA

results agree well with experiments done by Allison and coworkers in the Materials

Science Department of our laboratory. An interesting example of the verification of the

modeling has been the measurement by neutron diffraction (with Los Alamos National

Laboratory) of internal particle strains due to residual and applied stresses. On the other

hand, comparing FEA calculations with measurements of creep rates has not shown as good agreement, possibly demonstrating the significance of changes in matrix

microstructure in the composite relative to the unreinforced material. One of the current

limitations of the FEA approach is the restriction to unit cell models. It has not been

possible to analyze realistic particle arrays (e.g., random) in three dimensions, although

some work in two dimensions has been reported.

145
Green's Function BEM Research at Rutgers CMAS

and its Applications to the ModeUng of Mechanical and

Piezoelectrical Behaviors of Advanced Materials

M. Denda

Mechanical and Aerospace Engineering

Rutgers University

Piscataway, New Jersey

We are developing a series of boundary element methods for multiple crack

problems at the Center for Computational Modeling of Aircraft Structures (CMAS). The

methods can analyze multiple curvilinear cracks in finite bodies with elastic-plastic

deformation under cyclic loading conditions. The methods are based on micromechanics

and complex variable techniques that consist of three modules: complex variable

boundary element method (CVBEM), crack source method (CSM), and the plastic source method (PSM). A systematic use of the micromechanics tools such as the point force,

dislocation, their dipoles, and their continuous distributions enables us to build each

module independently that can be called upon as needed. The derivation of the Green's

functions and their analytical integrations are critically dependent on the use of complex

variables.

The CVBEM has been developed to accommodate the effects of the finite

boundary. It is formulated with the help of a physical interpretation of Somigliana's

identity which provides the basis of the use of the continuous distribution of point forces

and dislocation dipoles along the boundary of the finite bodies which is embedded in the

infinite matrix. The CSM has been developed to simulate a crack by the continuous

distribution of dislocation dipoles along a line. Meanwhile, the PSM has been developed

to represent the plastic deformation in terms of Green's functions; the method is essentially based on the continuous planar distribution of dislocations.

Although these techniques have been developed primarily for the damage

tolerance analyses of aging aircraft structures which contain multiple cracks, their

applications to other areas are extremely promising. Among them are the boundary

element method for anisotropic solids. For plane problems of anisotropic elasticity we are proposing to formulate the BEM based on the physical interpretation of Somigliana's

formula. The necessary governing equations are provided by the Stroh formalism which is essentially a formalism based on complex variables. The power of complex variable

techniques in the derivation of various Green's functions used in the BEM will be

demonstrated by examples including the half-plane, bimaterial plane, plane with a crack

or a hole, etc. Another promising area is the BEM for plane piezoelectricity. It turns out

that the piezoelectricity can be accommodated by extending the Stroh formalism to include the effect of electricity. Formulation of the BEM as well as the suggestion for the

derivation of the various Green's function for piezoelectricity will be given in the

presentation.

146
Electroelastic Green's Functions and their use in Piezoelectric

Material Modeling

M. Dunn

Center for Acoustics, Mechanics, and Materials

Department of Mechanical Engineering

University of Colorado

Boulder, CO 80309

We present recent progress on our efforts to model microstructural-level

electromechanical phenomena in piezoelectric ceramics using a Green's function based

approach. The phenomena studied include grain and sub-grain level electromechanical

interactions, microcracking, fracture mechanics, and the effective properties of

piezoelectric composites.

147
Materials Research and Development at Biosym

B. Eichinger

Biosym Technologies

San Diego, California

Techniques in computational physics and chemistry have advanced rapidly within

the last five years. Today, Computer Aided Molecular Design (CAMD) is widely

accepted and practiced within the pharmaceutical industry, and new drugs are emerging

that have been designed with the aid of these methods. Applications to other areas of

chemistry and physics, such as to the structure and properties of polymers and catalysts,

the band structure of metals and insulators, and the magnetic properties of metals, are now routine.

At the other end of the spectrum, computational methods are currently widely

utilized in the mechanical design and manufacturing management areas. Software for

thermal diffusion, stress, and flow analysis is widely applied in the materials industry.

However, comprehensive methods that unite the power of atomistic methods with the

practical advantages of engineering calculations do not currently exist. Biosym intends to bridge this gap.

To achieve this goal, we propose to establish a major new consortium project with

companies, research laboratories, and universities involved in materials development

production and/or utilization The project has four main objectives:

1 . To collect, validate, integrate and package the best predictive methods.

2. To produce, together with affiliated research groups, enhancements to software

which will extend current methods, strengthen their theoretical base and broaden

their applicability.

3. To combine these methods into an integrated materials design tool set suitable

for use by engineers and scientists.

4. To link these methods to currently used mechanical design and process

management models.

The project will be undertaken in several three year phases: The first phase is to begin in January 1995, with subsequent phases to be defined and established as members

and Biosym later agree. Work will be carried out at Biosym's Development Center in San

Diego by a dedicated staff of Materials Scientists and Software Engineers. The

opportunities for member companies and institutions and Biosym lie in the collaborative

nature of the project. By working with Consortium Members, Biosym consistently

develops software tools that address the specific needs of its Members. We have several

on-going projects with the National Laboratories, and others are being discussed. Biosym

also provides its academic and non-profit collaborators with the means by which their

scientific codes can be usefully applied to practical problems by non-experts.

148
Approximate Operators for Boundary Integral Equations in

Transient Elastodynamics

T. L. Geers and B. A. Lewis

Center for Acoustics, Mechanics and Materials

Dept. of Mechanical Engineering

University of Colorado

Boulder, Colorado

The integral equation for the three-dimensional transient dynamics of an isotropic

elastic medium has been known for about a century, and numerical methods for solving

the equation have been under development for about a quarter-century. However,

application has been limited to idealized problems because of the complexity of the

integral operators and the intense consumption of computational resources. The state of

development for anisotropic media is much further behind. The present situation

motivates the search for accurate approximate operators that facilitate straightforward

implementation and rapid computation. One such search is described herein, which

involves the formulation and evaluation of singly and doubly asymptotic operators for

unbounded domains. The former approach exactness at either early time or late time; the

latter approach exactness at both early and late time. Singly asymptotic operators yield

satisfactory results only in restricted circumstances, whereas doubly asymptotic operators

have proven to be quite robust.

149
Stress in Epitaxial Thin Films

Len Gray

Mathematical Sciences Section

Oak Ridge National Laboratory

Oak Ridge, Tennessee

Materials for electronics applications can be fabricated by depositing a thin layer

of one material on a substrate of a second. The boundary element method for

elastostatics is applied to model the stress state in a system for which the thin film is grown epitaxially. Boundary conditions at the material interface are developed to model

the epitaxial growth of the film. While most previous work has ignored the relaxation of

the substrate, these initial two-dimensional calculations demonstrate that the effect of the

substrate is significant. The ultimate aim of this work is to develop a dynamic model of

thin film growth. Example results for thin films of germanium deposited on a silicon

substrate will be presented.

150
A Hybrid Micro-Macro BEM for Ceramic Composites

Z. Q. Jiang and A. Chandra

Department of Aerospace and Mechanical Engineering

The University of Arizona

Tucson, Arizona

Many common engineering materials contain defects in the form of cracks, voids

or inclusions. The micro-scale interactions of these defects profoundly influence the

strength and life of these materials. This is especially true for relatively brittle materials

such as ceramics, intermetallics and ceramic matrix composites. The micro-scale features

are typically of the order of a few mms while the overall components are of the order of a few cms to a few meters. Thus, one must incorporate the micro-scale effects into a widely different macro-scale analysis.

Accordingly, the present work aims at developing a hybrid BEM formulation that can capture the effects of micro-scale features within a macro-scale computational

scheme. As a first attempt, particular attention is paid to crack-crack and inclusion-crack

interactions. Local analysis schemes capable of detailed representations of micro-features

of a problem are mtegrated with a macro-scale BEM technique capable of handling

complex geometries and realistic boundary conditions, the micro-scale effects are introduced into the macro-scale BEM analysis through an augmented fundamental

solution obtained from an integral equation representation of the micro-scale features.

Thus, the necessity for a fundamental solution, which has been traditionally viewed as a weakness of BEM, is turned into one of its strengths and is utilized as a conduit for

introducing micro-scale analytical advancements into macro-scale BEM computations.

The proposed hybrid micro-macro BEM formulation allows decomposition of the

complete problem into two sub-problems, one residing entirely at the micro-level and the

other at the macro-level. This allows for investigations of the effects of the

microstructural attributes while retaining the macro-scale geometric features and actual

boundary conditions for the component or structure under consideration. As a first attempt elastic fracture mechanics problems with interacting cracks and inclusions at close spacings are considered. The numerical results obtained from the hybrid BEM

formulation can easily be extended to investigate various micro-features (e.g., interfaces,

short or continuous fiber reinforcements, voids, inclusions) of composite materials on

macroscopic failure modes observed in numerous real-life components.

151
Stochastic Boundary Elements for Two-Dimensional Potential

Flow in Nonhomogeneous Media

I. Kaljevic

Ohio Aerospace Institute

Brook Park, Ohio

A stochastic boundary element formulation is presented for the treatment of twodimensional problems of steady state potential flow in nonhomogeneous media, that

involve a random operator in the governing differential equation. The randomness is introduced through the material parameter of the domain, which is described as a nonhomogeneous random field. The random field is discretized into a set of correlated

random variables, and a perturbation is applied to the differential equation of the

problem. This leads to differential equations for the potential and its first and second

order derivative, respectively, evaluated at the mathematical expectations of the random

variables resulting from the discretization of the random field. An approximate method is applied for the solution of these equations, by expressing the potential and its derivatives,

respectively, as the sums of functions of descending order. These solutions are introduced

into the differential equations and upon equating similar order terms, the sequences of

Poisson's equations are obtained for each order. These equations are solved using the

boundary element method to obtain the unknown boundary values of the response

variables and their respective first and second order derivatives which are then used to compute the desired response statistics. Quadratic, conforming boundary elements are used in the boundary integration, and four node quadrilateral cells are used in the domain

integration. Strongly singular terms of the boundary element matrices are obtained

indirectly, by applying a state of the unit potential over the entire contour of the object.

The singular domain integrals are calculated analytically. Analytical derivations for the

computation of the singular domain integrals are presented. Direct solution techniques are used to calculate the response variables and their derivatives, respectively. A

transformation of the correlated random variables into an uncorrelated set is performed to reduce the number of numerical operations by retaining a small number of transformed

random variables. A number of example problems are solved using the present

formulation and the results are compared with those obtained from Monte Carlo

simulations. A good agreement of the results is observed.

152
Boundary Element Method for Computer Aided Design

Guna Krishnasamy

PDA Engineering

Costa Mesa, CA

Computer aided design (CAD) is a process where a computer model of solids are

created. The computer model of a solid is usually a collection of surfaces which when

put together forms the boundary of the solid. Boundary element method (BEM) is a technique which uses the boundary of the solid and the associated boundary conditions to solve for the unknown displacements and traction's on the boundary. From this brief

description of CAD and BEM it is clear that the boundary of a solid plays an important

role in both of these areas. This makes BEM an attractive analysis tool on computer

models based on boundary representation.

The interface between CAD and BEM has posed several challenging problems.

Generally the elements used for BEM analysis need to be congruent. Such a restriction on

the elements will produce a large number of nodes and therefore make the BEM analysis

expensive, even for small problems. Approximation of the geometry by quadratic or cubic representation, as is usually done for BEM, can also require a large number of

elements to define the geometry. To overcome these problems several new ideas have

been researched at PDA Engineering. This includes (a) generating elements which are not congruent (b) the description of the geometry and the approximation of boundary

variables are separated (c) several analysis are performed so that additional nodes are introduced only in regions where the results vary significantly (adaptive analysis). These

ideas have been implemented in TEAM (Trimmed Element Analysis Method) a commercial engineering software based on BEM for CAD systems. This is developed at PDA Engineering.

When formulating the boundary integral equations for a linear elastostatic

problem, certain requirements on the smoothness of the boundary displacement is required. At the collocation point the displacement should be at least Holder continuous.

This requirement on the displacement is satisfied if the elements are congruent which is the case in the traditional BEM. Such a requirement on the mesh will make the analysis

very expensive. In TEAM non conforming elements are generated to start with and the

requirement on the continuity of the displacements at the collocation point is enforced

after the system of equations are generated.

The underlying math forms which describe the geometry in CAD and those that

describe the geometry for BEM are very different. The boundary in a computer model are usually composed of nurbs, splines and analytic surfaces which tend to describe the

boundaries exactly. On the other hand for BEM piece wise quadratic or cubic

isoparametric elements are used. Since computer models are usually made of complex

surfaces, it takes many boundary elements to represent the geometry and therefore

increases the size of the model. By uncoupling the function that describe the geometry

and those that describe the variation of the boundary variables, it is possible to have far

153
fewer geometric elements to describe the geometry. Moreover the number of variables

used to describe the boundary variables can be controlled since they are independent of

the geometry.

154
Overview of BEM Research in Mukherjee's Group at Cornell

Subrata Mukherjee

Department of Theoretical and Applied Mechanics

Cornell University

Ithaca, New York

Ongoing research in Prof. Subrata Mukherjee's group at Cornell, involving the

boundary element method (BEM), is briefly summarized below.

1. The Boundary Contour Method (BCM)

The conventional BEM, for linear problems such as potential theory or linear

elasticity, requires numerical evaluation of surface integrals for 3-D and line integrals for

2-D problems. A novel approach recently developed at Cornell, called the Boundary

Contour Method (BCM), achieves a further reduction in dimension. Through a very

interesting application of Stokes' theorem, the integrals are transformed so that only line

integrals have to be numerically evaluated for 3-D linear problems. 2-D problems require no numerical integration at all. Numerical results from the BCM have been obtained for 2-D and 3-D elasticity

problems. The idea is quite general and can be applied to other linear problems such as acoustics and elastodynamics.

2. Micro-electro-mechanical (MEM) structures Micro-electro-mechanical (MEM) systems consist of integrated movable

microstructures with electronics. Typical MEM structures consist of arrays of thin beams

or plates with cross-sections in the order of microns and lengths in the order of ten to a hundred microns. Sensors and actuators using MEM technology are already finding

diverse applications in the automotive and medical industries. Many other promising

applications of this technology are being currently developed.

A Cornell project on simulation and design of MEM structures involves

Professors Mukherjee and Noel MacDonald (EE) and Dr. Ramesh from Xerox

Corporation. A hybrid BEM/FEM approach is being implemented with the BEM being

used to solve the exterior electrostatic and the FEM the interior (conductors and

dielectrics-beams and plates) elastostatic problem. The project also involves sensitivity

analysis, optimal design and reliability analysis of MEM structures, in addition to simulation.

3. Optimal Design of Manufacturing Processes

The goal here is to carry out optimal design of processes - such as extrusion or sheet metal forming - rather than merely the design of product shapes. These problems

are extremely challenging and involve both material and geometric nonlinearities. The

starting point here is a nonlinear BEM formulation that was first proposed by Mukherjee

and Chandra about ten years ago. BEM simulation and sensitivity analysis of these large

155
deformation processes are being carried out. Sensitivities are rates of change of response

quantities such as stresses or displacements in an extruded product, with respect to design

variables such as geometrical parameters that determine the shape of the die. Integral

equations for sensitivities are obtained by direct analytical differentiation of the governing

BEM equations. These are then used in an optimization algorithm to optimize the

process, e.g. to determine optimal die shapes for extrusion.

156
Simulation of the Electrochemical Machining Process Using a 2D

Fundamental Singular Solution

A. Neid

Process Physics Laboratory

General Electric Company

Schenectady, New York

Electrochemical machining is a process for the removal of material from the

surface of an electrically conductive metal or alloy by anodic (workpiece) dissolution

when high current densities distributions are generated on the anode (workpiece).

The major features of the ECM process are illustrated in Fig. 1. hi this figure,

which shows two flat electrodes for simplification, metal is machined from the workpiece

(anode) due to the high current densities generated on it by the electric current flow

between the cathode and anode by ion transport. Even for such a simple geometry, a non- uniform current density distribution is produced on this workpiece creating a non-flat and

non-parallel final shape in the workpiece. Hydrogen gas is also generated by the process

which tends to collect at the cathode having an effect on the final shape of the workpiece.

In addition to the removal of heat from the tool and the workpiece, the electrolyte flow in

the narrow gap separating the electrodes sweeps away the dissolution products and the

hydrogen downstream.

Since the Electrochemical Machining (ECM) process is essentially a surface

phenomenon, a solution is readily obtained by employing classical potential theory. By

using a 2D Green's function or fundamental singular solution, the boundary value

problem for ECM process can be reduced to the solution of a boundary integral equation.

When the surfaces are discretized, a numerical solution for the potential and the

derivatives at the surface being machined can be derived. The boundary element method

(BEM) developed by Brebbia was used to solve the field equations for this nonconventional machining process. The electrochemical anodic reaction was furnished by

Faraday's Law, which provided the relationship for the rate of dissolution of the

workpiece. This 2D ECM process model was developed for the purpose of simulating

the ECM process. When this solution was employed, the airfoil shapes obtained during

the machining of compressor blades can be determined by computing the mass of metal

removed by the machining process.

Since the ECM process is essentially a moving surface phenomenon which

requires a precise determination of the 2D current density distribution on the surface of

the workpiece, the boundary element method was selected as the preferred simulation

tool. The boundary element method was shown to be a powerful computational tool for

determining the history of the current density and electrolyte velocity distribution on the

surface of the airfoil. Using the boundary integral equation method in conjunction with the appropriate

2D fundamental singular solution, many complicated geometries can be conveniently

represented by a simple 2D contour of the ECM process model consisting of the

157
boundaries enclosing the tool, electrolyte gap, insulation, and workpiece. The advantage

of this formulation is that only the boundary composed of the tool, the workpiece, and

insulators need be discretized to construct a 2D model without the necessity of

discretizing the interior of the electrolyte between the tool and workpiece as would be the case if the finite element method were employed. This analytical modeling method is most appealing since in the electrochemical machining process we are only interested in

the reactions and dissolution at the surface of the anode (workpiece). This feature has a

distinct advantage for simulating the electrochemical machining of airfoils using the

opposed electrode tooling method. Fig. 2 shows a typical mesh geometry using this

method for opposed tool machining the concave and convex sides of the airfoil

simultaneously. Note that only a surface mesh is required to describe the tooling,

electrolyte inlet and outlet, and the airfoil workpiece shape. This figure shows the initial large machining gap between the airfoil preform and the cathode tooling when the

process is started. Fig. 3 shows the final position of the tooling during the ECM

simulation. The machining gap at the final stages of machining is very small and

contains a rapidly flowing conductive electrolyte with a temperature distribution along

the flow path. Fig. 4 shows a typical result for the simulation just before the process is stopped. The solution shows that the current density varies significantly from the leading

to trailing edge of the airfoil. hi this paper, the theory and numerical application to machining 2D airfoils is described. : . -

158
3-D Boundary Element Formulation of Anisotropic Elasticity with

Gravity

E. Pani and B. Amadei^

Department of Civil Engineering

University of Colorado

Boulder, CO 80309

This paper presents a new 3-D boundary element formulation for linear and

anisotropic elastostatics with gravity. The material property within the 3-D finite domain

is transversely isotropic with any orientation. Green's functions for such anisotropic

media are derived in exact closed-form. Particular stresses and displacements related to the gravity body force are also presented for generalized anisotropic media. The boundary

integral equation with only weak singularities is formulated based on the principle of

superposition. No volumetric integration associated either with the Green's functions for

anisotropic media or with body forces is involved. This new 3-D boundary element

formulation is programmed, and several numerical examples are presented and checked

with existing closed-form solutions. It is found that even with very coarse discretization,

very good results can be obtained.

1 Research Associate

^Professor

159
Boundary Element Modeling for Nondestructive Evaluation using

Ultrasonic and Electromagnetic Fields and the Idea of Discretized

Green's Functions

F. J. Rizzo

Dept. of Aerospace Engineering and Engineering Mechanics

Iowa State University

Ames, Iowa

This presentation describes several computer models of nondestructive evaluation

(NDE) procedures to detect flaws and cracks in materials and structures, particularly

aircraft structures. The models depict inspection situations which use ultrasonic and

electromagnetic waves and eddy current fields as energy sources. Boundary elements are the main computational vehicle for the models; the character of the governing equations

and the regions over which the fields are defined are such that boundary elements are a good choice. The objectives of these models are to confirm data from physical

experiments and to provide criteria for the design of new configurations for physical

inspection. The goal is to improve our ability to detect and characterize flaws and defects

and ultimately improve criteria for safety and reliability of structures containing flaws and

defects.

Specifically, boundary integral formulations for scattering of ultrasonic waves from cracks and the perturbation of eddy current fields by cracks are presented, along

with formulations for scattering of electromagnetic waves from one or more volumetric

obstacles. Some comparison data are given to verify our models on test problems.

However, emphasis in this presentation is on viewing the steps in the formation

and solution of the algebraic equations, which arise following the boundary element

approximation to the boundary integrals, in a modular and systematic fashion. The idea

here is to decompose a scattering or eddy current field problem into pieces, which may be

discretized and solved independently. Some useful features of this idea are to take

maximum advantage of the limited memory of a single workstation, to do a problem in

steps on that workstation, or to do the decomposed-problem pieces on several

workstations at once as a kind of parallel processing.

More generally, it is well known that if the region-dependent Green's function is available for a given domain, linear differential operator and type of boundary data, the

solution to a boundary value problem for those boundary data is just a matter of

quadrature over the boundary. On the other hand, when only the free-space (singular)

region-independent part of the Green's function is known (as is usually the case), the

problem is, of course, much more difficult. Indeed, this is the situation which requires

computational methods as we know them; e.g., boundary discretization, element

integration for each collocation point, formation and solution of algebraic equations (or

matrix inversion), i.e., the boundary element method (BEM), to provide unprescribed

boundary data. Only after all this work comes the final (fairly simple) step of boundary

160
quadrature, and this step is comparable to a case when the region-dependent Green's

function is known.

However it is possible, for many problems, not only in NDE, to view several of

the products of matrices and matrix inverses which arise in the BEM as discretized

Green's functions. These may be profitably identified and stored for repeated use. Once

this is done, solution of a given problem for a particular domain should be

computationally equivalent to having the analytical Green's function for that domain,

from the beginning, and getting numerical values for boundary value problems from it. Problems involving several surfaces, but one of which is the same for each problem,

would be especially well-suited to the discretized Green's function concept.

Going further with this idea, suppose one is faced with a class of problems where

one region differs from another according to the value of a parameter, e.g., shape of a

hole, inclination of a surface-breaking crack, density of an inclusion, etc. We assume that

each shape will carry loads or disturb certain specified ultrasonic or electromagnetic

fields. It is possible and practical now, with current data-storage technology, to proceed

as follows, (a) Discretize the surfaces of many domains based on the mentioned

parameter, and form and invert and store the relevant boundary element matrices. Then at a later date, perhaps in another place, (b) pick a region, specify loads or nominal fields,

e.g., an incident wave, plus the coordinates of the solution point in the region, and call up

the precomputed matrices, thought of now as discretized Green's functions for that region.

Finally, do the matrix multiplication involving a precomputed matrix and the specified

data (i.e. the quadrature) and the problem is solved. Note that little if any knowledge of

computational schemes, let alone the boundary element method, is required for step (b).

The real work has been done in (a) in forming and identifying the discretized Green's

functions.

The advantages with this approach are that expertise in computational methods

and the labor of discretization plus the extensive computations are done in advance. In

effect, a kind of modern handbook or library of commonly useful configurations is formed. Thus the labor, expertise, and extensive computing which goes into the library are traded simply for storage space which is now available economically in great

quantities with fairly convenient and quick retrieval. In summary, this presentation is concerned with three things: (1) some specific

computational models for NDE using BEM, (2) modular computing with the idea of

discretized Green's functions, and (3) the idea of a handbook or library of what are in essence Green's functions in discretized form.

161
Fracture and Crack Growth of Polymeric Composites for Use in

Dynamic Applications

R. A. L. Rorrer

Advanced Materials Research

Gates Rubber Company

Denver, Colorado

Examination of the technical concerns and issues in failure of polymeric,

specifically elastomeric, composites used in dynamic applications. One of the dominant

failure mechanisms is related to the initiation and propagation of cracks in the composite,

which often lead to subsequent catastrophic failure. With multi-component composites,

the design issues center around material selection, engineering of the component

materials and interfaces, and geometry changes. The overall goal is to extend the useable

life of the composite by elimination or reduction of crack propagation as a dominant

failure mode. Examples of typical failure paths will be presented and discussed for an

elastomeric composite.

162
Inverse Elasticity Solutions Using Boundary Elements

S. Saigal and L. M. Bezerra

Department of Civil and Environmental Engineering

Carnegie Mellon University

Pittsburgh, Pennsylvania

Given a set of measurements, either at the boundary of or within an object, the

inverse elastostatics problems (lESP) of: (a) detection of a geometric flaw, and (b)

reconstruction of missing boundary conditions are considered. Such inverse problems

usually start with an initial guess of the missing information and proceed towards the

final configuration in a sequence of iterative steps. The prevalent formulations using

finite elements will require a remeshing of the object corresponding to the revised

configuration in each iteration leading to a computationally expensive and cumbersome

procedure. The boundary element approach presented here does not require such a remeshing. Further, the solutions for lESPs are generally driven by boundary data which

is predicted more accurately using boundary elements.

The inverse problem is written as an optimization problem with the objective

function being the sum of the squares of the differences between the measured

displacements and the displacements computed for the assumed configuration. The

geometric condition that the flaw lies within the domain of the object is imposed using

the internal penalty function approach in which the objective function is augmented by

the constraint using a penalty parameter. The missing traction distribution is assumed to

lie within a certain portion of the boundary of the object and this constraint is imposed

using a similar approach. The unknown geometric flaw and the unknown traction

distribution are defined in terms of load and geometric parameters to limit the number of

design variables for a faster convergence of the optimization procedure. A first-order

regularization procedure is implemented to modify the objective function in order to minimize the numerical fluctuations that may be caused in the numerical procedure due to errors in the experimental measurements. The constrained optimization problem is first transformed into an unconstrained optimization problem using the internal barrier penalty

functions approach. The design sensitivities required in the numerical optimization

procedure are obtained by the implicit differentiation approach.

A series of numerical examples involving the detection of circular and elliptical

flaws of various sizes and orientations, and the reconstruction of linear, parabolic, and

trigonometric boundary tractions are solved using the present approach. For the flaw

detection problems, the effectiveness of a two step procedure that employs fewer design

variables during the initial iterations and the complete set of design variables

during later iterations is demonstrated. The utility of the present approach in guiding a non-destructive evaluation in identifying the measurements that must be made to allow an

effective reconstruction of the missing data is shown. The present approach for

reconstruction of missing tractions considers the location of the traction distributions as

163
an unknown in addition to its magnitude and distribution. A number of example problems

are considered to demonstrate this feature. The effect of Gaussian errors in the sensors for

experimental measurements on the numerical predictions are also studied. A good

performance of the present optimization based approach is observed in the lESPs of both

flaw detection and boundary traction reconstruction.

164
Elastic Green's Functions for Anisotropic Solids

V. K. Tewary

Materials Reliability Division

National Institute of Standards and Technology

Boulder, Colorado

A solution of the time dependent Christoffel's equation for anisotropic solids and

the corresponding Cauchy problem will be described in terms of delta functions in the

space of time and slowness vectors. This solution is used to obtain an integral

representation of the three-dimensional elastodynamic and elastostatic Green's function

for anisotropic solids. The representation contains a product of two delta function and is, therefore, computationally very convenient. The relation between the delta function

representation and the Fourier/Laplace representation will be discussed. The application

of the Green's function to solution of static as well as dynamic boundary value problems

by the virtual force method will be described. The correspondence between the virtual

force method and the boundary element method will be discussed.

165
Lattice Green's Functions for the Study of Defects in Solids

R. Thomson

Materials Science and Engineering Laboratory

National Institute of Standards and Technology

Gaithersburg, Maryland

I will give a brief overview of the techniques for using lattice Green's functions

for the study of defects in solids [1] and in particular for the study of dislocations and

cracks in solids. A general summary of the technique will be given with some

characteristics of the numerical methods. In addition, I will briefly review our use of the

methods for the study of the ductility of solids.

[1] This method was developed by V. Tewary and is detailed in his review, Adv. in

Phys., Vol. 22, p. 757 (1973).

166

o U.S. GOVERNMENT PRINTING OFFICE: 1997-417-392/67919
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) devel- oped 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.

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 under 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

bimonthly for NIST by the American Chemical Society (ACS) and the American Institute 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 treatment 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 in support of the efforts of private-sector

standardizing organizations.

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 nongovernment). 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

325 Broadway

Boulder, Colorado 80303-3328

Official Business

Penalty for Private Use, $300