C. R. Acad. Sci. Paris, t. 327, Sbrie II b, p. 1215-1221, 1999 Mecanique des solides et des structures/Mechanics of solid and structures

A general boundary-only formula for crack shape sensitivity of integral functionals Marc BONNET

Laboratoire de mkanique des solides, CNRS UMR 7649, kcole polytechnique, 91128 Palaiseau cedex, France

E-mail: bonnet@lms.polytechnique.fr

(Recu le 19 avril 1999, accept6 le 25 mai 1999)

Abstract. This note presents, in the framework of three-dimensional linear elastodynamics in the time domain, a method for evaluating sensitivities of integral functionals to crack shapes, based on the adjoint state approach and resulting in a sensitivity formula expressed in terms of surface integrals (on the external boundary and the crack surface) and contour integrals (involving the direct and adjoint stress intensity factor distribu- tions on the crack front). This method is well-suited to boundary element treatments of e.g. crack reconstruction inverse problems. 0 1999 Academic des science&ditions scientifiques et medicales Elsevier SAS

elastodynamics / shape sensitivity I crack / inverse problem / adjoint state

D&iv&es de fonctionnelles in tdgrales dans des perturbations de fissures exprimces par intbgrales de front&e

RCsumC. Cette note propose, dans le cadre de l’e’lastodynamique lineaire tridimensionnelle en variable temporelle, une me’thode d’evaluation de la sensibilite de fonctionnelles integrales a des perturbations deJissures, qui repose sur la de’jinition d’un &at adjoint et permet d’exprimer la sensibilite’ en termes d’integrales de su$ace (frontier-e exteme etjssure) et de contour (contribution des distributions directes et adjointes de facteurs d’intensite’ de contraintes sur le front de fissure). Ce resultat est par exemple bien adapt6 aux techniques de reconstruction de jissures fond&es sur les elements de frontier-e. 0 1999 Academic des scienceskditions scientifiques et medicales Elsevier SAS

elastodynamique / variation de domaine /pssure /probEme inverse / &at adjoint

Version francaise abr&gke

l?valuer la sensibilite de fonctionnelles integrales par rapport a des perturbations du domaine est notamment utile dans des situations a domaine inconnu (problemes inverses), ajustable (optimisation) ou variable. Cette note propose une formulation de la sensibilite a des perturbations de fissures, developpee dans le cadre de l’elastodynamique lineaire tridimensionnelle en variable temporelle et fondee sur la methode de l’etat adjoint.

Note prCsent6e par Huy Duong BUI.

1287-4620/99/032701215 0 1999 Acadhie des sciences&ditions scientifiques et mkdicales Elsevier SAS Tow droits rCserv& 1215

M. Bonnet

On consider-e ainsi un corps Clastique borne Sz c lR3, de front&e exteme S et contenant une fissure r. Les deplacements U, deformations E et contraintes B verifient les equations de champ (1) et les conditions aux limites et initiales (2), l’ensemble constituant le problkme direct.

La sensibilite aux perturbations de fissure de fonctionnelles integrales de la for-me (3) est analysee en considerant des transformations de domaine qui laissent S fixe, c’est-a-dire de la forme x7 = x + ~0( x), oti q est un parametre cintmatique, et la vitesse de transformation 0(x ) verifie

B = 0 sur S. Si j =x7 + Vf.0 designe la derivee lagrangienne de f, les derivees d’integrales gtneri- ques sont donntes par (4).

Adoptant les methodes du controle optimal, le probleme direct est introduit comme contrainte explicite dans le lagrangien (5), dont la deride dans une perturbation de fissure est trouvee Cgale a (6). Notre choix de derivation lagrangienne resulte du fait que la d&-i&e partielle (Vu ),lr est singuliere comme de3” le long du front ar tandis que Vu* et Vu ont la m&me singularitt d-“2 (d : distance a Jr). Ensuite, l’etat adjoint ti,-, choisi de facon a annuler les termes en u* et p’, est defini par les conditions aux limites et jinales (8). Cela conduit a l’expression (9) de la sensibilite. Cette demiere expression pourrait &tre immediatement convertie en integrales de front&es a l’aide de I’identite (10) en l’absence de singular& en front de fissure. Ici, il est ntcessaire d’isoler un voisinage tubulaire D, = {XI dist (x, ar) I E} du front ar, d’utiliser les expressions bien connues (12) des parties singulihes des champs et d’evaluer la contribution non nulle de l’integrale sur Z;: quand e + 0. On obtient ainsi, resultat principal de cette note, l’expression (13) de la sensibilite en termes d’integra- les de surface et de contour.

L’equation (13) s’applique a des fissures et des vitesses de transformation regulieres assez g&n&ales. Elle est par exemple bien adaptee aux techniques de reconstruction de fissures fondees sur les elements de front&e. Soulignons qu’elle conceme des perturbation de configurations jixes de jssures, non des propagations. Bien entendu, elle s’applique, moyennant les modifications Cvidentes, a l’eastostatique et l’elastodynamique en domaine frequentiel. Par exemple, j( r) est l’energie potentielle a l’equilibre pour le choix particulier vu = - (pD.u )/2, vp = ( uD.p )/2 dans l’equation (3), l’etat adjoint &ant alors u” = ( l/2 ) U, soit Z?, = K, /2, etc. Dans (13), le facteur de ( 0.v ) ( s ) est bien, comme attendu, l’op- pose du taux de restitution d’energi G( s ), tandis que le facteur de ( 0.n ) ( s ) est l’equivalent 3-D de l’integrale J, [ 11, 121. L’integrale H (13) apparait Cgalement comme un cas particulier de (13).

1. Introduction

The need to compute the sensitivity of integral functionals to variations in shape parameters arises in many situations where a geometrical domain plays a primary role; shape optimisation and inverse problems are the most obvious, as well as possibly the most important, of such instances. In addition to numerical differentiation techniques, shape sensitivity evaluation can be based on either direct differ- entiation or the adjoint variable approach. This paper is focused on the latter.

The object of this note is to present a method for evaluating crack shape sensitivities of integral functionals, in the framework of three-dimensional linear elastodynamics in the time domain. This method is based on the adjoint state approach and results in a sensitivity formula expressed in terms of surface integrals (on the external boundary and the crack surface) and contour integrals (involving the direct and adjoint stress intensity factor distributions on the crack front). The sensitivity formula is thus well-suited to the use of boundary element methods (BEMs), which are quite often used for solving, for example, crack reconstruction inverse problems, ([ 1, 21).

1216

A general boundary-only formula for crack shape sensitivity of integral functionals

Consider an elastic body Sz G lR3 of finite extension, externally bounded by the closed surface S and containing a crack l7 The displacement u, strain E and stress ts are related by the field equations:

diva-pii=O ~=C:E e=$(Vu+Vru) inQ (1)

(C: fourth-order elasticity tensor). In addition, displacements and tractions are prescribed on the portions S,, and S,, = S - S,, of S, the crack surface r is stress-free and initial rest is assumed:

u=rP(onS,), P=$(onS,), P=O(onr), u=ti=O(in Q,att=O) (2)

where p = cr.12 is the traction vector, defined in terms of the outward unit normal n to Q. The above conditions define the direct problem.

Let us introduce the following generic objective function:

which is encountered for instance in minimisation-based algorithms for solving the inverse problem of crack detection - ( ur, pr ) refers to the solution of problem (l), (2) for a given crack configuration. A boundary-only expression for the derivative of ,J( ZJ with respect to crack perturbations is sought.

Any sufficiently small perturbation of r can be described by means of a domain transformation which does not affect the external boundary S, i.e. of the form xv = x + ~0( x), where r is a time-like parameter and the transformation velocity field 8(x ) is such that 8 = 0 on S. Denoting by f =fq + Vf.0 th e L agrangian derivative of some field variablef, the derivatives at q = 0 of integrals over generic domains V and surfaces Z have the well-known form:

-& s s V

fdV= (f* +fdiv 0) dV -& fdS = (f* +fdiv, 0) dS V s s z z

where divs 8 = div 0 - n.Vf2.n is the surface divergence of 8 (n: outward unit normal vector). Also, recall that (Vu )* = Vu’ - Vu.VB.

2. Adjoint problem and domain integral formulation

For the present purposes, it is convenient to use an optimal control approach, whereby the variables in the objective function J( u,p, r) are formally considered as independent ones and the direct problem (i.e. the fact that one actually has u = uy, p =pr) is treated as an explicit constraint. The following Lagrangian is thus introduced:

U(u, ii,p,f-A r> = J(u,P, T> + ss

T [a:VrZ+pii.ti]dVdt 0 D

T ss s (u-uD).@lsdt- T ss

T - p.ii ds dt - o

ss pD .il CL!3 dt

u 0 S” O 5 (5)

where (r&p”), the test functions of the direct problem in weak form, act as Lagrange multipliers.

1217

M. Bonnet

Using formulas (4), noticing that (C: Vu*): Vti = ~:VU* and ignoring the terms containing ti*,p* arising in this calculation (they merely reproduce the direct problem constraint and thus vanish), the total derivative of 9 with respect to a given domain perturbation is given by:

=@(u, ti,p,B, r> =

T

+ ff

o d[a:Vz~+pii.ti]divt3- [o.ViZ+B.Vu]:VB}dVdt

+ s

r[Vy/.B + I,Y div, 01 dS

The partial derivative ( Vu )ihhas generally a dP 3’2 singularity along the crack edge ar, while Vu* and Vu have the same d- derivatives u*.

singularity (d: distance to ar); this explains our using Lagrangian

Since the initial conditions u( . , 0) = ti(. , 0) = 0 hold for any location of the assumed defect, one should assume u’( . , 0) = ti*( . , 0) = 0 as well. Equation (6) is thus altered by taking into account the identity:

f T

zi*.P dt = ( zi*.ti - t.u*) If = T + 0 f

T

u*.ii dt (7) 0

Now, the multipliers zZ,p are chosen specifically so that all terms containing u* andp* in equation (6) combine to zero for any u* and p*. The weak formulation of an adjoint problem, of a form similar to the direct problem in (5) but with u*, p* now acting as test functions, is thus defined. Its adjoint solution ti, Pr is therefore found (by analogy to (5)) to solve equations (1) together with the following boundary and jinal conditions:

p=-g(onS,) ti=g(on&) p=O(onP) zZ=&=O(inQ,att=T) (8)

Finally, equation (6) allows one to express the derivative of J in terms of the direct and adjoint solutions:

= ff ,: .{ [ar: VZZ~ + pii,z~~] div 0 - [ a,.Vti, + gPVu,-] : V8)

+ f [Vy/.B+ y/div,@] dS I- (9)

1218

A general boundary-only formula for crack shape sensitivity of integral functionals

3. Sensitivity in terms of boundary integrals

The next step is to find an equivalent form of equation (9) involving only boundary integrals. Any elastodynamic states (u, a) and (u”, 5) satisfying initial and final rest conditions, respectively, verify the identity:

s ;I[ rx E + piii] div ~9 - [ a.Vti + B.Vu]: VS} dt (10)

=

s

T

div ([a: El--pJ.i] 8- [o.Vti +B.Vu].B) dt 0

However, the well-known singular behaviour of strains and stresses at the crack front prevents a direct application of the divergence formula to equation (9) for the entire cracked domain 0. To circumvent this difficulty, the body 52 is partitioned into D = 52, u (D, W), where D, = {xldist (x, ar) I E} for some sufficiently small E > 0 is a tubular neighbourhood of the crack front ar bounded by the tubular surface ZE. Further, put Sz, = Om, and r, = m,. Upon introducing this splitting into equation (9) applying the divergence formula for the contribution of Q2, and invoking boundary conditions (23, 8,), one obtains:

- &-ii,] ( 0.n ) - brVtir +PrVur] .fI} dS dt

T

+ ss o r[ar: VS, - pti,.&-] ( 0.n ) dS dt E

iir + pipa,] div 0 - [ ar.Vtir + B,Vu,] VS} dV dt

+ f

r[Vv.e + I// div, e] ds (11)

where n is the outward unit normal to Q2, and fl =f( x+) -f( x-) (discontinuity off across I). Now we seek the limiting form of equation (1 l), when E + 0. Recall the well-known expansions of

the mechanical fields near thefied crack front, isotropic elasticity being assumed (v: Poisson ratio, ~1: shear modulus):

1 UT = zy; ~[~~(s,f)cos~(3-4s-cos8)+K,(s,t)sin~(4v-l+3cos8)] +0(r)

=&-,e,s)+o(I)

1 Ue = z;r ~~[-K,(s,r)sin~(l-4v-3cos8)+K,,(s,r)c0s~(4v-5+3c0se)] +o(r)

2 &II( s9 t > u, = P v- &sin$+O(r)=z&r,8,s)+o(r)

and similarly for ti with stress intensity factors k r,,r,m; (r, e> denote local polar co-ordinates, attached to a point x(s) of 13r characterized by its arc length S, in the plane orthogonal to 6’I’ and

1219

M. Bonnet

emanating from x(s), and z is such that (r, 8, z) define cylindrical coordinates. The singular displacements defined by equations (12) produce direct and adjoint strains and stresses, denoted (.?,dF) and (.?,8). S’ mce by virtue of these expansions cr: E” = 0( l/r), the integral over D, vanishes in the limit (( dV = Y( 1 + 0( r) ) dr d0 dr ) in 0,). Moreover, it can be verified that [as: Vu”‘] = O(d), and hence that the integral over r, tends to the corresponding convergent integral over r, when E + 0. Finally, under mild smoothness assumptions on the closed curve ~?r and the velocity field 0, one has:

{br: VEr-P z&-.i&-] ( 8.n ) - [pj-.Vd, + j7~VU,] .e} ds dt

= ss s : ar R {[&vaS](e.t2)(~)- [pS.v~S+~S.vuS].e(s)}Ededsdt+0(Ez’2) --R

The integral in the right-hand side, which yields a finite contribution as the radius E of the tubular neighbourhood goes to zero, can be evaluated in a straightforward way using the expansions in (12). This last calculation results in the following expression for/(r), which constitutes the main result of this note:

where v( s ) denotes the unit outward normal to Jr lying in the tangent plane to r at x( s ). To ensure that formula (13) can actually be evaluated using only the boundary traces of the direct and adjoint solutions, the bilinear form C: Vzi must be expressed in terms of V,u, V,zZ, taking p =a = 0 into account in the process:

a:Vu”=p &div, udivsti+~(V,u+V~~):(Vsti+V$Z)-(n.V,u).(n.V,u”)

4. Discussion

Equation (13) holds true for general smooth crack shapes and transformation velocity fields. It provides an attractive computational tool, in conjunction with BEMs and gradient-based optimisation algorithms, for solving inverse or optimisation problems for unknown or adjustable cracks. Both the direct and the adjoint problems stated in this note can be solved by usual BEM techniques [3], including the evaluation of the direct and adjoint SIF distributions [4]. Other applications of shape sensitivity analysis using the adjoint state approach in connexion with BEMs can be found in e.g. [5-lo]. The present approach is certainly applicable to integral functionals somewhat more general than those of the form (3). Integral functionals defined in terms of domain integrals over the body R still lead to the boundary-only sensitivity expression (13) but body force distributions appear in the definition of the adjoint problem.

1220

A general boundary-only formula for crack shape sensitivity of integral functionals

It is important to stress that equation (13) provides the sensitivity of an integral functional to a perturbation of a jixed crack configuration, not to a crack propagation; hence the use of expansions (12) is valid for a crack which does not physically propagate.

Equation (13) is also applicable, with straightforward modifications, to elastostatics and elastody- namics in the frequency domain. For instance, in elastostatics, f(r) is the potential energy at equilibrium for the particular choice q,, = - (pD.u )/2, ul, = ( uDLp)/2 in equation (3). For this special case, the adjoint solution turns out to be ti = ( l/2) u, i.e. K, = K, /2, etc. In equation (13), the factor of ( 0.v ) ( s ) turns out to be, as expected, minus the energy release rate G( s ), i.e. minus the J1-integral, whereas the factor of ( 0.n ) ( s ) is the 3-D generalisation of the &integral [ 11, 121 Finally, with the choice S, = S, s, = (75 and qP =pDA - u.p, where fi,p are the boundary traces of a pre-selected auxiliary elastodynamic state with$nal homogeneous conditions, one finds that ti = u^ and that the factor of ( 8.v )( s ) in (3) is the 3-D generalization of the so-called H-integral [ 131.

References

[l] Mellings S.C., Aliabadi M.H., Flaw identification using the boundary element method, Int. J. Num. Meth. Eng. 38 (1995) 399419.

[2] Nishimura N., Application of boundary integral equation method to various crack determination problems, in: Aliabadi M.H. (ed.), Dynamic Fracture Mechanics, chap. 7, Comp. Mech. Publ., Southampton, 1995.

[3] Bonnet M., Boundary Integral Equations Methods for Solids and Fluids, John Wiley in press. [4] Fedelinski P., Aliabadi M.H., Rooke D.P., Boundary element formulations for the dynamic analysis of cracked structures,

Eng. Anal. Bound. Elem. 17 (1996) 45-56. [5] Aithal R., Saigal S., Shape sensitivity analysis in thermal problems using BEM, Eng. Anal. Bound. Elem. 15 (1995)

115-120. [6] Bonnet M., BIE and material differentiation applied to the formulation of obstacle inverse problems, Eng. Anal. Bound.

Elem. 15 (1995) 121-136. [7] Burczydski T., Application of BEM in sensitivity analysis and optimization, Comp. Mech. 13 (1993) 2944. [S] Burczydski T., Kane J.H., Balakrishna C., Shape design sensitivity analysis via material derivative - adjoint variable

approach, Int. J. Num. Meth. Eng. 38 (1995) 2839-2866. [9] Choi J.O., Kwak B.M., Boundary integral equation method for shape optimization of elastic structures, Int. J. Num. Meth.

Eng. 26 (1988) 1579-1595. [lo] Merit R.A., Differential and integral sensitivity formulations and shape optimization by BEM, Eng. Anal. Bound. Elem. 15

(1995) 181-188. [ll] Budiansky B., Rice J.R., Conservation laws and energy release rates, ASME J. Appl. Mech. 40 (1973) 201-203. [12] Bui, H.D., Mtcanique de la rupture fragile, Masson, 1978. [13] Bui H.D., Maigre H., Facteur d’intensitt dynamique des contraintes tir& des grandeurs micaniques globales, C. R. Acad.

Sci. Paris, strie II 306 (1988) 1213-1216.

1221