- Email: [email protected]
Inverse diffusion by boundary elements
EngineeringAnalysis with BoundaryElements15 (1995) 197-205 0955-7997(95)00018-6
ELSEVIER
Copyright © 1995 Elsevier Science Limited Printed in Great Britain. All rights reserved 0955-7997/95/$09.50
Inverse diffusion by boundary elements R. Pasquetti Laboratoire J.A. Dieudonn~, URA C N R S 168, Universit~ de Nice-Sophia Antipolis, Parc Valrose, 06108 Nice, France
& D. Petit Institut Universitaire des Systbme Thermiques Industriels, URA CNRS 1168, Centre de St Jdr6me, Avenue Escadrille Normandie Ni~men, 13397 Marseille, France
This paper considers the inverse problem which involves,for the steady or transient diffusion equation, the identification of boundary conditions from a given set of internal values. This problem is solved by using a boundary elementapproach, and an example of its application, concerned with a 2D test problem, is provided. Key words: Inverse problems, diffusion processes, boundary element method, regularization techniques.
Section 3 is devoted to the transient diffusion equation. As for the steady case, from the BIE, which now makes use of a time and space dependent fundamental solution, one deduces the basic discrete equations of the inverse problem. Moreover, for taking into account the damping and lagging effects of the diffusion operator, one utilizes future time steps. Then, one can produce a sequential algorithm, which only needs, at each time step, the straightforward resolution of an optimization problem similar to the one met in the steady case. However, as previously, the formulation of the functional involved in the optimization procedure depends on the path chosen for handling the inverse problem. In section 4, we conclude the paper with some further comments, in order to define how the boundary element approach presented in this paper can be handled in an efficient way.
1 INTRODUCTION
This paper is concerned with the ill-posed inverse problems which consist, for two- and three-dimensional diffusion, of the identification of boundary conditions from the function values at a set of internal points (see, for example, Ref. 1). The method is based on a boundary integral formulation, of which the main advantages are the usual ones of the BEM (boundary element method), especially no domain mesh is required and internal point positions can be freely defined. Obviously, the inverse method presented in this paper is not the only one practicable and quite different approaches are possible, generally based on inverse schemes (e.g. Refs 2, 3), optimization,4 optimal control s or filtering,6 or deconvolution techniques. 1'7 However, the present boundary element approach, first introduced in Ref. 8, is probably one of the most efficient. Section 2 is devoted to the steady problem. From the continuous formulation of the BIE (boundary integral equation), which makes use of the usual Green's function associated with the Laplace operator, one deduces first the basic discrete equations. Then, different ways for handling the inverse problem are discussed. All of them yield a regularized quadratic optimization problem, from which the explicit solutions can be easily produced. In order to check the efficiency of the method, we consider an academic 2D test problem for which an analytical solution is available.
2 BOUNDARY ELEMENT APPROACH FOR STEADY INVERSE DIFFUSION
Let us consider the steady diffusion equation: Au+g=0
(1)
in a 2D or 3D domain f~ of boundary F and with a given space dependent source term g. In order to produce the boundary integral formulation, this diffusion equation is 197
R. Pasquetti, D. Petit
198
rewritten as an integral equation with the help of the corresponding fundamental solution; for any point M of the closure of domain f~ the boundary integral equation (BIE) may read 9
dr =
qu* dr + J. gu*df2
(2)
with M: point of lP or f~, c = 1 if M is in f~ and c < 1 if M is on F (c = 0.5 i f F is smooth at M ) and q = Onu (On is the normal derivative operator). Function u*(M) is the usual space dependent Green's function associated with the Laplacian operator, and q*(M) is its normal derivative. The discrete form of eqn (2) is obtained by discretization of the boundary into N boundary elements (1 _< i _< N). Approaches of different approximation orders are then possible, but all of them finally lead to an equation of the form: Hu = Gq + s
(3)
where u and q are now the vectors the components of which are the collocation point values ui and qi, where H and G are matrices of dimension (N, N ) and where s is a vector associated with the source term. In the framework of direct problems, eqn (3) is associated with boundary conditions, of the Dirichlet, Neumann or mixed type when restricting ourselves to linear conditions. Then, rearranging system (3) leads to a standard system of N equations with N unknowns. Let us suppose now that (i) the boundary conditions are unknown on a part of the boundary and also that (ii) u is given at some internal points. Then, system (3) may only lead to a system of N equations with M unknowns, with M > N, but if one considers now the set of N' internal points where the u values are given, by expressing eqn (2) for each of them one an additive matrix equation similar to eqn (3):
A. Minimization of a quadratic norm of the residuals of the BIE
Let us introduce the residuals e and e' of eqns (3) and (4), respectively:
Hu = Gq + s + ~
(5)
u' + H'u = G'q + s' + ~'
(6)
and consider the functional J such as
N
N' e2 q- Z e;2 i=1 i=1
J = tel2 + 16-'12 = Z
(7)
Then, the inverse problem solution can be taken to be equal to the solution of the optimization problem: m i n J ( X ) , where all the unknown components of vectors u and q are gathered in vector X. After rearranging eqns (5) and (6) as
(8)
[ A A I X = [ B J + Ice]
where A, A' are (N, M ) , (N', M ) matrices, and B, B' are the corresponding vectors of dimension N and N', the solution of the optimization problem can be easily produced. Nevertheless, since the inverse diffusion problem is ill posed, such a solution is generally not smooth and so a regularization procedure is necessary to get satisfying results. 1° By introducing a regularization matrix R, associated with a norm or semi-norm of the solution, and a regularization coefficient #, used to adjust the amplitude of the regularization, one has to solve min [I,~'X - ~12 + #IRX[ 2 ]
(9)
with
[:l and
[."l
The solution of min J ( X ) reads
u' + H ' u = G'q + s'
(4) X = (~¢~"~¢ + # R ' R ) - ' . ~ " ~
where u' is the vector the components of which are the values u[ (1 < i < N ' ) of u at the internal points and H', G' are matrices of dimension ( N ' , N ) . Hereafter we assume that N + N ' > M, in order to be able to produce a solution to the inverse problem. Three different ways are now suggested for handling the inverse problem. They are based, respectively, on: • the minimization of a quadratic norm of the residuals of the BIE; • the minimization of a quadratic norm of the gaps between the given internal values and the cornputed ones; • the minimization of a quadratic norm of the gaps between all the given values (internal and boundary conditions) and the computed ones.
(10)
B. Minimization of a quadratic norm of the gaps between the given and computed internal values
The previous approach has the advantage of simplicity, but the minimization of the BIE residuals suffers from a lack of ground theory. On the other hand, from statistical considerations it is well justified to consider a quadratic norm of the gaps between the given and the computed internal values, when the measurement errors are supposed to be Gaussian and of zero mean value. 1,11 Let us introduce the gap vector, say s', between the computed values and the measurements: 41 = u' - ~/
(11)
Inverse diffusion by boundary elements where the circumttex is used for specifying the computed vector. The functional which has to be minimized simply reads N'
j = i ,i 2 = ~
eft2
(12)
i=l
Now, from eqn (4) one can say that the computed internal values read fi' = - H ' u + G'q + s'
(13)
Then, from eqns (11) and (13) one finishes again with eqn (6). Thus, the inverse problem is now the solution of min IA'X - B'I 2 AX = B
199
temperature than an internal one. Let us show that formulation B can be extended, assuming, for the sake of simplicity, that along F the boundary conditions either are unknown (on F0) or of Dirichlet type (on F1), in such a way that F = F0 u F1. The gap vector between the computed and prescribed values reads now
[;1] [u, 1 where the subscript 1 refers to I"1. The functional which has to be introduced reads
(21)
J = Iql 2 + I~'12
(14) F r o m eqns (3) and (4) one has
which is a quadratic optimization problem with linear constraints. As discussed in formulation A, in order to take out non-smooth solutions one must introduce a regularization term, and the problem which has to be solved then reads rain [[A'X - f i l E + AX=B
ulRX121
(15)
Such a problem has an explicit solution which can be produced by using the Lagrange multiplier technique. It consists in looking for the saddle-point of the quadratic functional:
J'(X, A)
fi = H - 1 G q + H - l s + c
(22)
where H - I is a generalized inverse of the singular matrix H and c is an additive constant vector. By taking into account eqn (20) and after rearranging, if necessary, eqn (3), one obtains
Iu°]=H u,
,23,
and u' = ( - H ' H - 1 G + G ' ) q - H ' H - l s + c' + sI + e I (24)
= IA'X - n'l 2 + ~ I R X I 2 + (A, AX - B)
(16) where (.,.) stands for the dot-product and ), for a Lagrange multiplier vector. At the optimum, the derivatives of J ' with respect to X and A vanish. This yields: X = M ( A ' t B ' - i1A t A)
(17)
AX - B = 0
(18)
where M = (A'tA ' + #RtR) -1. By substituting X of eqn (17) into eqn (18) one can express the Lagrange multiplier vector and finally, coming back to eqn (17), one obtains for the optimal value of X: X = M[A'tB ' - A t ( A M A t ) - I ( A M A ' t B ' - B)]
,
fi' = - H ' f i + G'q + s
(19)
Although more complicated than formulation A, formulation B can be used in a very efficient way since the solution is still explicit.
C. Minimization of a quadratic norm of the gaps between all the imposed values and the computed ones Although satisfactory from the theoretical point of view, formulation B supposes that the boundary conditions are perfectly known, but such an assumption is often not realistic. In the framework of heat conduction, it is, for example, more difficult to measure a superficial
where vector c' corresponds to the additive constant. Introducing the vector X of the unknowns (i.e. the components of vectors q and u0 and the additive constant value appearing in vectors c and c') results in the following optimization problem:
min[A:l
,2,,
A 0 X = 80 where the subscripts 0 and 1 refer again to F0 and r I and the superscript prime to the internal points and where the different matrices A0, AI, A t and vectors B0, B1, B' can be easily identified. As eqn (25) is quite similar to eqn (14), its solution is similar to the one described in formulation B. Example of application Let us come now to an example of application, based on the use of formulation A which is the one implemented in our software. Far from the industrial context (see, e.g. Ref. 12 for a 3D example of application issued from industry) we are going to consider the test problem introduced in Ref. 13, but slightly modified as justified later. Considering the Laplace equation in the square domain of Fig. 1, the aim is to recover the analytical solution u = x 2 - yZ, from Dirichlet conditions along
R. Pasquetti, D. Petit
200
(1,1)
(0.25,0.5) (0.75,0.5)
X r
Fig. 1. Schematic of the domain.
the three unitary sides x = 0, x = 1 and y = 0, and from given values o f u between the two points (x = 0.25, y = 0-5) and (x = 0.75, y = 0-5). Such a problem is an inverse one since the b o u n d a r y conditions are u n k n o w n along the side y = 1; clearly, along this side u = x 2 - 1 and q = Oyu = - 2 , f r o m the analytical solution. First, one has to mention that g o o d results cannot be obtained without regularization. This points out that the considered test problem is severely ill posed, since the F o r t r a n double precision appears here to be not sufficient. In Fig. 2 (a, b, c) are given the c o m p u t e d values o f u(y = 1) and q(y = 1) for three different values o f the coefficient # and when the L2(0 < x < 1,y = 1) n o r m o f the function q is used for the regularization, i.e. in discretized form: # ~
q~,
i, such as y = 1
(26)
i
0.5
0.5
t
i
I
I
I
I
i
I
I.E-16
0
0
i
i
i
i
i
i
I.E-15
i
i
u
-0.5
-0.5 -1
-1-2" 5 ~ ~ q ~
-1.5 -2
-2.5
-2.5
- 3
-3 0.2 0.5
I
0.4
t
,
I
0.6 I
,
I
I
I
I
0.5
1.E-15
•
i
,
,
0.4
I
|
1.E-14
0
0 -0.5
,
0.2
0.8
i
,
0.6
I
I
w
0.8
l
I
I
I
U
-0.5
/
-1
-1
-1.5
-1.5
-2
-2
q
-2.5
q
-2.5 -3
-3
o 0.5
o'.2
o'.4
l.E/14.
'
0
,
,
=
,
,
,
a U
0.5
-1
-1
-1.5
-1.5
-2
-2
q '
0',2
'
0'.4
0'.6
'
018
'
x (y=l) Fig. 2. Computed
I
I
I
I
o18 I
I
I
q
-2.5
-3
I
0'.8
0 -0.5
0
I
0'4
I.E-13
-0.5
-2.5
o'.2
o'.e
values of u(x,l) and q(x, I) with the
Euclidean norm for regularization: (a) # = 1076; (b) # = 10 is; (C) # ~- 10 -14.
-3
0
0'.2
' 0'.4
' o'.5
' o18
'
x (y=Z)
Fig. 3. Computed values of u(x,l) and qlx, 1) with the weighted norm for regularization: (a~3# = 10~5; (b) # = 10-14;
(c) ~ = lO
.
Inverse diffusion by boundary elements These results have been obtained with N = 160 (= 4 x 40) boundary elements, which is the most severe situation considered in Ref. 13, and N ' = 41 internal points. They demonstrate the efficiency of the regularization technique, which one permits to obtain an accurate solution for # = 10-15. Nevertheless, one observes that when the coefficient # is increased the variations of q become smoother, but quite false in the vicinity of the comers, where the regularization term tends to vanish the q value. In order to overcome this difficulty, without using the facility of a regularization of higher order (see the remark below), the following weighted norm has been tried: 2 2
i, such as y
# 2-, Pi qi, i
(27)
1
with the coefficients Pi calculated, at the middle of each 0.5
boundary element, from the parabola: p = 4(e + x)(1 + e - x)
0.5
|
1.E- 11
0
0
-0.5
-0.5
-1
-1
-1.5
-1,5
-2
-2 -2.5
q
-3
o'.,
0.5
(28)
where e (_> 0) has been taken as equal to half of the boundary element length. The results obtained when such a regularization term is used are presented in Fig. 3 (a, b, c), for three values of the regularization coefficient. One can observe the improvement of the solution, which is now smoother, except in the vicinity of the two comers where the regularization term acts in a weaker way, since eqn (28) yields a ratio of the maximum and minimum values of coefficients Pi equal to 10.5. In the considered test problem the internal values are supposed to be perfectly known, i.e. only limited by the Fortran double precision. However, in a realistic context these internal values are generally affected by measurement errors. In order to simulate measurements, the internal values have been altered by Gaussian errors.
1.E-13
-2.5
201
z
!
a
|
ole |
z
-3
o;e i
I
|
o12
0.5
z
t
o'.e
0'.4 |
i
i
z
olo |
|
i
I.E-I0
I.E-12 0
0
U
-0.5
-0.5
-1
-1
-1.5
-1.5 -2
-2
q
-2.5 -3
0'4
o'.a
0.5
J
t
t
i
o'e i
i
-3
ole i
J
i
o12
0.5
|
i
o14 t
|
ols i
z
ole i
t
z
! ,E-.09
I.E-II
0
q
-2.5
0
-0.5
-0.5
-1
-1
-1.5
-1.5 -2
-2
q -2.5
-2.5 -3
o
0'.2'
o'.,
' o'.6
olo
x (y=l) Fig. 4. Computed values of u(x, 1) and q(x, 1) with measurement errors such as 2a = 10-°: (a) # = 10-13; (b) -- 10 -12, (c) # = 10TM (weighted norm).
-3
q
o
' o'.2 ' o'.4
o'.s ' o18 '
x (y=l) Fig. 5. Computed values of u(x, 1) and q(x, 1) with measurement errors such as 2or= 10-~: (a) # = 10-ll; (b) # = 10-1°; (c) # = 10 -9 (weighted norm).
202
R. Pasquetti, D. Petit 0.5
for the regularization can yield artificially a very accurate result, by taking the value of coefficient # sufficiently high. In order to obtain a more realistic result, it is preferable to assume the boundary conditions unknown at y = 1, where q ¢ 0. Moreover, we have restricted ourselves to a regularization of zeroth order, since by using a regularization of higher order the same phenomenon occurs. Indeed, the L 2 (0 < x < 1, y = 1) of all the derivatives of q with respect to x vanish, owing to the constant value q - - - 2 of the flux density along the side 0 < x < 1, y = 1.
I
I.E-10
0 -0.5
/
-1 -1.5 -2 -2.5 -3 0 0.5
0.2 I
I
0.4 I
I
0.6 I
l
0.8 I
I
I
I.E-09
3 B O U N D A R Y E L E M E N T A P P R O A C H FOR TRANSIENT INVERSE DIFFUSION
0 -0.5 -1
Although transient inverse problems are much more difficult to solve than steady ones, the approach that we are going to describe for the transient case is very similar to the one used for the steady problems. Such an approach has already been described in Ref. 14. Briefly, let us recall that, from the transient diffusion equation:
-1.5 -2 -2.5
q
-3 0.5
o'2 I
I
0'4 I
I
46 I
i
0;8 I
i
I
1 - Ot u = Au + g
(29)
(/
1 .E-08
0
U
where t stands for time and a for the diffusivity, one can deduce the following BIE:
-0.5 -1
CUM,IF +
-1.5 -2 -2.5
auq* dF dt = to
F
aqu* dI" dt to
F
+ J ; i J ~ agu*df, d t + l n U o u * d f f t
(30)
q
-3 '
012
'
o'.4 ' o'.6
'
018
'
x (y=1)
Fig. 6. Computed values of (.;~1t and q(x 1) with measurement errors such as 2~r u : (a) # = 10 o: (b) # ----- 1 0 - 9 ; (C) p, = 10 -8 (weighted norm). The Gaussian distribution is supposed to be centred and the standard deviation cr is taken as equal to half of the measurement precision. The results obtained with 2or = 10 -6, l0 5 and 10 4 are respectively presented in Figs 4 (a, b, c), 5 (a, b, c) and 6 (a, b, c) for different values of the regularization coefficient and when using the weighted norm defined by eqns (27) and (28). They clearly demonstrate the high sensitivity of the inverse problem solution to the measurement errors, but one observes that a right choice of the regularization coefficient can lead to satisfying results. Remark In the problem introduced in Ref. 13 the boundary conditions are unknown along the unitary side 0 < x < 1, y = 0 (instead of y = 1 in the present study). However, along this side, the exact value of the flux density being equal to q = 0, the use of eqn (26) (or eqn (27))
where u0 is now the initial condition and where to and tF stand for the initial and final time. In such a BIE, the fundamental solution u* is a time and space dependent Green function: u*
(47rat)S/2 exp
H(T)
(31)
where H is the Heaviside function, ~- = tF -- t, S is the space dimension, and r is the distance to point M. The discrete form of eqn (30) is obtained by discretization of the time ( f : time index, 0 _
(32)
F = Gq F ÷ s F
where u g a n d q F are the vectors the components of which are u~ and qi at time tF; H and G are new matrices of dimension (N, N ) , as in eqn (3), and s F is a vector which involves the source term and the initial conditions as well as the u f a n d q f vectors for f < F. Considering now N ' internal points, one obtains an equation similar to eqn (4): u IF ÷ H ' u F
= G,qF
+ s tF
(33)
Inverse diffusion by boundary elements
203
u 'F is the vector the components of which are, at time
4.1 Regularization
tF, the u' ( l < i < N ' ) values of u at the internal points, and I-I', G' are matrices of dimension (N', N). As for the steady case, at this point one has to suppose that the number of equations is greater than the number of unknowns, equal to M, which is problem dependent:
The characteristic parameters of the regularization being numerous, some experience is needed from getting satisfying results. For us:
N+NI>M. Nevertheless, an approach similar to the one described for the steady case would be inefficient, owing to the damping and lagging effects of the transient diffusion operator. Efficient inverse methods make use of 'future time steps': in order to determine what happens at time tF one uses data at times tF+r, 0
H
= G u
u
R
1
Eq]
+ S
q
(34)
R
Lu;+R J
qF+R + S'
(35)
where H and G are block triangular matrices of dimension (N(R + 1), N(R + 1)), H ' and G ' block triangular matrices of dimension (N'(R + 1), N(R + 1)), S is a vector o f dimension N(R+ 1), and S' is a vector of dimension N ' ( R + 1). All of them can be easily identified. Equations (34) and (35) are now the equivalents of eqns (3) and (4) for the steady case, and thus can be handled in a very similar manner in order to solve the inverse problem at time t e. It must be mentioned that the computed values of U F ÷ r a n d qF+r, r ~ O, are simply rejected from the solution of eqns (34) and (35), their definitive values being calculated at time t F+'. For the regularization, from our experience it appears that a space regularization over all the future time steps leads to satisfactory results. This yields through substitution by the regularization matrix R of the steady problem, a block diagonal matrix R constituted with ( R + 1) blocks and such as: R = diag ( R , . . . , R). Two examples of application of the boundary element approach presented in this section can be found in Refs 8 and 14. Moreover, in Ref. 8 a comparison with the superposition method as described in Ref. 1 is given.
4 COMPLEMENTARY REMARKS In this section we will discuss the use of the present boundary element approach and more generally inverse diffusion problems.
• In discrete form, the regularization term usually expresses the Euclidean norm of the flux density q along the part of the boundary where the boundary conditions are unknown. Our numerical tests have clearly shown that the use of q as control variable is more efficient than the use of u. Nevertheless, in order to improve the results, it is sometimes interesting to choose a regularization of higher order. In this case, we use the Euclidean norm of a finite difference approximation of a tangential derivative of q. • The value of the regularization coefficient can be easily adjusted when the exact result is known, but this is not the case for real applications. Without information on the final result, the suitable procedure involves adjusting the value of # in such a way that the gaps between the computed values and the measurements are of the magnitude of the measurement precision. • In the framework of transient problems one has also to adjust the number of future time steps, which depends on the distances from the unknown boundary conditions to the measurement points. From our experience, in realistic applications the parameter R can generally be taken low, e.g. R < 3. • When a sharp corner occurs in the part of the boundary where the boundary conditions are unknown, it is generally difficult to obtain satisfying results in the vicinity of this corner. It is then worthwhile to use the local regularization technique introduced in Ref. 14. It consists of adding to the functional of the optimization problem, an additive local regularization term which expresses the finite difference approximation of one of the tangential derivatives of q.
4.2 Nonlinear diffusion The diffusion equations considered in this paper have been assumed as linear. Concerning nonlinearities in the domain, let us recall here that a nonlinear steady diffusion equation such as, with A for the conductivity: V.(A(u)Vu) + g = 0 may be handled by using the Kirchhoff transformation (see, e.g. Ref. 15). The same remark applies to the transient diffusion equation, but it should be mentioned that in the transient case, the use o f the Kirchhoff transformation can lead to approximations in the results, unless the diffusivity does not depend on u.
204
R. Pasquetti, D. Petit
Concerning nonlinearities at the boundary, one can proceed with the numerical techniques used for direct problems. For us, we utilize an iterative technique, based on local linearizations, i.e. specifics of each boundary element of the nonlinear boundary conditions. ~5 4.3 Characteristic parameters of the diffusion process One of the main difficulties of the inverse problem lies in the fact that the characteristic parameters of the diffusion process, especially the conductivity and the diffusivity, are generally known only approximately. This can lead to serious difficulties, owing to an approximation in these parameters which does not induce centred aleatory errors, as it is usually assumed. The influence of such approximations can be simply checked, by solving the inverse problem with altered values of the characteristic parameters. Thus, in Refs 12 and 14, the inverse problems are solved with a conductivity higher or lower than the conductivity used for the direct problem, and, roughly speaking, this induces a translation in the values of the flux densities, which become higher or lower than the nominal ones. Consequently, when the gaps between the computed values and the measurements are clearly not of the centred Gaussian aleatory type, then one may expect some lack of accuracy in the values of the characteristic parameters.
5 CONCLUSION We have been interested in the inverse problem which involves, for both the steady and transient diffusion equations, the identification of boundary conditions from internal measurements. The boundary element approach which has been suggested is applicable to 2D and 3D geometries. It can yield to space dependent (steady case) or time and space dependent (transient case) boundary conditions, without any assumption on the shape of the solution (constant, linearly variable, etc.). The approach uses regularization techniques, in order to overcome the ill-posed feature of the inverse problem, and, in the transient case, future time steps, owing to the lagging effects of the diffusion operator. Moreover, nonlinear boundary conditions can be supported for the classical superposition method, as well as variable conductivity. For transient problems, the algorithm is a sequential one, which is necessary for large industrial applications. In this paper we have especially focused on: • The quadratic functional, which has to be minimized for solving (at each time step in the transient case) the inverse problem. Three different approaches have been suggested. The first one
(formulation A) is strongly associated with the boundary element formulation, since it expresses the Euclidean norm of the residuals of the BIE. On the other hand, the two others (formulations B and C) result only from statistical considerations, by using the gaps between the computed and the measured values. Formulation B especially uses the functional introduced in Ref. 1 and it may be noticed that the improved functional of formulation C is not accessible to the classical superposition method. The regularization techniques, which are generally necessary to get an accurate solution. Considering the test problem introduced in Ref. 13, for comparing with the Laplace equation different inverse methods, a satisfactory result has been produced by using as regularization term the Euclidean norm of the normal derivative. Moreover, it has been shown that the use of a weighted norm can still improve this result and that the method remains efficient when measurement errors are introduced. Finally, some complementary remarks, taken from our experience of the solution of inverse diffusion problems by boundary elements, have been mentioned.
REFERENCES 1. Beck, J. V., Blackwell, B. & St. Clair, C. R. Inverse Heat Conduction, Ill Posed Problems. Wiley Interscience. New York, 1985. 2. Weber, C. F. Analysis and solution of the ill-posed inverse heat conduction problem. Int. J. Heat Mass Transfer, 1981, 24, 1783-92. 3. Raynaud, M. & Bransier, J. A new finite difference method for the nonlinear inverse heat conduction problem. Num. Heat Transfer, 1986, 9(1), 27-42. 4. Alifanov, O. M. et al. Boundary inverse heat conduction problem in parametric form. In 2nd Int. Conf. on Advanced Computer Methods in Heat Transfer, Milan, 1992, Vol. 1, pp. 479-89. 5. Guerrier, B. & Brnard, C. Identification of time and space dependent boundary conditions in a 2D thermal system. In 2nd Int. Conf. on Advanced Computer Methods in Heat Transfer, Milan, 1992, pp. 448-65.
6. Fort, C. Estimation de paramrtres et m&hode inverse en thermique. Application ~i la drtermination de la variation du flux parirtal dans une chambre de combustion. Thrse de doctorat, University of Poitiers, 1989. 7. Kurpisz, K. & Nowak, A. J. BEM approach to inverse heat conduction problems, Engng Anal. with Bound. Elem., 1992, 10, 291-7. 8. Pasquetti, R. & Le Niliot, C. Boundary element approach for inverse heat conduction problems: Application to a bidimensional transient numerical experiment. Num. Heat Transfer, Part B, 1991, 20, 169-89. 9. Wrobel, L. C. & Brebbia, C. A. The boundary element method for steady-state and transient heat conduction. In Numerical Methods in Thermal Problems, Vol. 1, ed. R. W. Lewis & K. Morgan. Pineridge Press, Swansea, 1979.
Inverse diffusion by boundary elements 10. Thikonov, A. & Arsenine, V. M~thode de Rdsolution des Problbmes Mal Posds. Editions de Moscou, 1976. 11. Eykhhoff, P. System Identification, Parameter and State Estimation. Wiley-Interscience, New York, 1974. 12. Petit, D., Debray, V., Le Niliot, C. & Pasquetti, R. Identification of local heat transfer coefficient using a boundary element formulation: In 2nd Int. Conf. on Advanced Computer Methods in Heat Transfer, Milan, 1992, Vol. 1, pp. 467-77. 13. Ingham, D. B. Improperly posed problems in heat transfer. In Boundary Element Methods in Heat Transfer, ed. L. C. Wrobel & C. A. Brebbia. Computational Mechanics
205
Publications, Elsevier Applied Science, London, 1992, pp. 269-94. 14. Pasquetti, R. & Petit, D. Inverse heat conduction problems with boundary elements: analysis of a comer effect. Engng Anal. with Bound. Elem., 1994, 13, 33-62 15. Pasquetti, R., Caruso, A. & Wrobel, L. C. Transient problems using time dependent fundamental solutions. Boundary Element Methods in Heat Transfer, ed. L. C. Wrobel & C. A. Brebbia. Computational Mechanics Publications, Elsevier Applied Science, London, 1992, pp. 33-62.
ELSEVIER
Copyright © 1995 Elsevier Science Limited Printed in Great Britain. All rights reserved 0955-7997/95/$09.50
Inverse diffusion by boundary elements R. Pasquetti Laboratoire J.A. Dieudonn~, URA C N R S 168, Universit~ de Nice-Sophia Antipolis, Parc Valrose, 06108 Nice, France
& D. Petit Institut Universitaire des Systbme Thermiques Industriels, URA CNRS 1168, Centre de St Jdr6me, Avenue Escadrille Normandie Ni~men, 13397 Marseille, France
This paper considers the inverse problem which involves,for the steady or transient diffusion equation, the identification of boundary conditions from a given set of internal values. This problem is solved by using a boundary elementapproach, and an example of its application, concerned with a 2D test problem, is provided. Key words: Inverse problems, diffusion processes, boundary element method, regularization techniques.
Section 3 is devoted to the transient diffusion equation. As for the steady case, from the BIE, which now makes use of a time and space dependent fundamental solution, one deduces the basic discrete equations of the inverse problem. Moreover, for taking into account the damping and lagging effects of the diffusion operator, one utilizes future time steps. Then, one can produce a sequential algorithm, which only needs, at each time step, the straightforward resolution of an optimization problem similar to the one met in the steady case. However, as previously, the formulation of the functional involved in the optimization procedure depends on the path chosen for handling the inverse problem. In section 4, we conclude the paper with some further comments, in order to define how the boundary element approach presented in this paper can be handled in an efficient way.
1 INTRODUCTION
This paper is concerned with the ill-posed inverse problems which consist, for two- and three-dimensional diffusion, of the identification of boundary conditions from the function values at a set of internal points (see, for example, Ref. 1). The method is based on a boundary integral formulation, of which the main advantages are the usual ones of the BEM (boundary element method), especially no domain mesh is required and internal point positions can be freely defined. Obviously, the inverse method presented in this paper is not the only one practicable and quite different approaches are possible, generally based on inverse schemes (e.g. Refs 2, 3), optimization,4 optimal control s or filtering,6 or deconvolution techniques. 1'7 However, the present boundary element approach, first introduced in Ref. 8, is probably one of the most efficient. Section 2 is devoted to the steady problem. From the continuous formulation of the BIE (boundary integral equation), which makes use of the usual Green's function associated with the Laplace operator, one deduces first the basic discrete equations. Then, different ways for handling the inverse problem are discussed. All of them yield a regularized quadratic optimization problem, from which the explicit solutions can be easily produced. In order to check the efficiency of the method, we consider an academic 2D test problem for which an analytical solution is available.
2 BOUNDARY ELEMENT APPROACH FOR STEADY INVERSE DIFFUSION
Let us consider the steady diffusion equation: Au+g=0
(1)
in a 2D or 3D domain f~ of boundary F and with a given space dependent source term g. In order to produce the boundary integral formulation, this diffusion equation is 197
R. Pasquetti, D. Petit
198
rewritten as an integral equation with the help of the corresponding fundamental solution; for any point M of the closure of domain f~ the boundary integral equation (BIE) may read 9
dr =
qu* dr + J. gu*df2
(2)
with M: point of lP or f~, c = 1 if M is in f~ and c < 1 if M is on F (c = 0.5 i f F is smooth at M ) and q = Onu (On is the normal derivative operator). Function u*(M) is the usual space dependent Green's function associated with the Laplacian operator, and q*(M) is its normal derivative. The discrete form of eqn (2) is obtained by discretization of the boundary into N boundary elements (1 _< i _< N). Approaches of different approximation orders are then possible, but all of them finally lead to an equation of the form: Hu = Gq + s
(3)
where u and q are now the vectors the components of which are the collocation point values ui and qi, where H and G are matrices of dimension (N, N ) and where s is a vector associated with the source term. In the framework of direct problems, eqn (3) is associated with boundary conditions, of the Dirichlet, Neumann or mixed type when restricting ourselves to linear conditions. Then, rearranging system (3) leads to a standard system of N equations with N unknowns. Let us suppose now that (i) the boundary conditions are unknown on a part of the boundary and also that (ii) u is given at some internal points. Then, system (3) may only lead to a system of N equations with M unknowns, with M > N, but if one considers now the set of N' internal points where the u values are given, by expressing eqn (2) for each of them one an additive matrix equation similar to eqn (3):
A. Minimization of a quadratic norm of the residuals of the BIE
Let us introduce the residuals e and e' of eqns (3) and (4), respectively:
Hu = Gq + s + ~
(5)
u' + H'u = G'q + s' + ~'
(6)
and consider the functional J such as
N
N' e2 q- Z e;2 i=1 i=1
J = tel2 + 16-'12 = Z
(7)
Then, the inverse problem solution can be taken to be equal to the solution of the optimization problem: m i n J ( X ) , where all the unknown components of vectors u and q are gathered in vector X. After rearranging eqns (5) and (6) as
(8)
[ A A I X = [ B J + Ice]
where A, A' are (N, M ) , (N', M ) matrices, and B, B' are the corresponding vectors of dimension N and N', the solution of the optimization problem can be easily produced. Nevertheless, since the inverse diffusion problem is ill posed, such a solution is generally not smooth and so a regularization procedure is necessary to get satisfying results. 1° By introducing a regularization matrix R, associated with a norm or semi-norm of the solution, and a regularization coefficient #, used to adjust the amplitude of the regularization, one has to solve min [I,~'X - ~12 + #IRX[ 2 ]
(9)
with
[:l and
[."l
The solution of min J ( X ) reads
u' + H ' u = G'q + s'
(4) X = (~¢~"~¢ + # R ' R ) - ' . ~ " ~
where u' is the vector the components of which are the values u[ (1 < i < N ' ) of u at the internal points and H', G' are matrices of dimension ( N ' , N ) . Hereafter we assume that N + N ' > M, in order to be able to produce a solution to the inverse problem. Three different ways are now suggested for handling the inverse problem. They are based, respectively, on: • the minimization of a quadratic norm of the residuals of the BIE; • the minimization of a quadratic norm of the gaps between the given internal values and the cornputed ones; • the minimization of a quadratic norm of the gaps between all the given values (internal and boundary conditions) and the computed ones.
(10)
B. Minimization of a quadratic norm of the gaps between the given and computed internal values
The previous approach has the advantage of simplicity, but the minimization of the BIE residuals suffers from a lack of ground theory. On the other hand, from statistical considerations it is well justified to consider a quadratic norm of the gaps between the given and the computed internal values, when the measurement errors are supposed to be Gaussian and of zero mean value. 1,11 Let us introduce the gap vector, say s', between the computed values and the measurements: 41 = u' - ~/
(11)
Inverse diffusion by boundary elements where the circumttex is used for specifying the computed vector. The functional which has to be minimized simply reads N'
j = i ,i 2 = ~
eft2
(12)
i=l
Now, from eqn (4) one can say that the computed internal values read fi' = - H ' u + G'q + s'
(13)
Then, from eqns (11) and (13) one finishes again with eqn (6). Thus, the inverse problem is now the solution of min IA'X - B'I 2 AX = B
199
temperature than an internal one. Let us show that formulation B can be extended, assuming, for the sake of simplicity, that along F the boundary conditions either are unknown (on F0) or of Dirichlet type (on F1), in such a way that F = F0 u F1. The gap vector between the computed and prescribed values reads now
[;1] [u, 1 where the subscript 1 refers to I"1. The functional which has to be introduced reads
(21)
J = Iql 2 + I~'12
(14) F r o m eqns (3) and (4) one has
which is a quadratic optimization problem with linear constraints. As discussed in formulation A, in order to take out non-smooth solutions one must introduce a regularization term, and the problem which has to be solved then reads rain [[A'X - f i l E + AX=B
ulRX121
(15)
Such a problem has an explicit solution which can be produced by using the Lagrange multiplier technique. It consists in looking for the saddle-point of the quadratic functional:
J'(X, A)
fi = H - 1 G q + H - l s + c
(22)
where H - I is a generalized inverse of the singular matrix H and c is an additive constant vector. By taking into account eqn (20) and after rearranging, if necessary, eqn (3), one obtains
Iu°]=H u,
,23,
and u' = ( - H ' H - 1 G + G ' ) q - H ' H - l s + c' + sI + e I (24)
= IA'X - n'l 2 + ~ I R X I 2 + (A, AX - B)
(16) where (.,.) stands for the dot-product and ), for a Lagrange multiplier vector. At the optimum, the derivatives of J ' with respect to X and A vanish. This yields: X = M ( A ' t B ' - i1A t A)
(17)
AX - B = 0
(18)
where M = (A'tA ' + #RtR) -1. By substituting X of eqn (17) into eqn (18) one can express the Lagrange multiplier vector and finally, coming back to eqn (17), one obtains for the optimal value of X: X = M[A'tB ' - A t ( A M A t ) - I ( A M A ' t B ' - B)]
,
fi' = - H ' f i + G'q + s
(19)
Although more complicated than formulation A, formulation B can be used in a very efficient way since the solution is still explicit.
C. Minimization of a quadratic norm of the gaps between all the imposed values and the computed ones Although satisfactory from the theoretical point of view, formulation B supposes that the boundary conditions are perfectly known, but such an assumption is often not realistic. In the framework of heat conduction, it is, for example, more difficult to measure a superficial
where vector c' corresponds to the additive constant. Introducing the vector X of the unknowns (i.e. the components of vectors q and u0 and the additive constant value appearing in vectors c and c') results in the following optimization problem:
min[A:l
,2,,
A 0 X = 80 where the subscripts 0 and 1 refer again to F0 and r I and the superscript prime to the internal points and where the different matrices A0, AI, A t and vectors B0, B1, B' can be easily identified. As eqn (25) is quite similar to eqn (14), its solution is similar to the one described in formulation B. Example of application Let us come now to an example of application, based on the use of formulation A which is the one implemented in our software. Far from the industrial context (see, e.g. Ref. 12 for a 3D example of application issued from industry) we are going to consider the test problem introduced in Ref. 13, but slightly modified as justified later. Considering the Laplace equation in the square domain of Fig. 1, the aim is to recover the analytical solution u = x 2 - yZ, from Dirichlet conditions along
R. Pasquetti, D. Petit
200
(1,1)
(0.25,0.5) (0.75,0.5)
X r
Fig. 1. Schematic of the domain.
the three unitary sides x = 0, x = 1 and y = 0, and from given values o f u between the two points (x = 0.25, y = 0-5) and (x = 0.75, y = 0-5). Such a problem is an inverse one since the b o u n d a r y conditions are u n k n o w n along the side y = 1; clearly, along this side u = x 2 - 1 and q = Oyu = - 2 , f r o m the analytical solution. First, one has to mention that g o o d results cannot be obtained without regularization. This points out that the considered test problem is severely ill posed, since the F o r t r a n double precision appears here to be not sufficient. In Fig. 2 (a, b, c) are given the c o m p u t e d values o f u(y = 1) and q(y = 1) for three different values o f the coefficient # and when the L2(0 < x < 1,y = 1) n o r m o f the function q is used for the regularization, i.e. in discretized form: # ~
q~,
i, such as y = 1
(26)
i
0.5
0.5
t
i
I
I
I
I
i
I
I.E-16
0
0
i
i
i
i
i
i
I.E-15
i
i
u
-0.5
-0.5 -1
-1-2" 5 ~ ~ q ~
-1.5 -2
-2.5
-2.5
- 3
-3 0.2 0.5
I
0.4
t
,
I
0.6 I
,
I
I
I
I
0.5
1.E-15
•
i
,
,
0.4
I
|
1.E-14
0
0 -0.5
,
0.2
0.8
i
,
0.6
I
I
w
0.8
l
I
I
I
U
-0.5
/
-1
-1
-1.5
-1.5
-2
-2
q
-2.5
q
-2.5 -3
-3
o 0.5
o'.2
o'.4
l.E/14.
'
0
,
,
=
,
,
,
a U
0.5
-1
-1
-1.5
-1.5
-2
-2
q '
0',2
'
0'.4
0'.6
'
018
'
x (y=l) Fig. 2. Computed
I
I
I
I
o18 I
I
I
q
-2.5
-3
I
0'.8
0 -0.5
0
I
0'4
I.E-13
-0.5
-2.5
o'.2
o'.e
values of u(x,l) and q(x, I) with the
Euclidean norm for regularization: (a) # = 1076; (b) # = 10 is; (C) # ~- 10 -14.
-3
0
0'.2
' 0'.4
' o'.5
' o18
'
x (y=Z)
Fig. 3. Computed values of u(x,l) and qlx, 1) with the weighted norm for regularization: (a~3# = 10~5; (b) # = 10-14;
(c) ~ = lO
.
Inverse diffusion by boundary elements These results have been obtained with N = 160 (= 4 x 40) boundary elements, which is the most severe situation considered in Ref. 13, and N ' = 41 internal points. They demonstrate the efficiency of the regularization technique, which one permits to obtain an accurate solution for # = 10-15. Nevertheless, one observes that when the coefficient # is increased the variations of q become smoother, but quite false in the vicinity of the comers, where the regularization term tends to vanish the q value. In order to overcome this difficulty, without using the facility of a regularization of higher order (see the remark below), the following weighted norm has been tried: 2 2
i, such as y
# 2-, Pi qi, i
(27)
1
with the coefficients Pi calculated, at the middle of each 0.5
boundary element, from the parabola: p = 4(e + x)(1 + e - x)
0.5
|
1.E- 11
0
0
-0.5
-0.5
-1
-1
-1.5
-1,5
-2
-2 -2.5
q
-3
o'.,
0.5
(28)
where e (_> 0) has been taken as equal to half of the boundary element length. The results obtained when such a regularization term is used are presented in Fig. 3 (a, b, c), for three values of the regularization coefficient. One can observe the improvement of the solution, which is now smoother, except in the vicinity of the two comers where the regularization term acts in a weaker way, since eqn (28) yields a ratio of the maximum and minimum values of coefficients Pi equal to 10.5. In the considered test problem the internal values are supposed to be perfectly known, i.e. only limited by the Fortran double precision. However, in a realistic context these internal values are generally affected by measurement errors. In order to simulate measurements, the internal values have been altered by Gaussian errors.
1.E-13
-2.5
201
z
!
a
|
ole |
z
-3
o;e i
I
|
o12
0.5
z
t
o'.e
0'.4 |
i
i
z
olo |
|
i
I.E-I0
I.E-12 0
0
U
-0.5
-0.5
-1
-1
-1.5
-1.5 -2
-2
q
-2.5 -3
0'4
o'.a
0.5
J
t
t
i
o'e i
i
-3
ole i
J
i
o12
0.5
|
i
o14 t
|
ols i
z
ole i
t
z
! ,E-.09
I.E-II
0
q
-2.5
0
-0.5
-0.5
-1
-1
-1.5
-1.5 -2
-2
q -2.5
-2.5 -3
o
0'.2'
o'.,
' o'.6
olo
x (y=l) Fig. 4. Computed values of u(x, 1) and q(x, 1) with measurement errors such as 2a = 10-°: (a) # = 10-13; (b) -- 10 -12, (c) # = 10TM (weighted norm).
-3
q
o
' o'.2 ' o'.4
o'.s ' o18 '
x (y=l) Fig. 5. Computed values of u(x, 1) and q(x, 1) with measurement errors such as 2or= 10-~: (a) # = 10-ll; (b) # = 10-1°; (c) # = 10 -9 (weighted norm).
202
R. Pasquetti, D. Petit 0.5
for the regularization can yield artificially a very accurate result, by taking the value of coefficient # sufficiently high. In order to obtain a more realistic result, it is preferable to assume the boundary conditions unknown at y = 1, where q ¢ 0. Moreover, we have restricted ourselves to a regularization of zeroth order, since by using a regularization of higher order the same phenomenon occurs. Indeed, the L 2 (0 < x < 1, y = 1) of all the derivatives of q with respect to x vanish, owing to the constant value q - - - 2 of the flux density along the side 0 < x < 1, y = 1.
I
I.E-10
0 -0.5
/
-1 -1.5 -2 -2.5 -3 0 0.5
0.2 I
I
0.4 I
I
0.6 I
l
0.8 I
I
I
I.E-09
3 B O U N D A R Y E L E M E N T A P P R O A C H FOR TRANSIENT INVERSE DIFFUSION
0 -0.5 -1
Although transient inverse problems are much more difficult to solve than steady ones, the approach that we are going to describe for the transient case is very similar to the one used for the steady problems. Such an approach has already been described in Ref. 14. Briefly, let us recall that, from the transient diffusion equation:
-1.5 -2 -2.5
q
-3 0.5
o'2 I
I
0'4 I
I
46 I
i
0;8 I
i
I
1 - Ot u = Au + g
(29)
(/
1 .E-08
0
U
where t stands for time and a for the diffusivity, one can deduce the following BIE:
-0.5 -1
CUM,IF +
-1.5 -2 -2.5
auq* dF dt = to
F
aqu* dI" dt to
F
+ J ; i J ~ agu*df, d t + l n U o u * d f f t
(30)
q
-3 '
012
'
o'.4 ' o'.6
'
018
'
x (y=1)
Fig. 6. Computed values of (.;~1t and q(x 1) with measurement errors such as 2~r u : (a) # = 10 o: (b) # ----- 1 0 - 9 ; (C) p, = 10 -8 (weighted norm). The Gaussian distribution is supposed to be centred and the standard deviation cr is taken as equal to half of the measurement precision. The results obtained with 2or = 10 -6, l0 5 and 10 4 are respectively presented in Figs 4 (a, b, c), 5 (a, b, c) and 6 (a, b, c) for different values of the regularization coefficient and when using the weighted norm defined by eqns (27) and (28). They clearly demonstrate the high sensitivity of the inverse problem solution to the measurement errors, but one observes that a right choice of the regularization coefficient can lead to satisfying results. Remark In the problem introduced in Ref. 13 the boundary conditions are unknown along the unitary side 0 < x < 1, y = 0 (instead of y = 1 in the present study). However, along this side, the exact value of the flux density being equal to q = 0, the use of eqn (26) (or eqn (27))
where u0 is now the initial condition and where to and tF stand for the initial and final time. In such a BIE, the fundamental solution u* is a time and space dependent Green function: u*
(47rat)S/2 exp
H(T)
(31)
where H is the Heaviside function, ~- = tF -- t, S is the space dimension, and r is the distance to point M. The discrete form of eqn (30) is obtained by discretization of the time ( f : time index, 0 _
(32)
F = Gq F ÷ s F
where u g a n d q F are the vectors the components of which are u~ and qi at time tF; H and G are new matrices of dimension (N, N ) , as in eqn (3), and s F is a vector which involves the source term and the initial conditions as well as the u f a n d q f vectors for f < F. Considering now N ' internal points, one obtains an equation similar to eqn (4): u IF ÷ H ' u F
= G,qF
+ s tF
(33)
Inverse diffusion by boundary elements
203
u 'F is the vector the components of which are, at time
4.1 Regularization
tF, the u' ( l < i < N ' ) values of u at the internal points, and I-I', G' are matrices of dimension (N', N). As for the steady case, at this point one has to suppose that the number of equations is greater than the number of unknowns, equal to M, which is problem dependent:
The characteristic parameters of the regularization being numerous, some experience is needed from getting satisfying results. For us:
N+NI>M. Nevertheless, an approach similar to the one described for the steady case would be inefficient, owing to the damping and lagging effects of the transient diffusion operator. Efficient inverse methods make use of 'future time steps': in order to determine what happens at time tF one uses data at times tF+r, 0
H
= G u
u
R
1
Eq]
+ S
q
(34)
R
Lu;+R J
qF+R + S'
(35)
where H and G are block triangular matrices of dimension (N(R + 1), N(R + 1)), H ' and G ' block triangular matrices of dimension (N'(R + 1), N(R + 1)), S is a vector o f dimension N(R+ 1), and S' is a vector of dimension N ' ( R + 1). All of them can be easily identified. Equations (34) and (35) are now the equivalents of eqns (3) and (4) for the steady case, and thus can be handled in a very similar manner in order to solve the inverse problem at time t e. It must be mentioned that the computed values of U F ÷ r a n d qF+r, r ~ O, are simply rejected from the solution of eqns (34) and (35), their definitive values being calculated at time t F+'. For the regularization, from our experience it appears that a space regularization over all the future time steps leads to satisfactory results. This yields through substitution by the regularization matrix R of the steady problem, a block diagonal matrix R constituted with ( R + 1) blocks and such as: R = diag ( R , . . . , R). Two examples of application of the boundary element approach presented in this section can be found in Refs 8 and 14. Moreover, in Ref. 8 a comparison with the superposition method as described in Ref. 1 is given.
4 COMPLEMENTARY REMARKS In this section we will discuss the use of the present boundary element approach and more generally inverse diffusion problems.
• In discrete form, the regularization term usually expresses the Euclidean norm of the flux density q along the part of the boundary where the boundary conditions are unknown. Our numerical tests have clearly shown that the use of q as control variable is more efficient than the use of u. Nevertheless, in order to improve the results, it is sometimes interesting to choose a regularization of higher order. In this case, we use the Euclidean norm of a finite difference approximation of a tangential derivative of q. • The value of the regularization coefficient can be easily adjusted when the exact result is known, but this is not the case for real applications. Without information on the final result, the suitable procedure involves adjusting the value of # in such a way that the gaps between the computed values and the measurements are of the magnitude of the measurement precision. • In the framework of transient problems one has also to adjust the number of future time steps, which depends on the distances from the unknown boundary conditions to the measurement points. From our experience, in realistic applications the parameter R can generally be taken low, e.g. R < 3. • When a sharp corner occurs in the part of the boundary where the boundary conditions are unknown, it is generally difficult to obtain satisfying results in the vicinity of this corner. It is then worthwhile to use the local regularization technique introduced in Ref. 14. It consists of adding to the functional of the optimization problem, an additive local regularization term which expresses the finite difference approximation of one of the tangential derivatives of q.
4.2 Nonlinear diffusion The diffusion equations considered in this paper have been assumed as linear. Concerning nonlinearities in the domain, let us recall here that a nonlinear steady diffusion equation such as, with A for the conductivity: V.(A(u)Vu) + g = 0 may be handled by using the Kirchhoff transformation (see, e.g. Ref. 15). The same remark applies to the transient diffusion equation, but it should be mentioned that in the transient case, the use o f the Kirchhoff transformation can lead to approximations in the results, unless the diffusivity does not depend on u.
204
R. Pasquetti, D. Petit
Concerning nonlinearities at the boundary, one can proceed with the numerical techniques used for direct problems. For us, we utilize an iterative technique, based on local linearizations, i.e. specifics of each boundary element of the nonlinear boundary conditions. ~5 4.3 Characteristic parameters of the diffusion process One of the main difficulties of the inverse problem lies in the fact that the characteristic parameters of the diffusion process, especially the conductivity and the diffusivity, are generally known only approximately. This can lead to serious difficulties, owing to an approximation in these parameters which does not induce centred aleatory errors, as it is usually assumed. The influence of such approximations can be simply checked, by solving the inverse problem with altered values of the characteristic parameters. Thus, in Refs 12 and 14, the inverse problems are solved with a conductivity higher or lower than the conductivity used for the direct problem, and, roughly speaking, this induces a translation in the values of the flux densities, which become higher or lower than the nominal ones. Consequently, when the gaps between the computed values and the measurements are clearly not of the centred Gaussian aleatory type, then one may expect some lack of accuracy in the values of the characteristic parameters.
5 CONCLUSION We have been interested in the inverse problem which involves, for both the steady and transient diffusion equations, the identification of boundary conditions from internal measurements. The boundary element approach which has been suggested is applicable to 2D and 3D geometries. It can yield to space dependent (steady case) or time and space dependent (transient case) boundary conditions, without any assumption on the shape of the solution (constant, linearly variable, etc.). The approach uses regularization techniques, in order to overcome the ill-posed feature of the inverse problem, and, in the transient case, future time steps, owing to the lagging effects of the diffusion operator. Moreover, nonlinear boundary conditions can be supported for the classical superposition method, as well as variable conductivity. For transient problems, the algorithm is a sequential one, which is necessary for large industrial applications. In this paper we have especially focused on: • The quadratic functional, which has to be minimized for solving (at each time step in the transient case) the inverse problem. Three different approaches have been suggested. The first one
(formulation A) is strongly associated with the boundary element formulation, since it expresses the Euclidean norm of the residuals of the BIE. On the other hand, the two others (formulations B and C) result only from statistical considerations, by using the gaps between the computed and the measured values. Formulation B especially uses the functional introduced in Ref. 1 and it may be noticed that the improved functional of formulation C is not accessible to the classical superposition method. The regularization techniques, which are generally necessary to get an accurate solution. Considering the test problem introduced in Ref. 13, for comparing with the Laplace equation different inverse methods, a satisfactory result has been produced by using as regularization term the Euclidean norm of the normal derivative. Moreover, it has been shown that the use of a weighted norm can still improve this result and that the method remains efficient when measurement errors are introduced. Finally, some complementary remarks, taken from our experience of the solution of inverse diffusion problems by boundary elements, have been mentioned.
REFERENCES 1. Beck, J. V., Blackwell, B. & St. Clair, C. R. Inverse Heat Conduction, Ill Posed Problems. Wiley Interscience. New York, 1985. 2. Weber, C. F. Analysis and solution of the ill-posed inverse heat conduction problem. Int. J. Heat Mass Transfer, 1981, 24, 1783-92. 3. Raynaud, M. & Bransier, J. A new finite difference method for the nonlinear inverse heat conduction problem. Num. Heat Transfer, 1986, 9(1), 27-42. 4. Alifanov, O. M. et al. Boundary inverse heat conduction problem in parametric form. In 2nd Int. Conf. on Advanced Computer Methods in Heat Transfer, Milan, 1992, Vol. 1, pp. 479-89. 5. Guerrier, B. & Brnard, C. Identification of time and space dependent boundary conditions in a 2D thermal system. In 2nd Int. Conf. on Advanced Computer Methods in Heat Transfer, Milan, 1992, pp. 448-65.
6. Fort, C. Estimation de paramrtres et m&hode inverse en thermique. Application ~i la drtermination de la variation du flux parirtal dans une chambre de combustion. Thrse de doctorat, University of Poitiers, 1989. 7. Kurpisz, K. & Nowak, A. J. BEM approach to inverse heat conduction problems, Engng Anal. with Bound. Elem., 1992, 10, 291-7. 8. Pasquetti, R. & Le Niliot, C. Boundary element approach for inverse heat conduction problems: Application to a bidimensional transient numerical experiment. Num. Heat Transfer, Part B, 1991, 20, 169-89. 9. Wrobel, L. C. & Brebbia, C. A. The boundary element method for steady-state and transient heat conduction. In Numerical Methods in Thermal Problems, Vol. 1, ed. R. W. Lewis & K. Morgan. Pineridge Press, Swansea, 1979.
Inverse diffusion by boundary elements 10. Thikonov, A. & Arsenine, V. M~thode de Rdsolution des Problbmes Mal Posds. Editions de Moscou, 1976. 11. Eykhhoff, P. System Identification, Parameter and State Estimation. Wiley-Interscience, New York, 1974. 12. Petit, D., Debray, V., Le Niliot, C. & Pasquetti, R. Identification of local heat transfer coefficient using a boundary element formulation: In 2nd Int. Conf. on Advanced Computer Methods in Heat Transfer, Milan, 1992, Vol. 1, pp. 467-77. 13. Ingham, D. B. Improperly posed problems in heat transfer. In Boundary Element Methods in Heat Transfer, ed. L. C. Wrobel & C. A. Brebbia. Computational Mechanics
205
Publications, Elsevier Applied Science, London, 1992, pp. 269-94. 14. Pasquetti, R. & Petit, D. Inverse heat conduction problems with boundary elements: analysis of a comer effect. Engng Anal. with Bound. Elem., 1994, 13, 33-62 15. Pasquetti, R., Caruso, A. & Wrobel, L. C. Transient problems using time dependent fundamental solutions. Boundary Element Methods in Heat Transfer, ed. L. C. Wrobel & C. A. Brebbia. Computational Mechanics Publications, Elsevier Applied Science, London, 1992, pp. 33-62.