- Email: [email protected]
Three-dimensional flaw identification using inverse analysis
Int. J. Engng Sci. Vol. 34, No. 4, pp. 453-469, 1996
Pergamon, 0020-7225(95)00125-5
Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0020-7225/96 $15.00+ 0.00
THREE-DIMENSIONAL FLAW IDENTIFICATION INVERSE ANALYSIS
USING
S. C. MELLINGS and M. H. A L I A B A D I t Wessex Institute of Technology, Ashurst Lodge, Ashurst, Southampton SO40 7AA, U.K. Abstract--In this paper, a three-dimensional crack identification method is developed. The Dual Boundary Element Method (DBEM) is used to compute potential and flux values for known crack identities and first order optimization is used to minimize the error due to approximations made to the true crack identity. Design sensitivities are computed with a derivative formulation of the Dual Boundary Element Method. The accuracy of these sensitivity values is investigated by a comparison with values obtained using a series of finite difference approximations. The flaw identification method developed is then demonstrated for the identification of a simple internal crack.
1. INTRODUCTION In maintenance of structures, components are often tested in order to detect any flaws. The techniques used can give an approximate size of the flaw and if the flaw is considered too large, the component is usually discarded. However, this procedure causes some components which do not present a risk of failure to be unnecessarily wasted. On the other hand, some smaller cracks may present a greater risk of failure due to their position. If a technique could be employed to identify the actual flaw shape and size, then it may be possible to make a more reliable predictio:a for the lifetime of the component. The methods used to detect internal flaws are known as Non-Destructive Testing (NDT) methods. In most of these methods, a field is applied to the component and responses at various sensor points are measured. The presence of the flaw will cause a distortion in the field which can be measured and used to approximate the flaw size. Common fields used in such NDT methods are electromagnetics, ultrasonics, X-rays or electrostatics. However, if this information could be included in some numerical technique a more accurate identity for the flaw could be produced. This could then be used to predict the life expectancy of the component and eliminate unnecessary wastage. One method of using this information is known as an inverse analysis technique. Inverse analysis is the general name given to problems where not all of the information, that is usually required to analyse the system, is known a priori. For example, a specimen subjected to a known load may have unknown material parameters or the load required to produce a known distortion may be unknown. In these problems the prescribed boundary conditions are insufficient to solve the problem using standard analysis techniques and an iterative solution is usually required. An initial guess is made for the unknown parameters which are subsequently altered to achieve 'better' results. The flaw identification method is a type of shape identification problem. In these shape identification problems the boundary of the specimen is usually required to satisfy some extra conditions, usually to improve the efficiency of the component. Problems of this type are often used in the design of structures and examples of this type of problem are given in Ref. [1]. In the identification of a flaw, the unknown boundary is an internal surface and is chosen so as to minimize the error in responses at the sensor points. Some research papers have been presented for suclh a problem, for example Friedman and Vogelius [2] have shown that a single crack can be uniquely identified from measurement of two boundary conditions. Later their work was extended to many cracks by Bryan and Vogelius [3]. A review of BEM methods for #To whom all correspondence should be addressed. E$34/4-F
453
454
s . c . MELLINGS and M. H. ALIABADI
flaw identification problems can be found in [4]. Various systems can be used to define this error, for instance, the specimen could be subjected to a static force and the displacements or stresses measured at various points. Then, the error between these measurements and those expected with a known crack position can be computed. Other systems can be used, for instance a dynamic force or a potential field could be applied. The values at the sensor points that are expected from known crack positions need to be computed numerically and this can be done with a variety of analysis methods. For instance, domain methods such as the Finite Difference or Finite Element Methods can be used. Here a new formulation of the Boundary Element Method, known as the Dual Boundary Element Method ( D B E M ) will be used. In the flaw identification method, the restriction of the discretization to the boundary is beneficial, since the crack shape must be changed at each optimization step. Use of this type of domain method would make it necessary to involve a remeshing routine after each step in the procedure. In the method developed, however, this type of remeshing is only required after the iteration procedure is complete and then only if a more accurate result is required. For the majority of the flaw identification methods previously developed with the Boundary Element Method, it has been necessary to divide the domain along the line of any crack. In this work, the new D B E M formulation is used to eliminate the need for this sub-division. The flaw identification method presented here has previously been formulated for two-dimensional analysis with both potential and elastostatic systems, [6, 7], but in this paper the method is extended to three-dimensional potential systems. A similar formulation has also been used by Bryan and Vogelius [8] for two-dimensional problems. In their formulation they transformed a perfectly insulating crack into a perfectly conducting crack by utilizing the conjugate harmonic potential.
2. F L A W I D E N T I F I C A T I O N
METHOD
It is assumed that measurements of either potential or flux values are available at a finite number of points in the cracked specimen. These points will be referred to as sensor points or nodes and may be either boundary or internal points. It is necessary to create a model for the crack shape, using design variables. Many models may be used, however too many design variables will cause the analysis to be slow. This will also increase the chances of finding false solutions, due to the presence of more local minima in the error function. Hence initial designs are chosen for their simplicity, though these may be refined later when an approximation to the crack has been found. Using the chosen design variables, an initial crack shape is selected and using Dual Boundary Element analysis, with the same boundary conditions as for the experimental set-up, the potential or flux values are computed at the sensor points. Unless the design for the initial flaw exactly matches the actual flaw, these two sets of readings at the sensor points, will differ. The size of these differences indicates the size of error in the guessed crack shape and minimization of this error should minimize the error in the crack identity. The error at the sensor points in the specimen can be quantified into an error function, F(Z), defined by V Z .N~ ~ (Au"(Z)) 2 F(Z) ro (1) The term A u ' ( Z ) is the error in the response at the nth sensor, Ns is the number of sensor points and Z is the vector of design variables. The term F ° is a value used in scaling the error function, the value of which is taken to make the error for the initial flaw to be unity. This is a nonlinear function, constrained so the crack remains inside the body. In the work presented here, only internal cracks are to be identified and, in this case, the use of constrained
3D flaw identification
455
analysis can be avoided by scaling of the changes made to the design variable vector. Hence a nonlinear, constrained optimization method can be used in the identification method and since evaluation of the error function is computationally expensive a first order optimization method will be used. In these first order methods, it is necessary to evaluate the gradient of the error function, VF, defined by
VF(Z)=(OF(Z)~
(2)
\ OZm ]
where Zm is the mth component of the design vector. In order to evaluate this function, derivatives of the potential and flux values at the sensor points are required. These values are known as design sensitivities and in the method presented here they are to be evaluated with the derivative DBEM. A selection of :first order optimization methods are available to minimize this error function, two of which are the Steepest Descent and BFGS Quasi-Newton methods. The latter of these is used in most of the examples presented here, due to its superior convergence rates and greater reliability. For details of these methods see [9, 10].
3. DUAL B O U N D A R Y ELEMENT METHOD The boundary element method is based on an integral equation transformation of the governing partial differential equations. In potential work, the governing equation is the Laplace equation V2u(x) = 0 for x in ~ (3) where f~ is the domain of the body and u(x) is the potential function. Boundary conditions are prescribed on the boundary of the body F and these can be written as u(x) = if(x)
for x on F 1
q(x) = q(x)
for x on F 2
(4)
where F 1 and F 2 are portions of the boundary F such that F = F 1 + F 2. The flux function, q(x), is the normal derivative of the potential function and tT(x) and ~(x) are the prescribed potential and flux values. From equation (3), use of a weighted residual technique gives an integral equation u(X t) +
fr
q*(X', x)u(x)
dr(x) = fr u*(X', x)q(x) dr(x).
(5)
In this equation X' is an internal source point, x is a boundary point and the terms u* and q* are known as the fundamental solutions, given for three-dimensional analysis by
u*(X', x) = q*(X', x)
1 4a'r 0u(x) On
1
Or
4 ~ 2 On
where r = Ix - X'I. This identity (5;) is valid for all points in the body, f~, including the boundary, F. Taking the internal point, X', to the boundary produces the boundary integral equation, c(x')u(x') + fr q*(x', x)u(x) dF(x) = fr u*(x', x)q(x) dF(x)
(6)
where c(x') is a jump term arising from taking the first integral to the boundary. For the
456
S.C.
MELLINGS
a n d M. H . A L I A B A D I
derivation of these equations, see Brebbia and Dominguez [11]. In order to evaluate the unknown potential and flux terms, it is necessary to discretize the boundary F into elements. In an element, the geometry as well as the potential and flux functions are assumed to have known variations. In this work, these variations are assumed to be quadratic and 8-noded quadrilateral elements are used. Transformation of the geometry in these elements to a local twodimensional system allows all the terms, in the integral equation, to be transformed into the local coordinate system. These terms can then be evaluated with the use of a standard Gaussian integration technique. c(x')u(x') + ~
3,=1 n = l
u 3"~; '
I
f_llq*(x ' , x(~, n))M'(~, n)J3"(x(~:, n)) dE dn =
q3" 3,=1 n = l
u*(x', x(~, r/))M"(~, r/)J3"(x(~¢, r/)) d~¢ dr/. 1
(7)
1
The terms M'(~, ~7) are the quadrilateral shape functions at the nth local node, J3"(x(~, 7/)) is the Jacobian of transformation and N is the number of boundary elements. In the Boundary Element Method, this equation is used at each node on the boundary, forming a set of linear equations in terms of the unknown potential or flux results at each node. This set of equations can be solved using a standard method, thus giving the unknown values.
3.1 Cracked bodies The standard application of the potential equation (7) to a body containing a crack does however cause difficulties. If equation (7) was used on both crack surfaces, two identical equations would be formed and the resulting matrix equation would be singular. The usual way of overcoming this problem is to sub-divide the body into two domains along the line of a crack. This allows each sub-region of the body to be computed independently and since the integration paths are now different, the system of equations is no longer singular. However, this introduces extra nodes along the division line, thus increasing the size of the system of linear equations. A new formulation, known as the dual Boundary Element Method, has been developed to overcome this problem without the use of sub-division. On one crack surface the standard, potential equation (6) is to be replaced by a new independent flux equation. This new equation is computed from the normal derivative of the potential equation. Using equation (5) for the potential value at an internal point X', the potential derivative at X', with respect to the coordinate direction, i, can be evaluated to be u,i(X') +
frq*(X',x)u(x ) dF(x)=
£ u*(X', x)q(x)dF(x)
(8)
where
1
u*(X', x) = 4zr----5 r,i
q*(X',x)=
1 (30r -n~).
4/t'r 3
c3-nr'i
It is necessary, as for the external boundary to evaluate the equations on a crack. An extra term is now produced due to the singularity on the opposite crack face. For the potential equation on a smooth crack boundary, the equation becomes
,1 u(x ,) + 1 u(x ") + f, q*(x', x)u(x) dr(x)
=
u*(x',x)q(x) dr(x)
(9)
3D flawidentification
457
where x" is the coincident point on the opposite crack surface. Similarly for the derivative potential equation on a smooth crack t!
1 u i(x')+ ~u i(x ) + fr q*(x', x ) u ( x ) d r ( x ) = fr u*(x', x)q(x)dr(x) 2
(lO)
where f represents the Cauchy principle value integral and ~ represents the Hadamard principle value integral. Since q(x)= u.i(x')ni(x'), multiplication of (10) by the normals at the source node will ~{ive the flux equation on a crack to be q(x') - ~ q(x") + ni(x')
q*(x', x)u(x) dr(x) = ni(x')
u*,i(xt, x)q(x) dr(x).
(11)
For a flux free crack q(x') = q(x") = 0. Discretization of equations (9) and (11) into boundary elements gives 1
1 ,, u(x , ) + 5"(x )+ ~
u "n
~_~J'+'
q*(x', x(~:, r/))M~(~, r/)Jr(x(sc, r/)) d~ dr/
3,=1 n = l
1
=
--1
u rn
u*(x', x(~, r/))M"(~, r/)J'(x(~, r/)) d~ dr/
"/=1 n = l
~q(x)-
1
q(x")+ni(x') ~ ~u~'n~ V=I n = l
q*(x',x(~,r/))Mn(~:,r/)JV(x(~,r/))d~dr/
a-1
= n;(x') ~= =mn=l qVn
(12)
1
1
f,;, 1
i u*(x', x(~:, r/))Mn(/.?, r/)JV(x(~, r/)) d~: dr/.
(13)
The Dual Boundary Element Method comprises of the use of both the potential equation, as given by equations (7) and (12), and the flux equation, as in equation (13). The potential equation (7) is to be used on the outer boundary of the problem where no problems arise due to coincident points. On one crack face the potential equation (12) is used, whilst the flux equation (13) is used on the other crack face. This gives, on application of the boundary conditions, a system of linear equations that may be solved for the unknown potential or flux values. This can be written in the matrix form as AX = F
(14)
where the vector F is formed from the application of the prescribed boundary potential and flux values, X is the vector of unknown potential or flux values and A is a coefficient matrix computed from the remaining terms. In the flaw identification method the Dual Boundary Element Method is used to compute the sensor values at tlhe boundary points, and if values are required at internal points, equations (5) and (8) are used to compute the values.
4. DESIGN SENSITIVITIES As mentioned previously, it is necessary in the minimization of the error function, to evaluate the design sensitivities. These are the sensitivities of the potential and flux values to changes to the design variables and are found by taking the derivatives, with respect these changes. There are three methods of computing these sensitivities; the finite difference approximation method, the material derivative method and the implicit, or direct, differentiation method.
458
S.C. MELLINGS and M. H. ALIABADI
In the first of these methods, values are computed using a finite difference approximation to a derivatives, for example a forward difference approximation is given by 0f(Z)
3Zm
-
f(Z
"~
~Zm) - f ( Z ) ~Zm
(15)
where 8z is a small forward difference step. This method gives an approximation that is simple to evaluate, but is computationally expensive and the accuracy is often very poor. The material derivative method is a variational method, a description of which can be found in Refs [14] and [15]. The third method mentioned above is the one used here and involves the direct differentiation of the boundary element equation(s) with respect to the design variable zm. The equations used in the Dual Boundary Element Method (7), (12) and (13) give implicit relationships for the potential and flux nodal values. Direct differentiation of these discretized equations would give new equations in terms of the potential, flux, potential derivative and flux derivative nodal values. Differentiation of the potential equation, on a standard boundary (7), gives c,m(x')u(x') + c(x')u,m(x') +
u,V~ ~'=1 n = l
+
+
J--I
J--I
u 3"n 3,=1 n = l
=
q (x, x(~:, r/))J,Vm(X(e, r/))m"(sc, r/) dgdr/
J--1 J--I
q,~
u*(x', x(sc, r/))M~(~, r/)J3"(x(sc, 7/)) d¢ dr/
y=l n=l
1
+
1
u*m(x', x(~, r/))J3"(x(~:, r/))Mn(~, 77) dsc dr/
q3" 3"=1 n = l
+
1
q,m(X, X(¢, r/))J3'(x(~, r/))Mn(~, 7-1)dg dr~
U3"n 3"=1 n = l
q*(x', x(~, r/))M"(sc, r/)J3"(x(~:, 7/)) d~ dr/ 1
1 ~ 1
q3"n 3"=1 n = l
u ( x , x(¢, r/))J3",m(X(g, r/))mn(¢, 7/) de dr/. J-1
(16)
J--1
On a smooth crack the derivative of the potential dual boundary element equation (12) gives
x)
1 , 1 ,,) + ~.~ - - Uyn ~U,m(X )+~U,m(X 3"=1 n = l
+
u rn 3"=1 n = l
+
f['f-' f['f[" 1
1
1
1
~ u 3"n 3"=1 n = l
=
q.~ v=l n=l
+ +
1
"-1
.!-I
q*(x',
X(¢, r/))M"(¢,
r/)J3"(x(e,
r/))dCd,
1
q*,,(x', x(~, r/))J3"(x(~, r/))Mn(sc, r/) d~ dr/
q*(x', x(sc, r/))JV,,(x(g, r/))Mn(g, r/) d~: dr/
u*m(x', x(¢, r/))J3"(x(¢, r/))m"(g, r/)de dr/
q3"n y=l n=l
1
u*(x', x(sc, r/))Mn(~:, r/)J3"(x(~:, 77)) d~: dr/ 1
q3"n 3"=1 n = l
f-l f_.
u (x, 1
-I
71)
dr/.
(17)
3D flaw identification
459
Finally, the flux Dual Boundary Element equation (13) differentiates to give
u z"
~q,.,(x )-~q~m(X")+ ni,.,(x') 3,=1 n = l
-[- ni(X t) ]~_~
I'~Y,~
~t=l n = l
+ ni(x') [~
1
f_~' q*(x', x(E, ~?))M"(E, ~7)JV(x(E, ~/)) dE d~? 1
u v"
"r=l n=l
q*(x', x(E, ~)M"(E, rt)JV(x(E, r/)) dE d*/ 1
1
, q,~m(X', X(E, n))JV(x(E, r/))M"(E, 77) dE dT/ 1 ~l
+ ni(x') /..~ "•
U Vn
T=I n=l
n=l
+ ni(x' ) .
1
+ ni(x') , , 3,=1 n = l
1
1
1
1
1
1
u*.,(x', x(E, n))JV(x(E, rl))Mn(E, "q) dE drl
qV-
-r=l n=l
+ n/(x")
1
f_~' u*(x', x(E, 7/))M"(E, ~?)JV(x(E, 7/))dE dr/
q,~
"),=1 n = l
f_, q*(x', x(E, r/))J,~m(x(E, ~))M"(E, ~l) dE d~? f - ' u*(x', x(E, ~))M"(E, r/)JV(x(~:, n))dE d~/
qV.
= ni, m(X' ) : ~ y=l
1
qV-
u*(x', x(E, r/))J,~.,(x(E, 7/))M"(E, r/) dE dr/.
(18)
In these equations the derivatives of the fundamental solutions u*, q*, u* and q*m are given by 1
u*,.(x', x) -
(19)
4xr2 r,m
1 (3Or ) q*,.(x', x) =4zrr ~ ~~nr, m -- [rlnl],m U*m(X',X) -q*.,(x',x)= 4zr ~ 1 (3Or ~n
(20)
1 4 ~ i(3rirm -- ro. )
(21) )
[5F'ir'm - ri'm] - 3[rtnt]'mri - 3nirm + ni'mr "
(22)
The potential and flux values can be computed using the DBEM as described earlier. Hence these equations implicitly relate the potential and flux derivatives. Derivative boundary conditions can be found from differentiation of the boundary conditions prescribed earlier, since it is assumed that they are independent of the crack shape and position. Hence the derivative boundary conditions are given by u,m(x) = 0
x E F~
q,m(x) = 0
x ~ F 2.
(23)
These derivative equations can now be used to provide a system of equations for each design variable. Use of the boundary conditions and the computed boundary potential and flux values allow the system to be written in matrix form as A X , , . = F,m - A . , , X .
(24)
In this system, X,,. is the vector of unknown flux and potential derivatives, (F,,,, - A,.,X) is the vector formed by the application of boundary conditions and A is the matrix of coefficients identical to that used in the DBEM analysis. In the flaw identification method, the system of
460
S.C.
MELLINGS
a n d M. H . A L I A B A D I
equations is solved with the LU Decomposition method, in which the A matrix is decomposed into an upper and a lower triangular matrix. The system is then solved using forward and backward substituion. Since the A matrix used here is identical to that in the previous matrix system (14) the decomposition can be preserved, rather than the matrix itself, and this can then be reused for each of the derivative analyses.
5. H Y P E R - S I N G U L A R I N T E G R A T I O N In the implementation of both the DBEM and the derivative formulation it is necessary to take special care when r ~ 0 (i.e. x ~ x ' ) , since the fundamental solutions q*, u*, q* and u'~,t and their derivatives all contain terms of the form O(1/r"), where n > 0. In the usual boundary element method, the rigid body condition is used to compute the weakly singular integrals that occur when the source and field points coincide. In the DBEM, however, this technique cannot be used on the crack, due to the presence of two equal and opposite singular terms, from each crack face, which cancel in the summation of row terms. Hence, some form of hyper-singular integration is required. The crack is assumed to be traction free, so all of the terms on the right-hand side of the equations (7)-(13) and (16)-(18) disappear. In the DBEM equations, the terms remaining on the left-hand side of the equation are given by the weakly singular kernel q* and the hyper-singular kernal q,*, i.e. -I
I 1=
f_ f_
q*(x', X(¢, r/))M"(~:, r/).D'(x(~:, 7/)) de dr/
(25)
, q*(x', x(~, r/))M"(s~, r/)Jr(x(~, r/)) d~ dr/.
(26)
1
~l
12 =
In the derivative equations, two extra integrals are required, these being the derivatives of I l and 12 given by '-I
J' =
f_ f_
[q*(x', x(~:, r/))JY(s¢, 7/) + q*(x', x(sc, r/))J~(x(¢, r/))]M(¢, r/) d~ dr/
(27)
[q*,.(x', x(~, r/))J~(~, r/) + q*(x', x(~, ~/))J,~,.(x(~, r/))]M(~, ~1) d~ dr/.
(28)
1
~l
j2 =
1
A transformation is made to a local polar coordinate system centred on the singular point x'. This allows the singularity to be removed by the exclusion of a small region of radius e, which is reduced by taking the limit as e ~ 0. I ~= lim e~O JO
12 = l i m c] TM ri~(o) e---*0 do
F~(p, O) dp dO
(29)
F2(p, 0) dp dO.
(30)
Jo~(e,O)
-'a(e,O)
In these expressions a(p, O) is a function representing the exclusion around the singularity and ~(0) the function for the external contour of Fv in polar coordinates, see Mikhilin [13]. The integrands are then given by
F'(p, e)=q*(p, e)M"(p, e)J~(p, e)p
(31)
0)= q~(p, O)M'(p, e)J='(p, e)p.
(32)
F2(p,
3D flawidentification
461
Similarly for the derivative integrals, the integrals can be transformed to polar coordinates giving
G~(p, 0) dp dO
(33)
G2(p, O) dp dO
(34)
G~(p, 0) = [q* M"J "/+ q*M"J~,.]p
(35)
G2(p, 0) = ['~*. L"I , t m M " J
(36)
jt = lim e-*O "0
" a(e,O)
j2 = lim e--*0 J 0
-'~(e.O)
with the integrands given by
~ + t'l ~ *,i M " J v, m J1^ l ~"
It is possible to rewrite the functions F 1, F 2, G I and G 2 in terms of a singular term, that can be integrated analytically, and a non-singular term, which can be computed using standard Gaussian integration, using ( F~(p, 0 ) = Fl(p, e)
FL~ F ~)_ ~ ( 0 ) 0 ) + - P
(37)
462
S.C. MELLINGS and M. H. ALIABADI
However, since the components of the normal can be written as ni = gdJ then Jn~ = g , hence this is evaluated in the expansion:
Ogk gk(p'O)=gk(rl)+ 0 (0¢;1 ~ cos O+-d-'~zslnO)+O(o2)" = gO + pg~(O) + 0(02).
(43)
Similarly the expansion of the shape functions can be written as
{OM" OM" O) O(02 ) M"(p,O)=Mn(n)+O o~ c°s 0 + - ~ - 2 sin + = M~ + pMT(O) + 0(02).
(44)
Manipulation of these terms, see Guggiani [12], allows the terms F 1 and F 2 to be expanded, giving Fl(p, 0) -
1 Or Mnj~ p 4zr On
- 4zl [F~-;(0)+O(1)] F2(p, 0 ) =
(45)
1 [3 Or J ~ r i - n i J ' ] M ' p 4~r 3 k On " 1 [ FZ_z(0)
- 4x[--~+
F~I(0)
]
0 + O(1)j.
(46)
The terms in the derivative integrals can be found in a similar way. The derivative of the distance components, ri,m is given by
rg,,.(p, O) = pAi,m(O ) -t- p2Bi,m(O ) -t- O(p 3)
(47)
where
Ai, m(O ) = Xi, lm COS 0 -b Xi,2m sin 0 Bi.,.(O) = Xi, ll,. cos 2 0 + 2X;,12,. cos 0 sin 0 + X~,22.. sin 2 0. This gives the derivative distance components to be
r m(p, O) = OAm(O) + 02Bm(O) + O(03)
(48)
where
Am(O) --
A,(0). A,.,,(0) A(O)
B,,( O)
ak,.,(O)Bk(O) + A~,(O)Bt,,.,(O) A(O)
Ak(O)Bk(O) A,.,(O)" "A'O ~, ))
The derivative of the Jacobian components is also required, with gk,,,, = J,,,,nt, + Jn~,.,,,:
{ Ogk,,. sin 0) + O(p z) gk,,,(O, O) = gk,m(rl) + 0~, 01~, COS 0 + Ogk.m 0~2 = gk,,O + Pgkm,(O) + O(02).
(49)
3D flawidentification
463
These terms can then be added in to give the integrands required, i.e. Gl(p, 0 ) = ~ ( 3 0 r r 4 7 0 .3
On ,m
-[rtnll,,,)M"JVp
- 4zl[G[p(0)+O(1)] G2(p, 0 ) =
1 (3 ~Or [5r, ir, m 4rtr------~
- ri,m] - 3 [ r l n t ] , m r , i -
(50)
3nir m + ni,,,r)M"JVp
1 [ G2_2(0) G2_I(0) ] =-4---~1_ p2 ~ - - p ~-0(1) .
(51)
6. RESULTS In the method presented here, it is necessary to examine both the accuracy of the design sensitivities as well as the effectiveness of the optimization procedure used. In the case of the sensitivity values., a comparison can be made with a series of finite difference approximations which will give an idea of the accuracy of the method. The performance of the flaw identification car, be tested using a known shape for the actual crack. The responses, which are usually required experimentally, can then be computed with the Dual Boundary Element Method with the known crack identity. The flaw found in the identification procedure can then be compared to the known flaw configuration.
6.1 Sensitivity values In this comparison, it is assumed that the cylinder contains a fiat, horizontal, penny shaped crack of unknown radius and location. Four design variables are used, these being the radius of the crack, R, and the three coordinate positions of the centre of the crack, i.e. x, y and z. Comparisons are to be made with the approximations given by the forward difference steps 8z,, = 0.1, 0.01, 0.001 and 0.0001. These approximations allow the percentage errors at each node, xn, of the ,cylinder mesh to be defined by
s(x') -
u f w d t x n ) -- U id ( x n ~ ,m \ ), ,m,~ ' * 1 0 0 % . ufWdtx~ ,it'/ k ]
(52)
In this expression, u f w d ( x n ) represents the forward difference approximation computed with the approximation formula given by equation (15) and uid(x ") is the value computed by the direct differentiation method, presented in this paper. A percentage error value is computed from the error at all N nodes as
E
VX e2(x") N
(53)
The value of this error function E is plotted on the graph for various sizes of the forward difference step, iSz,, and for each of the four design variables, x, y, z and R. As can be seen, in Fig. 1, the approximations were accurate and improved as the forward difference steps became smaller.
464
S.C. MELLINGS and M. H. ALIABADI i
4
i
i
i
i
i
ii
I
i
i
i
i
i
zxXC~ordin~t~o~~rackC~ntr~
ir
i
i
i
'~X~
V Y Coordinate of Crack Centre
/
z coo ,u..t ot c .ok C,r t ,
o
i
r
i
i
ii
/ ~
/
/
\
[] Radius of Cr
>
3
0
QJ
2
cO
0
1
r~
0
10-4
I
I
2
3
I
I
~
t
i irk~
t
t
10-a
2
S
Logarithmic
I
i
I
i
, ,
~
~ 2
10-s
Forward
Difference
3
10-t
Step Size
Fig. 1. Comparison of design sensitivities with forward difference approximations.
6.2 Inverse analysis 6.2.1 Identification of a penny shaped crack. Consider a cylinder of height 2 units and radius 1 unit, with a temperature of +100 units at the upper end and - 1 0 0 units at the lower end, with all the other sides being flux free. The cylinder actually contains a fiat, circular crack, of radius R = 0.5 with the crack centre at the origin (x = y = z = 0), as in Fig. 2. In the identification procedure, it is assumed that the crack is circular and fiat. The initial guesses made for the parameters of this crack are for a radius of R --- 0.1 and the central coordinates of x = 0.25, y =0.25, z =0.5. Potential measurements are taken at boundary sensor points which are positioned at the boundary nodes of the discretized cylinder where flux is zero. In the identification of this crack the BFGS Quasi-Newton Method was used to minimize the
u=+lO0
/
Y
/
q=O
/ /
X
Actual Initial
Crack j
u=-lO0 Fig. 2. Geometry of initial and actual cracks.
Crack
3D flaw identification
465
8th
Final Crack Initial Crack
lOth I t e r a t i o 5th Iteration
Fig. 3. Crack identities during iteration procedure.
error. Figure 3 shows the initial and final crack identities as well as those at the 4th and 7th iterations. The convergence both of the design variables (Fig. 4) and the sensor error normalized with respect to the initial value (Fig. 5) show that the crack can be identified successfully using only 20 iterations.
0.5,~ 0.4
l>
0.3
/
I V Y Coordinate of Crack Centre I
#
I 0 Z Coordinateof Crack Centre I
r=.=4
c~ a~ 0.2
t> la0
.,=.4
~D O.t
0.0
--0.|
I
o
I
z
I
I
4
I
I
e
I
l
a
I
I
1o
I
I
I
m
I
14
I
I
te
I
l
I
la
Iteration Number Fig. 4. Convergence of design variables in the identification of a penny shaped crack.
zo
466
S.C. MELLINGS and M. H. ALIABADI 100
•
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
0 ,--] tO-t 0 I-i I..i "~
e lO-e
of--I
a lO-a 0 Z 2
10-4 0
~
4
6
8
10
12
14
16
18
20
ILeration N u m b e r Fig. 5. Convergence of the normalized error in the identification of a penny shaped crack. F i g u r e 6 shows a c o m p a r i s o n b e t w e e n the f o u r o p t i m i z a t i o n m e t h o d s m e n i t i o n e d earlier. It can b e seen, that, a l t h o u g h the r a t e of c o n v e r g e n c e for the s i m p l e r m e t h o d s is l o w e r , all t h e p r o c e d u r e s u s e d i d e n t i f i e d the c r a c k a c c u r a t e l y . 6.2.2 Identification o f an elliptical crack. In the s e c o n d e x a m p l e , t h e a c t u a l c r a c k is the s a m e as in t h e p r e v i o u s e x a m p l e , b u t n o w t h e initial c r a c k is a s s u m e d to b e elliptical, with t h e c e n t r e of the elliptical c r a c k (xo, Yo, zo) is to b e i d e n t i f i e d a l o n g with t h e elliptical p a r a m e t e r s a a n d b defined by (X-Xo) 2 (y-yo) 2 + 1 and z = Zo. (54) b2 02
100
- q÷eer~es~ ih~scent Method
F
I
•
Conjugate
Gradient
Method
1 0-4
0
4
8
12
16
20
24
28
32
56
Iterotion Number Fig. 6. Comparison of the convergence of the normalized error for various optimization methods.
3D flaw identification
467
z u=+lO0
q=O
,21/~'--I
~
.....
Initial
I\Actual
Crack
Crack
u=-100 Fig. 7. Initial and final geometry for the identification of an elliptical crack.
The initial design is chosen with the design variables initially being a = 0.4,
b = 0.1,
Xo = -0.25,
Yo-- 0.25
and
Zo = -0.5
and flux values are measured at the sensor points. In this example the sensor points were located at the two ends of the cylinder where temperature was prescribed (see Fig. 7).
Initial Crack Final Crack
3rd Iteration 8th
\
Fig. 8. Crack identities during the identification of an elliptical crack.
468
S.C. MELLINGS and M. H. ALIABADI 0.5 0.4 0.3 0.2 ~"
0.1
~
0.0 -0.!
~
-0.2
l" El]~Parameter-b[
•~ -o.3 1 / I~
-0.4 L
]
? CentraICoordinate-x ~, C e n t r a l Coordinate - y I * Central Coordinate - z
-o.5 [ / -0.6
-0.7
0
4
8
12
16
20
24
28
32
[
J I
36
4
Iteration N u m b e r Fig. 9. Convergence of design variables in the identificationof an ellipticalcrack.
For this problem, the crack is again identified successfully, with the convergence of the design variables and the error normalized with respect to the initial value are shown in Figs 9 and 10. Figure 8 shows the final position of the crack as given by the interation procedure, along with the initial and actual crack identities.
7. CONCLUSION In this paper, a new method has been presented for the identification of cracks inside three-dimensional structures. The method is found to be highly efficient with the reuse of the
10
A
0
,
~
,
I
I
,
,
,
,
,
,
,
,
I
I
,
2
0 ,=1 t~ o
lO-t
It.1 8
lO-e
0 Z I
10-~
I
0
4
I
s
I
i
tz
I
i
t6
I
I
zo
I
I
z4
I
i
~
I
~
I
I
se
Iteration N u m b e r Fig. 10. Convergence of the normalized error in the identificationof an ellipticalcrack.
3D flaw identification
469
DBEM system matrix for each of the derivative analyses. The derivatives computed by this method are found to be highly accurate and the method is successful in the identification of an internal crack.
REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
G. N. VANDERPLAATS, A I A A J. 20, 992 (1981). A. FRIEDMAN and M. VOGELIUS, Indiana Univ. Math. J. 38, 527 (1989). K. BRYAN and M. VODELIUS, SlAM J. Math. Anal. 23, 950 (1992). M. TANAKA and Y. MASUDA, Engng Anal 3, 138. N. NISHUIMURA and S. KOBAYASHI, Int. J. Numer. Meth. Engng 7, 59 (1990). S. C. MELLINGS and M. H. ALIABADI, Engng Anal, Boundary Elements 12, 275 (1993). S. C. MELLINGS and M. H. AL1ABADI, Flaw identification using the boundary element method. To appear in hzt. J. Numer. Meth. Engng. K. BRYAN and M. VOGELIUS, Int. J. Engng Sci. 32, 579 (1994). R. FLETCHER, A review of methods for unconstrained optimization. In Optimization (Edited by R. FLETCHER), Academic Press, New York (1969). D. M. GREIG, Optimisation. Longman, London (1980). C. A. BREBBIA and J. DOMINQUEZ, Boundary Elements--An Introductory Course. Computational Mechanics Publications, McGraw-Hill, New York (1989). M. GUIGGIANI, G. KRISHNASAMY, T. J. RUDOLPHI and F. J. RIZZO, A general algorithm for the numerical solution of hypersingular boundary element integral equations. To appear in A S M E J. Appl. Mech. S. G. MIKHILIN, Multidimensional Singular Integral Equations. Pergamon Press, London (1965). A. LONGO, J. UNZUETA, E. SCAEIDT, A. A L V A R E Z and J. J. ANZA, Adv. Engng Software 16, 135 (1993). K. K. CHOI and E. J. HAUG, J. Structural Mech. U , 231 (1983). A. J. KASSAB, F. A. MOSLEHY and A. B. DARYAPURKAR, Nondestructive detection of cavities by and inverse elastostatics boundary element method. To appear in Engng Analysis Boundary Elements. (Received 24 October 1994; accepted 22 November 1995)
£S 3414-G
Pergamon, 0020-7225(95)00125-5
Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0020-7225/96 $15.00+ 0.00
THREE-DIMENSIONAL FLAW IDENTIFICATION INVERSE ANALYSIS
USING
S. C. MELLINGS and M. H. A L I A B A D I t Wessex Institute of Technology, Ashurst Lodge, Ashurst, Southampton SO40 7AA, U.K. Abstract--In this paper, a three-dimensional crack identification method is developed. The Dual Boundary Element Method (DBEM) is used to compute potential and flux values for known crack identities and first order optimization is used to minimize the error due to approximations made to the true crack identity. Design sensitivities are computed with a derivative formulation of the Dual Boundary Element Method. The accuracy of these sensitivity values is investigated by a comparison with values obtained using a series of finite difference approximations. The flaw identification method developed is then demonstrated for the identification of a simple internal crack.
1. INTRODUCTION In maintenance of structures, components are often tested in order to detect any flaws. The techniques used can give an approximate size of the flaw and if the flaw is considered too large, the component is usually discarded. However, this procedure causes some components which do not present a risk of failure to be unnecessarily wasted. On the other hand, some smaller cracks may present a greater risk of failure due to their position. If a technique could be employed to identify the actual flaw shape and size, then it may be possible to make a more reliable predictio:a for the lifetime of the component. The methods used to detect internal flaws are known as Non-Destructive Testing (NDT) methods. In most of these methods, a field is applied to the component and responses at various sensor points are measured. The presence of the flaw will cause a distortion in the field which can be measured and used to approximate the flaw size. Common fields used in such NDT methods are electromagnetics, ultrasonics, X-rays or electrostatics. However, if this information could be included in some numerical technique a more accurate identity for the flaw could be produced. This could then be used to predict the life expectancy of the component and eliminate unnecessary wastage. One method of using this information is known as an inverse analysis technique. Inverse analysis is the general name given to problems where not all of the information, that is usually required to analyse the system, is known a priori. For example, a specimen subjected to a known load may have unknown material parameters or the load required to produce a known distortion may be unknown. In these problems the prescribed boundary conditions are insufficient to solve the problem using standard analysis techniques and an iterative solution is usually required. An initial guess is made for the unknown parameters which are subsequently altered to achieve 'better' results. The flaw identification method is a type of shape identification problem. In these shape identification problems the boundary of the specimen is usually required to satisfy some extra conditions, usually to improve the efficiency of the component. Problems of this type are often used in the design of structures and examples of this type of problem are given in Ref. [1]. In the identification of a flaw, the unknown boundary is an internal surface and is chosen so as to minimize the error in responses at the sensor points. Some research papers have been presented for suclh a problem, for example Friedman and Vogelius [2] have shown that a single crack can be uniquely identified from measurement of two boundary conditions. Later their work was extended to many cracks by Bryan and Vogelius [3]. A review of BEM methods for #To whom all correspondence should be addressed. E$34/4-F
453
454
s . c . MELLINGS and M. H. ALIABADI
flaw identification problems can be found in [4]. Various systems can be used to define this error, for instance, the specimen could be subjected to a static force and the displacements or stresses measured at various points. Then, the error between these measurements and those expected with a known crack position can be computed. Other systems can be used, for instance a dynamic force or a potential field could be applied. The values at the sensor points that are expected from known crack positions need to be computed numerically and this can be done with a variety of analysis methods. For instance, domain methods such as the Finite Difference or Finite Element Methods can be used. Here a new formulation of the Boundary Element Method, known as the Dual Boundary Element Method ( D B E M ) will be used. In the flaw identification method, the restriction of the discretization to the boundary is beneficial, since the crack shape must be changed at each optimization step. Use of this type of domain method would make it necessary to involve a remeshing routine after each step in the procedure. In the method developed, however, this type of remeshing is only required after the iteration procedure is complete and then only if a more accurate result is required. For the majority of the flaw identification methods previously developed with the Boundary Element Method, it has been necessary to divide the domain along the line of any crack. In this work, the new D B E M formulation is used to eliminate the need for this sub-division. The flaw identification method presented here has previously been formulated for two-dimensional analysis with both potential and elastostatic systems, [6, 7], but in this paper the method is extended to three-dimensional potential systems. A similar formulation has also been used by Bryan and Vogelius [8] for two-dimensional problems. In their formulation they transformed a perfectly insulating crack into a perfectly conducting crack by utilizing the conjugate harmonic potential.
2. F L A W I D E N T I F I C A T I O N
METHOD
It is assumed that measurements of either potential or flux values are available at a finite number of points in the cracked specimen. These points will be referred to as sensor points or nodes and may be either boundary or internal points. It is necessary to create a model for the crack shape, using design variables. Many models may be used, however too many design variables will cause the analysis to be slow. This will also increase the chances of finding false solutions, due to the presence of more local minima in the error function. Hence initial designs are chosen for their simplicity, though these may be refined later when an approximation to the crack has been found. Using the chosen design variables, an initial crack shape is selected and using Dual Boundary Element analysis, with the same boundary conditions as for the experimental set-up, the potential or flux values are computed at the sensor points. Unless the design for the initial flaw exactly matches the actual flaw, these two sets of readings at the sensor points, will differ. The size of these differences indicates the size of error in the guessed crack shape and minimization of this error should minimize the error in the crack identity. The error at the sensor points in the specimen can be quantified into an error function, F(Z), defined by V Z .N~ ~ (Au"(Z)) 2 F(Z) ro (1) The term A u ' ( Z ) is the error in the response at the nth sensor, Ns is the number of sensor points and Z is the vector of design variables. The term F ° is a value used in scaling the error function, the value of which is taken to make the error for the initial flaw to be unity. This is a nonlinear function, constrained so the crack remains inside the body. In the work presented here, only internal cracks are to be identified and, in this case, the use of constrained
3D flaw identification
455
analysis can be avoided by scaling of the changes made to the design variable vector. Hence a nonlinear, constrained optimization method can be used in the identification method and since evaluation of the error function is computationally expensive a first order optimization method will be used. In these first order methods, it is necessary to evaluate the gradient of the error function, VF, defined by
VF(Z)=(OF(Z)~
(2)
\ OZm ]
where Zm is the mth component of the design vector. In order to evaluate this function, derivatives of the potential and flux values at the sensor points are required. These values are known as design sensitivities and in the method presented here they are to be evaluated with the derivative DBEM. A selection of :first order optimization methods are available to minimize this error function, two of which are the Steepest Descent and BFGS Quasi-Newton methods. The latter of these is used in most of the examples presented here, due to its superior convergence rates and greater reliability. For details of these methods see [9, 10].
3. DUAL B O U N D A R Y ELEMENT METHOD The boundary element method is based on an integral equation transformation of the governing partial differential equations. In potential work, the governing equation is the Laplace equation V2u(x) = 0 for x in ~ (3) where f~ is the domain of the body and u(x) is the potential function. Boundary conditions are prescribed on the boundary of the body F and these can be written as u(x) = if(x)
for x on F 1
q(x) = q(x)
for x on F 2
(4)
where F 1 and F 2 are portions of the boundary F such that F = F 1 + F 2. The flux function, q(x), is the normal derivative of the potential function and tT(x) and ~(x) are the prescribed potential and flux values. From equation (3), use of a weighted residual technique gives an integral equation u(X t) +
fr
q*(X', x)u(x)
dr(x) = fr u*(X', x)q(x) dr(x).
(5)
In this equation X' is an internal source point, x is a boundary point and the terms u* and q* are known as the fundamental solutions, given for three-dimensional analysis by
u*(X', x) = q*(X', x)
1 4a'r 0u(x) On
1
Or
4 ~ 2 On
where r = Ix - X'I. This identity (5;) is valid for all points in the body, f~, including the boundary, F. Taking the internal point, X', to the boundary produces the boundary integral equation, c(x')u(x') + fr q*(x', x)u(x) dF(x) = fr u*(x', x)q(x) dF(x)
(6)
where c(x') is a jump term arising from taking the first integral to the boundary. For the
456
S.C.
MELLINGS
a n d M. H . A L I A B A D I
derivation of these equations, see Brebbia and Dominguez [11]. In order to evaluate the unknown potential and flux terms, it is necessary to discretize the boundary F into elements. In an element, the geometry as well as the potential and flux functions are assumed to have known variations. In this work, these variations are assumed to be quadratic and 8-noded quadrilateral elements are used. Transformation of the geometry in these elements to a local twodimensional system allows all the terms, in the integral equation, to be transformed into the local coordinate system. These terms can then be evaluated with the use of a standard Gaussian integration technique. c(x')u(x') + ~
3,=1 n = l
u 3"~; '
I
f_llq*(x ' , x(~, n))M'(~, n)J3"(x(~:, n)) dE dn =
q3" 3,=1 n = l
u*(x', x(~, r/))M"(~, r/)J3"(x(~¢, r/)) d~¢ dr/. 1
(7)
1
The terms M'(~, ~7) are the quadrilateral shape functions at the nth local node, J3"(x(~, 7/)) is the Jacobian of transformation and N is the number of boundary elements. In the Boundary Element Method, this equation is used at each node on the boundary, forming a set of linear equations in terms of the unknown potential or flux results at each node. This set of equations can be solved using a standard method, thus giving the unknown values.
3.1 Cracked bodies The standard application of the potential equation (7) to a body containing a crack does however cause difficulties. If equation (7) was used on both crack surfaces, two identical equations would be formed and the resulting matrix equation would be singular. The usual way of overcoming this problem is to sub-divide the body into two domains along the line of a crack. This allows each sub-region of the body to be computed independently and since the integration paths are now different, the system of equations is no longer singular. However, this introduces extra nodes along the division line, thus increasing the size of the system of linear equations. A new formulation, known as the dual Boundary Element Method, has been developed to overcome this problem without the use of sub-division. On one crack surface the standard, potential equation (6) is to be replaced by a new independent flux equation. This new equation is computed from the normal derivative of the potential equation. Using equation (5) for the potential value at an internal point X', the potential derivative at X', with respect to the coordinate direction, i, can be evaluated to be u,i(X') +
frq*(X',x)u(x ) dF(x)=
£ u*(X', x)q(x)dF(x)
(8)
where
1
u*(X', x) = 4zr----5 r,i
q*(X',x)=
1 (30r -n~).
4/t'r 3
c3-nr'i
It is necessary, as for the external boundary to evaluate the equations on a crack. An extra term is now produced due to the singularity on the opposite crack face. For the potential equation on a smooth crack boundary, the equation becomes
,1 u(x ,) + 1 u(x ") + f, q*(x', x)u(x) dr(x)
=
u*(x',x)q(x) dr(x)
(9)
3D flawidentification
457
where x" is the coincident point on the opposite crack surface. Similarly for the derivative potential equation on a smooth crack t!
1 u i(x')+ ~u i(x ) + fr q*(x', x ) u ( x ) d r ( x ) = fr u*(x', x)q(x)dr(x) 2
(lO)
where f represents the Cauchy principle value integral and ~ represents the Hadamard principle value integral. Since q(x)= u.i(x')ni(x'), multiplication of (10) by the normals at the source node will ~{ive the flux equation on a crack to be q(x') - ~ q(x") + ni(x')
q*(x', x)u(x) dr(x) = ni(x')
u*,i(xt, x)q(x) dr(x).
(11)
For a flux free crack q(x') = q(x") = 0. Discretization of equations (9) and (11) into boundary elements gives 1
1 ,, u(x , ) + 5"(x )+ ~
u "n
~_~J'+'
q*(x', x(~:, r/))M~(~, r/)Jr(x(sc, r/)) d~ dr/
3,=1 n = l
1
=
--1
u rn
u*(x', x(~, r/))M"(~, r/)J'(x(~, r/)) d~ dr/
"/=1 n = l
~q(x)-
1
q(x")+ni(x') ~ ~u~'n~ V=I n = l
q*(x',x(~,r/))Mn(~:,r/)JV(x(~,r/))d~dr/
a-1
= n;(x') ~= =mn=l qVn
(12)
1
1
f,;, 1
i u*(x', x(~:, r/))Mn(/.?, r/)JV(x(~, r/)) d~: dr/.
(13)
The Dual Boundary Element Method comprises of the use of both the potential equation, as given by equations (7) and (12), and the flux equation, as in equation (13). The potential equation (7) is to be used on the outer boundary of the problem where no problems arise due to coincident points. On one crack face the potential equation (12) is used, whilst the flux equation (13) is used on the other crack face. This gives, on application of the boundary conditions, a system of linear equations that may be solved for the unknown potential or flux values. This can be written in the matrix form as AX = F
(14)
where the vector F is formed from the application of the prescribed boundary potential and flux values, X is the vector of unknown potential or flux values and A is a coefficient matrix computed from the remaining terms. In the flaw identification method the Dual Boundary Element Method is used to compute the sensor values at tlhe boundary points, and if values are required at internal points, equations (5) and (8) are used to compute the values.
4. DESIGN SENSITIVITIES As mentioned previously, it is necessary in the minimization of the error function, to evaluate the design sensitivities. These are the sensitivities of the potential and flux values to changes to the design variables and are found by taking the derivatives, with respect these changes. There are three methods of computing these sensitivities; the finite difference approximation method, the material derivative method and the implicit, or direct, differentiation method.
458
S.C. MELLINGS and M. H. ALIABADI
In the first of these methods, values are computed using a finite difference approximation to a derivatives, for example a forward difference approximation is given by 0f(Z)
3Zm
-
f(Z
"~
~Zm) - f ( Z ) ~Zm
(15)
where 8z is a small forward difference step. This method gives an approximation that is simple to evaluate, but is computationally expensive and the accuracy is often very poor. The material derivative method is a variational method, a description of which can be found in Refs [14] and [15]. The third method mentioned above is the one used here and involves the direct differentiation of the boundary element equation(s) with respect to the design variable zm. The equations used in the Dual Boundary Element Method (7), (12) and (13) give implicit relationships for the potential and flux nodal values. Direct differentiation of these discretized equations would give new equations in terms of the potential, flux, potential derivative and flux derivative nodal values. Differentiation of the potential equation, on a standard boundary (7), gives c,m(x')u(x') + c(x')u,m(x') +
u,V~ ~'=1 n = l
+
+
J--I
J--I
u 3"n 3,=1 n = l
=
q (x, x(~:, r/))J,Vm(X(e, r/))m"(sc, r/) dgdr/
J--1 J--I
q,~
u*(x', x(sc, r/))M~(~, r/)J3"(x(sc, 7/)) d¢ dr/
y=l n=l
1
+
1
u*m(x', x(~, r/))J3"(x(~:, r/))Mn(~, 77) dsc dr/
q3" 3"=1 n = l
+
1
q,m(X, X(¢, r/))J3'(x(~, r/))Mn(~, 7-1)dg dr~
U3"n 3"=1 n = l
q*(x', x(~, r/))M"(sc, r/)J3"(x(~:, 7/)) d~ dr/ 1
1 ~ 1
q3"n 3"=1 n = l
u ( x , x(¢, r/))J3",m(X(g, r/))mn(¢, 7/) de dr/. J-1
(16)
J--1
On a smooth crack the derivative of the potential dual boundary element equation (12) gives
x)
1 , 1 ,,) + ~.~ - - Uyn ~U,m(X )+~U,m(X 3"=1 n = l
+
u rn 3"=1 n = l
+
f['f-' f['f[" 1
1
1
1
~ u 3"n 3"=1 n = l
=
q.~ v=l n=l
+ +
1
"-1
.!-I
q*(x',
X(¢, r/))M"(¢,
r/)J3"(x(e,
r/))dCd,
1
q*,,(x', x(~, r/))J3"(x(~, r/))Mn(sc, r/) d~ dr/
q*(x', x(sc, r/))JV,,(x(g, r/))Mn(g, r/) d~: dr/
u*m(x', x(¢, r/))J3"(x(¢, r/))m"(g, r/)de dr/
q3"n y=l n=l
1
u*(x', x(sc, r/))Mn(~:, r/)J3"(x(~:, 77)) d~: dr/ 1
q3"n 3"=1 n = l
f-l f_.
u (x, 1
-I
71)
dr/.
(17)
3D flaw identification
459
Finally, the flux Dual Boundary Element equation (13) differentiates to give
u z"
~q,.,(x )-~q~m(X")+ ni,.,(x') 3,=1 n = l
-[- ni(X t) ]~_~
I'~Y,~
~t=l n = l
+ ni(x') [~
1
f_~' q*(x', x(E, ~?))M"(E, ~7)JV(x(E, ~/)) dE d~? 1
u v"
"r=l n=l
q*(x', x(E, ~)M"(E, rt)JV(x(E, r/)) dE d*/ 1
1
, q,~m(X', X(E, n))JV(x(E, r/))M"(E, 77) dE dT/ 1 ~l
+ ni(x') /..~ "•
U Vn
T=I n=l
n=l
+ ni(x' ) .
1
+ ni(x') , , 3,=1 n = l
1
1
1
1
1
1
u*.,(x', x(E, n))JV(x(E, rl))Mn(E, "q) dE drl
qV-
-r=l n=l
+ n/(x")
1
f_~' u*(x', x(E, 7/))M"(E, ~?)JV(x(E, 7/))dE dr/
q,~
"),=1 n = l
f_, q*(x', x(E, r/))J,~m(x(E, ~))M"(E, ~l) dE d~? f - ' u*(x', x(E, ~))M"(E, r/)JV(x(~:, n))dE d~/
qV.
= ni, m(X' ) : ~ y=l
1
qV-
u*(x', x(E, r/))J,~.,(x(E, 7/))M"(E, r/) dE dr/.
(18)
In these equations the derivatives of the fundamental solutions u*, q*, u* and q*m are given by 1
u*,.(x', x) -
(19)
4xr2 r,m
1 (3Or ) q*,.(x', x) =4zrr ~ ~~nr, m -- [rlnl],m U*m(X',X) -q*.,(x',x)= 4zr ~ 1 (3Or ~n
(20)
1 4 ~ i(3rirm -- ro. )
(21) )
[5F'ir'm - ri'm] - 3[rtnt]'mri - 3nirm + ni'mr "
(22)
The potential and flux values can be computed using the DBEM as described earlier. Hence these equations implicitly relate the potential and flux derivatives. Derivative boundary conditions can be found from differentiation of the boundary conditions prescribed earlier, since it is assumed that they are independent of the crack shape and position. Hence the derivative boundary conditions are given by u,m(x) = 0
x E F~
q,m(x) = 0
x ~ F 2.
(23)
These derivative equations can now be used to provide a system of equations for each design variable. Use of the boundary conditions and the computed boundary potential and flux values allow the system to be written in matrix form as A X , , . = F,m - A . , , X .
(24)
In this system, X,,. is the vector of unknown flux and potential derivatives, (F,,,, - A,.,X) is the vector formed by the application of boundary conditions and A is the matrix of coefficients identical to that used in the DBEM analysis. In the flaw identification method, the system of
460
S.C.
MELLINGS
a n d M. H . A L I A B A D I
equations is solved with the LU Decomposition method, in which the A matrix is decomposed into an upper and a lower triangular matrix. The system is then solved using forward and backward substituion. Since the A matrix used here is identical to that in the previous matrix system (14) the decomposition can be preserved, rather than the matrix itself, and this can then be reused for each of the derivative analyses.
5. H Y P E R - S I N G U L A R I N T E G R A T I O N In the implementation of both the DBEM and the derivative formulation it is necessary to take special care when r ~ 0 (i.e. x ~ x ' ) , since the fundamental solutions q*, u*, q* and u'~,t and their derivatives all contain terms of the form O(1/r"), where n > 0. In the usual boundary element method, the rigid body condition is used to compute the weakly singular integrals that occur when the source and field points coincide. In the DBEM, however, this technique cannot be used on the crack, due to the presence of two equal and opposite singular terms, from each crack face, which cancel in the summation of row terms. Hence, some form of hyper-singular integration is required. The crack is assumed to be traction free, so all of the terms on the right-hand side of the equations (7)-(13) and (16)-(18) disappear. In the DBEM equations, the terms remaining on the left-hand side of the equation are given by the weakly singular kernel q* and the hyper-singular kernal q,*, i.e. -I
I 1=
f_ f_
q*(x', X(¢, r/))M"(~:, r/).D'(x(~:, 7/)) de dr/
(25)
, q*(x', x(~, r/))M"(s~, r/)Jr(x(~, r/)) d~ dr/.
(26)
1
~l
12 =
In the derivative equations, two extra integrals are required, these being the derivatives of I l and 12 given by '-I
J' =
f_ f_
[q*(x', x(~:, r/))JY(s¢, 7/) + q*(x', x(sc, r/))J~(x(¢, r/))]M(¢, r/) d~ dr/
(27)
[q*,.(x', x(~, r/))J~(~, r/) + q*(x', x(~, ~/))J,~,.(x(~, r/))]M(~, ~1) d~ dr/.
(28)
1
~l
j2 =
1
A transformation is made to a local polar coordinate system centred on the singular point x'. This allows the singularity to be removed by the exclusion of a small region of radius e, which is reduced by taking the limit as e ~ 0. I ~= lim e~O JO
12 = l i m c] TM ri~(o) e---*0 do
F~(p, O) dp dO
(29)
F2(p, 0) dp dO.
(30)
Jo~(e,O)
-'a(e,O)
In these expressions a(p, O) is a function representing the exclusion around the singularity and ~(0) the function for the external contour of Fv in polar coordinates, see Mikhilin [13]. The integrands are then given by
F'(p, e)=q*(p, e)M"(p, e)J~(p, e)p
(31)
0)= q~(p, O)M'(p, e)J='(p, e)p.
(32)
F2(p,
3D flawidentification
461
Similarly for the derivative integrals, the integrals can be transformed to polar coordinates giving
G~(p, 0) dp dO
(33)
G2(p, O) dp dO
(34)
G~(p, 0) = [q* M"J "/+ q*M"J~,.]p
(35)
G2(p, 0) = ['~*. L"I , t m M " J
(36)
jt = lim e-*O "0
" a(e,O)
j2 = lim e--*0 J 0
-'~(e.O)
with the integrands given by
~ + t'l ~ *,i M " J v, m J1^ l ~"
It is possible to rewrite the functions F 1, F 2, G I and G 2 in terms of a singular term, that can be integrated analytically, and a non-singular term, which can be computed using standard Gaussian integration, using ( F~(p, 0 ) = Fl(p, e)
FL~ F ~)_ ~ ( 0 ) 0 ) + - P
(37)
462
S.C. MELLINGS and M. H. ALIABADI
However, since the components of the normal can be written as ni = gdJ then Jn~ = g , hence this is evaluated in the expansion:
Ogk gk(p'O)=gk(rl)+ 0 (0¢;1 ~ cos O+-d-'~zslnO)+O(o2)" = gO + pg~(O) + 0(02).
(43)
Similarly the expansion of the shape functions can be written as
{OM" OM" O) O(02 ) M"(p,O)=Mn(n)+O o~ c°s 0 + - ~ - 2 sin + = M~ + pMT(O) + 0(02).
(44)
Manipulation of these terms, see Guggiani [12], allows the terms F 1 and F 2 to be expanded, giving Fl(p, 0) -
1 Or Mnj~ p 4zr On
- 4zl [F~-;(0)+O(1)] F2(p, 0 ) =
(45)
1 [3 Or J ~ r i - n i J ' ] M ' p 4~r 3 k On " 1 [ FZ_z(0)
- 4x[--~+
F~I(0)
]
0 + O(1)j.
(46)
The terms in the derivative integrals can be found in a similar way. The derivative of the distance components, ri,m is given by
rg,,.(p, O) = pAi,m(O ) -t- p2Bi,m(O ) -t- O(p 3)
(47)
where
Ai, m(O ) = Xi, lm COS 0 -b Xi,2m sin 0 Bi.,.(O) = Xi, ll,. cos 2 0 + 2X;,12,. cos 0 sin 0 + X~,22.. sin 2 0. This gives the derivative distance components to be
r m(p, O) = OAm(O) + 02Bm(O) + O(03)
(48)
where
Am(O) --
A,(0). A,.,,(0) A(O)
B,,( O)
ak,.,(O)Bk(O) + A~,(O)Bt,,.,(O) A(O)
Ak(O)Bk(O) A,.,(O)" "A'O ~, ))
The derivative of the Jacobian components is also required, with gk,,,, = J,,,,nt, + Jn~,.,,,:
{ Ogk,,. sin 0) + O(p z) gk,,,(O, O) = gk,m(rl) + 0~, 01~, COS 0 + Ogk.m 0~2 = gk,,O + Pgkm,(O) + O(02).
(49)
3D flawidentification
463
These terms can then be added in to give the integrands required, i.e. Gl(p, 0 ) = ~ ( 3 0 r r 4 7 0 .3
On ,m
-[rtnll,,,)M"JVp
- 4zl[G[p(0)+O(1)] G2(p, 0 ) =
1 (3 ~Or [5r, ir, m 4rtr------~
- ri,m] - 3 [ r l n t ] , m r , i -
(50)
3nir m + ni,,,r)M"JVp
1 [ G2_2(0) G2_I(0) ] =-4---~1_ p2 ~ - - p ~-0(1) .
(51)
6. RESULTS In the method presented here, it is necessary to examine both the accuracy of the design sensitivities as well as the effectiveness of the optimization procedure used. In the case of the sensitivity values., a comparison can be made with a series of finite difference approximations which will give an idea of the accuracy of the method. The performance of the flaw identification car, be tested using a known shape for the actual crack. The responses, which are usually required experimentally, can then be computed with the Dual Boundary Element Method with the known crack identity. The flaw found in the identification procedure can then be compared to the known flaw configuration.
6.1 Sensitivity values In this comparison, it is assumed that the cylinder contains a fiat, horizontal, penny shaped crack of unknown radius and location. Four design variables are used, these being the radius of the crack, R, and the three coordinate positions of the centre of the crack, i.e. x, y and z. Comparisons are to be made with the approximations given by the forward difference steps 8z,, = 0.1, 0.01, 0.001 and 0.0001. These approximations allow the percentage errors at each node, xn, of the ,cylinder mesh to be defined by
s(x') -
u f w d t x n ) -- U id ( x n ~ ,m \ ), ,m,~ ' * 1 0 0 % . ufWdtx~ ,it'/ k ]
(52)
In this expression, u f w d ( x n ) represents the forward difference approximation computed with the approximation formula given by equation (15) and uid(x ") is the value computed by the direct differentiation method, presented in this paper. A percentage error value is computed from the error at all N nodes as
E
VX e2(x") N
(53)
The value of this error function E is plotted on the graph for various sizes of the forward difference step, iSz,, and for each of the four design variables, x, y, z and R. As can be seen, in Fig. 1, the approximations were accurate and improved as the forward difference steps became smaller.
464
S.C. MELLINGS and M. H. ALIABADI i
4
i
i
i
i
i
ii
I
i
i
i
i
i
zxXC~ordin~t~o~~rackC~ntr~
ir
i
i
i
'~X~
V Y Coordinate of Crack Centre
/
z coo ,u..t ot c .ok C,r t ,
o
i
r
i
i
ii
/ ~
/
/
\
[] Radius of Cr
>
3
0
QJ
2
cO
0
1
r~
0
10-4
I
I
2
3
I
I
~
t
i irk~
t
t
10-a
2
S
Logarithmic
I
i
I
i
, ,
~
~ 2
10-s
Forward
Difference
3
10-t
Step Size
Fig. 1. Comparison of design sensitivities with forward difference approximations.
6.2 Inverse analysis 6.2.1 Identification of a penny shaped crack. Consider a cylinder of height 2 units and radius 1 unit, with a temperature of +100 units at the upper end and - 1 0 0 units at the lower end, with all the other sides being flux free. The cylinder actually contains a fiat, circular crack, of radius R = 0.5 with the crack centre at the origin (x = y = z = 0), as in Fig. 2. In the identification procedure, it is assumed that the crack is circular and fiat. The initial guesses made for the parameters of this crack are for a radius of R --- 0.1 and the central coordinates of x = 0.25, y =0.25, z =0.5. Potential measurements are taken at boundary sensor points which are positioned at the boundary nodes of the discretized cylinder where flux is zero. In the identification of this crack the BFGS Quasi-Newton Method was used to minimize the
u=+lO0
/
Y
/
q=O
/ /
X
Actual Initial
Crack j
u=-lO0 Fig. 2. Geometry of initial and actual cracks.
Crack
3D flaw identification
465
8th
Final Crack Initial Crack
lOth I t e r a t i o 5th Iteration
Fig. 3. Crack identities during iteration procedure.
error. Figure 3 shows the initial and final crack identities as well as those at the 4th and 7th iterations. The convergence both of the design variables (Fig. 4) and the sensor error normalized with respect to the initial value (Fig. 5) show that the crack can be identified successfully using only 20 iterations.
0.5,~ 0.4
l>
0.3
/
I V Y Coordinate of Crack Centre I
#
I 0 Z Coordinateof Crack Centre I
r=.=4
c~ a~ 0.2
t> la0
.,=.4
~D O.t
0.0
--0.|
I
o
I
z
I
I
4
I
I
e
I
l
a
I
I
1o
I
I
I
m
I
14
I
I
te
I
l
I
la
Iteration Number Fig. 4. Convergence of design variables in the identification of a penny shaped crack.
zo
466
S.C. MELLINGS and M. H. ALIABADI 100
•
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
0 ,--] tO-t 0 I-i I..i "~
e lO-e
of--I
a lO-a 0 Z 2
10-4 0
~
4
6
8
10
12
14
16
18
20
ILeration N u m b e r Fig. 5. Convergence of the normalized error in the identification of a penny shaped crack. F i g u r e 6 shows a c o m p a r i s o n b e t w e e n the f o u r o p t i m i z a t i o n m e t h o d s m e n i t i o n e d earlier. It can b e seen, that, a l t h o u g h the r a t e of c o n v e r g e n c e for the s i m p l e r m e t h o d s is l o w e r , all t h e p r o c e d u r e s u s e d i d e n t i f i e d the c r a c k a c c u r a t e l y . 6.2.2 Identification o f an elliptical crack. In the s e c o n d e x a m p l e , t h e a c t u a l c r a c k is the s a m e as in t h e p r e v i o u s e x a m p l e , b u t n o w t h e initial c r a c k is a s s u m e d to b e elliptical, with t h e c e n t r e of the elliptical c r a c k (xo, Yo, zo) is to b e i d e n t i f i e d a l o n g with t h e elliptical p a r a m e t e r s a a n d b defined by (X-Xo) 2 (y-yo) 2 + 1 and z = Zo. (54) b2 02
100
- q÷eer~es~ ih~scent Method
F
I
•
Conjugate
Gradient
Method
1 0-4
0
4
8
12
16
20
24
28
32
56
Iterotion Number Fig. 6. Comparison of the convergence of the normalized error for various optimization methods.
3D flaw identification
467
z u=+lO0
q=O
,21/~'--I
~
.....
Initial
I\Actual
Crack
Crack
u=-100 Fig. 7. Initial and final geometry for the identification of an elliptical crack.
The initial design is chosen with the design variables initially being a = 0.4,
b = 0.1,
Xo = -0.25,
Yo-- 0.25
and
Zo = -0.5
and flux values are measured at the sensor points. In this example the sensor points were located at the two ends of the cylinder where temperature was prescribed (see Fig. 7).
Initial Crack Final Crack
3rd Iteration 8th
\
Fig. 8. Crack identities during the identification of an elliptical crack.
468
S.C. MELLINGS and M. H. ALIABADI 0.5 0.4 0.3 0.2 ~"
0.1
~
0.0 -0.!
~
-0.2
l" El]~Parameter-b[
•~ -o.3 1 / I~
-0.4 L
]
? CentraICoordinate-x ~, C e n t r a l Coordinate - y I * Central Coordinate - z
-o.5 [ / -0.6
-0.7
0
4
8
12
16
20
24
28
32
[
J I
36
4
Iteration N u m b e r Fig. 9. Convergence of design variables in the identificationof an ellipticalcrack.
For this problem, the crack is again identified successfully, with the convergence of the design variables and the error normalized with respect to the initial value are shown in Figs 9 and 10. Figure 8 shows the final position of the crack as given by the interation procedure, along with the initial and actual crack identities.
7. CONCLUSION In this paper, a new method has been presented for the identification of cracks inside three-dimensional structures. The method is found to be highly efficient with the reuse of the
10
A
0
,
~
,
I
I
,
,
,
,
,
,
,
,
I
I
,
2
0 ,=1 t~ o
lO-t
It.1 8
lO-e
0 Z I
10-~
I
0
4
I
s
I
i
tz
I
i
t6
I
I
zo
I
I
z4
I
i
~
I
~
I
I
se
Iteration N u m b e r Fig. 10. Convergence of the normalized error in the identificationof an ellipticalcrack.
3D flaw identification
469
DBEM system matrix for each of the derivative analyses. The derivatives computed by this method are found to be highly accurate and the method is successful in the identification of an internal crack.
REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
G. N. VANDERPLAATS, A I A A J. 20, 992 (1981). A. FRIEDMAN and M. VOGELIUS, Indiana Univ. Math. J. 38, 527 (1989). K. BRYAN and M. VODELIUS, SlAM J. Math. Anal. 23, 950 (1992). M. TANAKA and Y. MASUDA, Engng Anal 3, 138. N. NISHUIMURA and S. KOBAYASHI, Int. J. Numer. Meth. Engng 7, 59 (1990). S. C. MELLINGS and M. H. ALIABADI, Engng Anal, Boundary Elements 12, 275 (1993). S. C. MELLINGS and M. H. AL1ABADI, Flaw identification using the boundary element method. To appear in hzt. J. Numer. Meth. Engng. K. BRYAN and M. VOGELIUS, Int. J. Engng Sci. 32, 579 (1994). R. FLETCHER, A review of methods for unconstrained optimization. In Optimization (Edited by R. FLETCHER), Academic Press, New York (1969). D. M. GREIG, Optimisation. Longman, London (1980). C. A. BREBBIA and J. DOMINQUEZ, Boundary Elements--An Introductory Course. Computational Mechanics Publications, McGraw-Hill, New York (1989). M. GUIGGIANI, G. KRISHNASAMY, T. J. RUDOLPHI and F. J. RIZZO, A general algorithm for the numerical solution of hypersingular boundary element integral equations. To appear in A S M E J. Appl. Mech. S. G. MIKHILIN, Multidimensional Singular Integral Equations. Pergamon Press, London (1965). A. LONGO, J. UNZUETA, E. SCAEIDT, A. A L V A R E Z and J. J. ANZA, Adv. Engng Software 16, 135 (1993). K. K. CHOI and E. J. HAUG, J. Structural Mech. U , 231 (1983). A. J. KASSAB, F. A. MOSLEHY and A. B. DARYAPURKAR, Nondestructive detection of cavities by and inverse elastostatics boundary element method. To appear in Engng Analysis Boundary Elements. (Received 24 October 1994; accepted 22 November 1995)
£S 3414-G