- Email: [email protected]
Reconstruction of the velocity distribution in conducting melts from induced magnetic field measurements
-_ __ l!iB
J-J&
ELSEYIER
Computer Physics Communications Computer Physics Communications 86 (1995) 255-263
Reconstruction of the velocity distribution in conducting melts from induced magnetic field measurements D.V. Berkov ‘, N.L. Gorn 2 Inst.fiir Werkstoffiissenschaften-6,
Uniu. Erlangen-Niimberg,
Martensstrasse 7, D-91058 Erlangen, Germany
Received 10 August 1994;revised 20 December 1994
Abstract A possibility to reconstruct the velocity field in a conducting melt without using tracer particles or any mechanical contact to the melt is presented. To recover the flow distribution the liquid under study is placed in an external magnetic field and the induced magnetic field due to the currents appearing in such a liquid is measured. The problem of velocity field reconstruction from these measurement data is closely related to inverse problems of potential theory. It is shown that in contrast to the latter problem the reconstruction of the velocity field can be performed uniquely if there exists a possibility to measure the induced field for different applied field configurations. A stable solution of the corresponding integral equation (arising from a Poisson equation for the induced magnetic field), which is known as an ill-conditioned problem, is obtained by means of the method of singular value analysis. Demonstration results for the case of a two-dimensional flow pattern in an infinitely long pipe are presented.
1. Introduction
Knowledge of the velocity field in a liquid is important in many technical applications. The determination of a flow distribution is a complicated problem, especially in applications where the liquid is not transparent for visible,. IR or UV light. This is unfortunately the case for such important technologies as semiconductor crystal growth or metallurgical processes where liquid metals are involved [Il. Various methods using either tracer particles [2,3] or sensors placed in the melt [4,5] were proposed. The first group requires a very precise choice of the tracer particle density for each liquid and assumes that there exists a method allowing to observe the movement of these particles (X-ray absorption, ultrasound beam frequency shift, etc.). By methods which use probes imbedded in the melt only local information about the liquid velocity in the vicinity of the sensor is available. Both groups face serious additional difficulties if the convection in a hot and (or) aggressive liquid should be studied.
r Permanent address: Inst. of Chem. Physics, Chernogolovka, ’ Permanent
address:
Moscow
Aircraft
OOlO-4655/95/$09.50 0 1995 Elsevier SSDI OOlO-4655(95)0001 l-9
Institute,
Science
Russia.
Russia.
B.V. All rights reserved
D.V. Berkov, N.L. Gorn/ Computer Physics Communications86 (1995)255-263
256
Recently a new method was introduced [6], where neither tracer particles nor direct probes to the liquid are used. In this method the conducting liquid is placed in an external magnetic field which leads to the appearance of electrical currents in the flowing melt. These currents generate a magnetic field, and this induced field is measured outside the fluid volume. In [6] the authors could detect the unsteady flow behavior (non-stationary convection) occurring due to a hydrodynamical instability in the presence of a vertical temperature gradient. Thus experimental evidence was presented that the magnetic field produced by a conducting liquid in a sufficiently weak external magnetic field (so that the flow itself is not disturbed) can be measured with sufficient accuracy. The next step obviously would be an attempt to reconstruct the flow distribution in the melt from the measurement of the induced magnetic field. This question was also discussed in [6], where a certain connection between the flow configuration and the induced magnetic field was established using numerical simulations. Arguments for the possibility of the unique flow recovering were presented, but this possibility was not rigorously proved. In this paper we present a method which allows a complete reconstruction of the velocity field in the conducting liquid from the induced magnetic field outside the melt. In Section 2 the basic integral equation connecting the fluid velocity and the induced magnetic field is derived and a possibility to solve this equation uniquely is analyzed using its relation to inverse problems of potential theory. The regularization method required for the solution of this ill-conditioned integral equation is demonstrated in detail in Section 3 on the simple 1D example of a radius-dependent flow in an infinitely long pipe. Finally a generalization of the method to 2D and 3D problems is presented and demonstrated for a 2D velocity distribution.
2. Integral equation for the velocity distribution The integral equation connecting the velocity distribution u(r) of the conducting melt with the conductivity cr and the induced magnetic field b(r) can be obtained directly from the general magnetohydrodynamic equation for the total magnetic field B [7],
0B c2 - - = [VX [u XB]] + - - A B , 0t 470"
(1)
where B = Bo(r) + b(r, t) is the sum of the external time-independent field Bo(r) and the induced field b(r, t). Outside the magnet creating the external field we have AB 0 = 0 and (1) converts to 0b
c2
Ot = [VX [,, Xao] ] + [Vx [,, xb]] + 4--~ab. If we are interested in studying steady flow, the induced field also becomes time-independent, so that Ob/Ot = 0. The last term on the right-hand side representing the "diffusion" of the induced field is of the order c2Ab/4"rrtr = c2b/4~rtrl 2 = fd, where l denotes the characteristic length scale of the flow distribution. The term describing the "convective" changes in the induced field can be estimated as [Vx [u x b]] ~-ub/1 =fc, so that the ratio of these terms is the so-called magnetic Reynolds number f d f d = 4~rcrul/c 2=- Re m. For values of physical quantities typical for the crystal growth processes (tr= 1016... 1017 s -1, u = 1... 10 cm/s, l = 1... 10 cm), Re m -~ 10-2... 10 -3 << 1, so that the convection term can be neglected and we obtain the final differential equation for the induced field [6], Ab = -
4"rro" c-----T-[WX[u x a 0 ] ].
(2)
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
257
It can also be shown [6], that the same equation (2) is valid for a physical situation, where the so-called magnetic Prandtl number is small: Pm = 4~rtrk/c z << 1 (k denotes the thermal diffusivity of the melt). For melts with a relatively low diffusivity (so that k << ul) this condition may be satisfied for a wider range of applications than the inequality Re m << 1. The solution of the Poisson-like equation (2) is given by the Biot-Savart formula, so that for non-magnetic materials we obtain
~r t [VX [u(r') XBo(r')] ] Ir'-rl
b(r)=~-Jv
dV'.
(3)
To reconstruct the velocity distribution in the conducting liquid from the induced magnetic field we have to solve the integral equation (IE) (3). It is easy to show that if we can perform measurements only for one single configuration of the external field and if the induced field can be measured only outside the fluid volume, this reconstruction can in general not be performed uniquely. Indeed, the solution of (3) is equivalent to the solution of the inverse problem of potential theory [8]: find the charge (mass) distribution density inside a body from the electrostatic (gravitational) potential measured outside it. It is well known, that this problem does not have a unique solution in the general case. For example, a sphere with a radially symmetric charge distribution produces an outside potential fb = Q/r, depending only on its total charge Q, so that no other information about the charge distribution density is available in this case. The mathematical background of this problem can be understood as follows [9]: outside a body the potential ~b is a harmonic function (A~b = 0, in our case Ab = 0). In the 3D case a harmonic function (which tends to zero at the infinity) outside some closed surface is completely defined by its values on this surface (Dirichlet problem). So if we measure the potential on some closed surface outside the charged body, we have all the information about the potential in the whole space, and no other additional information is available. But the values of the potential on the closed surface provide only an amount of 2D information from which a truly 3D charge distribution obviously cannot be recovered uniquely. In our case, however, we can change the configuration of the external magnetic field B0(r), and measuring the induced field b(r) for different B0-configurations we can (at least in principle) obtain the necessary amount of 3D information. For example, if it would be possible to apply the external field only on the given point Bo(r)= B o ~ ( r - ro), we could reconstruct the velocity at this point measuring the induced field somewhere outside. Scanning with this "point" external field through the fluid volume, we could obtain the complete flow picture.
3. Flow reconstruction for a 1D problem
It is impossible, of course, to create a 3-functional magnetic field or even a field which exists only in a small region of free space. To find out how an actual reconstruction can be performed, we proceed with a simplified version of the problem. We consider a pipe which is assumed to be infinitely long in the z-direction with the liquid flowing along it, so that the velocity does not depend on the z-coordinate and can be written as u(r) = uz(x, y) ez. Furthermore we assume that the external magnetic field lies in the x-y plane and does not depend on the z-coordinate also: Bo(r) = Bx(X, y) ex + By(x, y) ey. It is easy to show that in this situation only a z-component of the induced field exists, so that b(r) = bz(x, y) ez. Rewriting the rotation of the vector product in (3) as
[vx [u Xno]] =
(no
. v),,
- (,, . V)no
258
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
(here we have used the Maxwell equation VB = 0 and the incompressibility of the liquid Vu = 0), we obtain for the geometry explained above
" [B Ou~ where l o g ( r - r') 2 represents the potential for the 2D problems and the integration is performed over the cross section of the pipe. The task is to reconstruct the function uz(x, y) from the external field values bz(x, y). One should realize, of course, that in a real physical experiment such a flow would not induce any external field at all, because all electric currents in the liquid in this case would flow perpendicular to the pipe axis. These currents would produce corresponding closing currents in the pipe walls to fulfil the discontinuity condition for the current density. Hence there would be only closed current loops in the x - y plane, each forming a solenoid infinitely extended in the z-direction, which produces no external field. Here we ignore the magnetic field creating by the closing wall currents mentioned above. So the main purpose of these simplest liquid flow examples discussed below is to demonstrate how the missing information necessary for the unique flow recovering can in principle be extracted and optimally used. For comparison with any actual experimental results a complete 3D model should be solved. For our example, analogously to the 3D case, the induced field is completely determined by its values on some closed curve around the pipe, and hence we have at our disposal a 1D information amount. This means, in turn, that with a single external field configuration only 1D information about the flow can be extracted from the measured data. If we know, for example, in advance that the velocity depends on the radial coordinate only, u z = uz(r) (for the circular pipe), then such a flow can be reconstructed uniquely from the induced field measured for a single external field configuration. Below we consider this simplest but important case in detail to demonstrate the method solving integral equations like (3) and (4). For this problem it is convenient to rewrite the IE (4) in cylindrical coordinates,
bro(~p) =
f:dr'
j0f2~rd~b'(B~ cos ~b' + By sin ¢ ' ) l o g [ r 2 + r 2 - 2 r o r '
r' 0--u-U0r'
cos(tb - 4¢)],
(5)
where R denotes the pipe radius and the notation br0 means that measurements of the induced field are performed on a circle with radius r0(> R) around the pipe. After introducing the notation Ou/Or' = g(r') and
Kro(r', q~) = r ' f : ~ d ¢ ' (B x cos ~b'+ B r sin qS')log[r 2 + r ' 2 - 2ror' cos(~b - q~')],
(6)
Eq. (5) takes the form of a Fredholm integral equation of the first kind for the function g(r') with the kernel Kro(r', ~b): br0(~b ) = / : d r '
K~o(r' , 6) g(r').
(7)
It is well known that the solution of this IE is a so-called ill-conditioned problem in the Hadamard sense [10]. This means that small errors in the input data (unavoidable measurement errors in bro(q~)) can lead to arbitrarily large changes in the solution. The mathematical nature of this instability can be explained using eigenvalues Ai and eigenfunction h i of the integral operator/~ appearing in IE (7), which can be written as b =/~g. The eigenfunctions of K, which can be found as the solutions of the IE /?,hi = Aihi, form a complete orthogonal system of basis functions on the given interval (in our case (0, R)). Expanding b(~b) and g(r) in this system as b(q5) = 2CBihi(cb) and g(r) = EGihi(r) and substituting these expansions into the original IE, we obtain
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
259
the relation between the expansion coefficients as B i = AiG i. Hence the coefficients G i required to calculate g ( r ) are G i = Bg/Ai. The eigenvalues of the integral o p e r a t o r / { with a "smooth" kernel (see [11] for more details) decay usually very rapidly so that all Ai except the first few eigenvalues are very small. For this reason even quite small errors in the coefficients B~ due to measurement errors in b(6), being divided by very small A~, lead to large errors in coefficients G i especially for large i, for which the corresponding eigenfunctions are rapidly oscillating. For this reason strong unphysical oscillations occur in the recovered function g ( r ) , which does not allow to extract any meaningful information from the input data. Obviously some smoothing procedure, or regularization of the solution, is required. Among various techniques described in the literature [10-13] we chose the so-called singular value decomposition (SVD) because this method does not contain adjustable parameters in contrast to the methods using some regularization functionals or to the maximum entropy algorithm. Another important advantage of the SVD is the possibility to optimize the experimental setup in advance so that maximum information can be extracted from the measurement data for the given measurement accuracy. The first step in the numerical solution of (7) is the discretization of the problem: applying some numerical quadrature formula to the integral in (7) we obtain the system of linear algebraic equations b =/
(8)
Here b is the M-vector consisting of M induced field values b i ==-b(6g) measured with the experimental errors tr~ (in the following we assume that all values are measured with the same accuracy, i.e., tri = or, but the procedure described below can be extended also on the case of different ~). The N-vector g contains N values g i - g ( r j ) of the function to be recovered in the grid points rj used by the approximation of the integral (7), and K is an M x N matrix obtained from the same integral by numerical quadrature. It is reasonable to have more data values than grid points where the function is to be recovered, so that M > N and (8) can be solved in the least-square sense only: find a vector g which minimizes the deviation X 2 = IIb - Kg II z To obtain this solution we perform the SVD of the system matrix decomposing it as K = U S V T, where U and V are M X M and N x N orthogonal matrices (i.e., U T = U -~) and S is an M x N rectangular matrix which has non-zero elements only on its diagonal beginning in the upper left corner. These elements s i (i = 1, N ) are called singular values of the matrix K and are very useful by analyzing least-square fitting problems [12,13]. This decomposition allows us to rewrite (8) as b = u S V T g . Introducing new vectors /~ = U T b and = V X g , we obtain a diagonal system /1 = S~ with the solution gi = [~i/si • The function to be recovered can be found as g = V~. For an ill-conditioned problem some singular values of K are very small. In fact, they are closely related to the eigenvalues of the integral operator/(" described above: for the symmetric operator, converted to the symmetric matrix K these singular values are simply the corresponding eigenvalues of the matrix. To eliminate the instabilities connected with the small singular values the following procedure is applied [12,13]: we compute g = V~ taking gi = [~i/si for s i > Smin and gi = 0, if the corresponding s i < Smin. The choice of the threshold Smin depends on the concrete problem and is extensively discussed in the literature [8,11-13]. The only reduced quantity we have at our disposal to determine the Smin-value is the relative measurement error tr, so the natural choice would be Smi, ~ ~r. Some additional factors depending on the number of measurement data are discussed in a literature [12,13], but we obtained a good reconstruction quality simply choosing Smin = tr. We also found, that this threshold can be lowered if the number of measurement data is increased; further study of this problem is in progress. Now it is clear that we have to choose the experimental setup so that the system matrix K obtained for our experimental conditions has as large relative singular values (values normalized to Sm,x = max{s~}) as possible. The more sufficiently large singular values are available, the more oscillations in g ( r ) can be
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
260
Table 1 Normalized singular values of the matrix K for various experimental methods (i)
(ii)
(iii)
(iv)
External field
Gaussian
Gaussian
Dipolar
Dipolar
Measurement method
radial
circular
radial
circular
1 2 3 4 5
1.00 7.79 x 10- 3 1.35 x 10- 4 8.18X 10 -6 6.93 X 10 8
1.00 8.48 x 9.89 x 1.94x 3.09 x
1.00 4.55 x 1.96 x 5.53x 3.93 x
10 -1 10- 3 10 -3 10-4
10 -2 10 - 2 10 -4 10-s
1.00 3.01 x 9.34 x 2.94x 9.09 x
10-1 10- 2 10 -2 10-3
reconstructed. In our case the system matrix K depends both on the measurement geometry and on the external field configuration, which allows us to search for the optimal experimental strategy for the flow reconstruction. First of all, for our concrete problem (7) it should be pointed out that no reconstruction with a homogeneous external field is possible. Substituting in the kernel (6), for example, B x = B 0 = const., By = 0, we immediately obtain the result 2~'BoC 2 bz(r ) =
- - , r
where C2 denotes the second moment of the velocity derivative. This means that b z ( r ) does not depend on the polar angle ~b and decays as r-1 no matter which flow distribution exists. Performing the SVD analysis of this experimental situation, we would find out that all singular values of the K-matrix except the first one are zeros. Hence an inhomogeneous field is required. To find out which experimental setup is optimal for the flow reconstruction, we performed an analysis of four experimental situations choosing between Gaussian or dipolar external fields and radial or circular measurements. A Gaussian field was chosen as B~ = - 0 A / a y , By = a A / ~ x , A = A o e x p { - [(x - x 0 ) 2 + y 2 ] / 2 A 2 } with x 0 = 0.5 and A = 0.1; dipolar field means that the field of a 2D magnetic dipole which was located near the pipe wall with the orientation perpendicular to the pipe radius was inserted. Radial measurements stand for the induced field measurements performed along the continuation of the pipe radius (for r > R); circular measurements are performed on a circle around the pipe. Calculations were performed with M = 64 measurement points, and N = 32 grid points for the flow reconstruction. The first five singular values for all strategies are listed in Table 1. Clearly circular measurements for the dipolar external field are superior to the other configurations. And indeed, comparing the flow reconstruction for the last two methods for various flow configurations, we found that the setup (iv) is qualitatively better - see Fig. 1. In both cases the induced field was calculated from the given flow profile via the integral (7); then a random Gaussian noise with the amplitudes listed in the figure caption was added to the computed bz-values, and the flow reconstruction using the SVD algorithm described above was performed.
4. The 2D case and possible generalization of the reconstruction method for 3D problems If the flow distribution in the pipe is not one-dimensional (i.e., u z = u z ( r , ~b)), then such a flow cannot be recovered uniquely from the measurements performed with a single external field configuration for the reasons explained above. Hence we have to measure the induced field for different B0-configura-
D. V. Berkov, 1V.L. Gorn / Computer Physics Communications 86 (1995) 255-263
261
010
II' r
(I')
i
000
-0.10
-0.20
-0.30
r/R -0.40 0.0
0.2
0.4
o3o010020 "'(r) r 0.00
0.6
0.8
1.0
;i( • m~mm
,0.I0
r/R -0.20 00
0.2
04
0.6
0.8
10
Fig. 1. Two examples of the flow reconstruction for different measurement methods. Solid line: velocity derivative u'r to be reconstructed. Line with filled squares: reconstruction from radial measurements. Line with circles: reconstruction from circular measurements. The amplitude of the Gaussian random noise added to the computed values of the induced field was ~r = 1% for (a) and ~r = 0.5% for (b).
tions. The simplest possibility for the dipolar external field is to rotate the magnet creating this field around the pipe and to measure the induced field on a circle around the pipe for each magnet location. For 2D and 3D problems it is convenient to expand the function to be recovered over some suitable complete orthogonal basis is the corresponding functional space. In our problem of a flow in a cylindrical pipe the most natural possibility is to represent the fluid velocity as
uz(r, ~b) = E
E [a~l) sin(k~b) + a ~ ) cos(k~b)] Jk
,
(9)
k=0l=l
where Jk(x) is the Bessel function of order k and fl~t is its lth root. This particular choice of the radial-dependent functions ensures that the boundary condition uz(R, ~b) = 0 is satisfied automatically. Substituting this expansion into the IE (7), we obtain a system of linear algebraic equations connecting the measured values of the induced field with the expansion coefficients a,t, which also should be solved in the least-square sense. The maximal orders of basis functions kmax and lmax should be chosen in agreement with the experimental precision and with the desired resolution. Usually the SVD analysis of the obtained system gives a good guide: the number of basis functions for each spatial dimension should not be much larger than the number of "good" singular values of the system matrix.
D.V. Berkov, N.L. Gorn/ ComputerPhysicsCommunications 86 (1995) 255-263
262
0.4 0.2 o
.~
~'-o.2 -0.4 -0.6 0.5
o
o Y
-1
-1
x
(b)
0'] 0.2
"~'0.2 -0.4
-0.~ 0.5
...._~.i~
1 0.5
0 Y
-1
-1
x
Fig. 2. 2D flow reconstruction. (a) The velocity distribution to be reconstructed. (b) Reconstruction for a Gaussian random noise amplitude o- = 1%.
Representation (9) has several advantages, which are important for a numerical implementation of the algorithm. First of all, in the IE (3) for the general case
r ( B o " V)u - ( u . V ) B o dV b(r) = ~Jv ]r' r]
(10)
both the velocity components and their derivatives are present. Using (9) these derivatives can be expressed analytically via the same basis functions and expansion coefficients a, which makes it possible to treat terms with and without derivatives in the same manner. The analytical representation (9) also makes it easy to adjust the number of variables according to the available amount of information. And in addition, it makes it possible to fulfil the boundary condition or other constraints imposed on the solution in advance. An example of the velocity reconstruction for the 2D case is given in Fig. 2. The induced field was calculated from the given velocity profile
Uz(r, 4)) =
-
~
[sin 4) + cos(2~b)]
(11)
(which is represented in Fig. 2a) for Nm~a~= 12 angular positions of the magnet creating the external dipolar field. If no random noise was added to the induced field, a good reconstruction was achieved proving that a unique flow reconstruction from the measurements of the induced magnetic field can in principle be performed. When a random Gaussian noise with tr = 1% is added to the input data, the
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
263
quantitative agreement becomes poorer, but the main features of the flow are still reproduced quite well (Fig. 2b). The amplitude of the random noise was chosen according to experimental results published recently [6], where it was shown that the m e a s u r e m e n t accuracy within several percents can be achieved. The reconstruction of a flow in the 3D case can be performed in principle in the same way. A possible experimental procedure include m e a s u r e m e n t s of the induced field on the (closed) surface surrounding the vessel with the fluid for several configurations of the external magnetic field to provide the necessary amount of information. However, the expansion (9) in 3D should be performed over a function system which fulfils not only the no-slip boundary conditions on the surface u ( r ~ S) = 0, but also the constraint Vu -- 0 due to the incompressibility of the liquid (for our example this was always true, because the only existing velocity component u z did not depend on the z-coordinate). For this purpose one has to choose a suitable function basis for each concrete geometry, as described, for example, in [14].
5. Conclusion The possibility of a complete reconstruction of the flow velocity distribution from the measurements of the induced magnetic field is shown. The validity of the suggested method is demonstrated on 1D and 2D problems. The optimal experimental procedure among the tested ones involves measurements of the induced field on the circle (2D case) or on the surface (general 3D problems) surrounding the vessel. The magnet creating the external field should be rotated around the vessel, and the complete m e a s u r e m e n t set has to be performed for each magnet position. We found that for a satisfactory flow reconstruction the n u m b e r of measured induced field values for each magnet position should be a few times larger than the n u m b e r of grid points (1D) or basis functions chosen for the velocity field representation (in 2D and 3D) if the relative m e a s u r e m e n t s error is or ~ 1%.
Acknowledgements The authors greatly acknowledge many valuable discussions with Prof. Dr. A. H u b e r t and Dr. S. Meshkov. They thank Dr. J. Baumgartl for clarifying some magnetohydrodynamical questions and Prof. Dr. G. Miiller for supporting this work. Financial support of the Deutsche Forschungsgemeinschaft is gratefully acknowledged.
References [1] G. Miiller, Convection and inhomogenities in crystal growth from the melt, in: Crystals, Vol. 12 (Springer, Berlin-Heidelberg, 1988). 12] K. Kakimoto, M. Eguchi, H. Watanabe and T. Hibija, J. Cryst. Growth 88 (1988) 365. [3] Y. Takeda, Nucl. Techn. 79 (1987) 120. [4] C.B. Reed, B.F. Picologlou, P.V. Dauzvardis and J.L. Bailey, Fusion Techn. 10-3 (1986) 813. [5] S. Horani and L. Krebs, in: Experimental Heat Transfer, Fluid Mechanics and Thermodynamics, eds. R. Shah, E. Ganic and K. Yang (Elsevier, Amsterdam, 1988) p. 279. [6] J. Baumgartl, A. Hubert and G. Miiller, Phys. Fluids A 5 (1993) 3280. [7] P.H. Roberts, An Introduction to Magnetohydrodynamics(Longman, London, 1967). [8] V.G. Romanov, Inverse Problems of Mathematical Physics (VNU, Utrecht, 1987). [9] O.D. Kellogg, Foundation of Potential Theory (Springer, Berlin, 1967). [10] A.N. Tichonov, Solutions of Ill-posed Problems (Winston, Washington, 1977). [11] M. Bertero, C. De Mol and E.R. Pike, Inverse Problems (GB) 1 (1985) 301; 4 (1988) 573. [12] W.H. Press, S.A. Teukolsky, V.H. Vetterling and B.P. Flannery, Numerical Recipes in FORTRAN, 2nd Ed. (Cambridge Univ. Press, Cambridge, 1992). [13] G.E. Forsythe, M.A. Malcolm and C.B. Moler, Computer Methods for Mathematical Computations (Prentice-Hall, Englewood Cliffs, NJ, 1977). [14] P.M. Morse and H. Feshbach, Methods of Theoretical Physics, Part II (McGraw-Hill, New York, 1961).
J-J&
ELSEYIER
Computer Physics Communications Computer Physics Communications 86 (1995) 255-263
Reconstruction of the velocity distribution in conducting melts from induced magnetic field measurements D.V. Berkov ‘, N.L. Gorn 2 Inst.fiir Werkstoffiissenschaften-6,
Uniu. Erlangen-Niimberg,
Martensstrasse 7, D-91058 Erlangen, Germany
Received 10 August 1994;revised 20 December 1994
Abstract A possibility to reconstruct the velocity field in a conducting melt without using tracer particles or any mechanical contact to the melt is presented. To recover the flow distribution the liquid under study is placed in an external magnetic field and the induced magnetic field due to the currents appearing in such a liquid is measured. The problem of velocity field reconstruction from these measurement data is closely related to inverse problems of potential theory. It is shown that in contrast to the latter problem the reconstruction of the velocity field can be performed uniquely if there exists a possibility to measure the induced field for different applied field configurations. A stable solution of the corresponding integral equation (arising from a Poisson equation for the induced magnetic field), which is known as an ill-conditioned problem, is obtained by means of the method of singular value analysis. Demonstration results for the case of a two-dimensional flow pattern in an infinitely long pipe are presented.
1. Introduction
Knowledge of the velocity field in a liquid is important in many technical applications. The determination of a flow distribution is a complicated problem, especially in applications where the liquid is not transparent for visible,. IR or UV light. This is unfortunately the case for such important technologies as semiconductor crystal growth or metallurgical processes where liquid metals are involved [Il. Various methods using either tracer particles [2,3] or sensors placed in the melt [4,5] were proposed. The first group requires a very precise choice of the tracer particle density for each liquid and assumes that there exists a method allowing to observe the movement of these particles (X-ray absorption, ultrasound beam frequency shift, etc.). By methods which use probes imbedded in the melt only local information about the liquid velocity in the vicinity of the sensor is available. Both groups face serious additional difficulties if the convection in a hot and (or) aggressive liquid should be studied.
r Permanent address: Inst. of Chem. Physics, Chernogolovka, ’ Permanent
address:
Moscow
Aircraft
OOlO-4655/95/$09.50 0 1995 Elsevier SSDI OOlO-4655(95)0001 l-9
Institute,
Science
Russia.
Russia.
B.V. All rights reserved
D.V. Berkov, N.L. Gorn/ Computer Physics Communications86 (1995)255-263
256
Recently a new method was introduced [6], where neither tracer particles nor direct probes to the liquid are used. In this method the conducting liquid is placed in an external magnetic field which leads to the appearance of electrical currents in the flowing melt. These currents generate a magnetic field, and this induced field is measured outside the fluid volume. In [6] the authors could detect the unsteady flow behavior (non-stationary convection) occurring due to a hydrodynamical instability in the presence of a vertical temperature gradient. Thus experimental evidence was presented that the magnetic field produced by a conducting liquid in a sufficiently weak external magnetic field (so that the flow itself is not disturbed) can be measured with sufficient accuracy. The next step obviously would be an attempt to reconstruct the flow distribution in the melt from the measurement of the induced magnetic field. This question was also discussed in [6], where a certain connection between the flow configuration and the induced magnetic field was established using numerical simulations. Arguments for the possibility of the unique flow recovering were presented, but this possibility was not rigorously proved. In this paper we present a method which allows a complete reconstruction of the velocity field in the conducting liquid from the induced magnetic field outside the melt. In Section 2 the basic integral equation connecting the fluid velocity and the induced magnetic field is derived and a possibility to solve this equation uniquely is analyzed using its relation to inverse problems of potential theory. The regularization method required for the solution of this ill-conditioned integral equation is demonstrated in detail in Section 3 on the simple 1D example of a radius-dependent flow in an infinitely long pipe. Finally a generalization of the method to 2D and 3D problems is presented and demonstrated for a 2D velocity distribution.
2. Integral equation for the velocity distribution The integral equation connecting the velocity distribution u(r) of the conducting melt with the conductivity cr and the induced magnetic field b(r) can be obtained directly from the general magnetohydrodynamic equation for the total magnetic field B [7],
0B c2 - - = [VX [u XB]] + - - A B , 0t 470"
(1)
where B = Bo(r) + b(r, t) is the sum of the external time-independent field Bo(r) and the induced field b(r, t). Outside the magnet creating the external field we have AB 0 = 0 and (1) converts to 0b
c2
Ot = [VX [,, Xao] ] + [Vx [,, xb]] + 4--~ab. If we are interested in studying steady flow, the induced field also becomes time-independent, so that Ob/Ot = 0. The last term on the right-hand side representing the "diffusion" of the induced field is of the order c2Ab/4"rrtr = c2b/4~rtrl 2 = fd, where l denotes the characteristic length scale of the flow distribution. The term describing the "convective" changes in the induced field can be estimated as [Vx [u x b]] ~-ub/1 =fc, so that the ratio of these terms is the so-called magnetic Reynolds number f d f d = 4~rcrul/c 2=- Re m. For values of physical quantities typical for the crystal growth processes (tr= 1016... 1017 s -1, u = 1... 10 cm/s, l = 1... 10 cm), Re m -~ 10-2... 10 -3 << 1, so that the convection term can be neglected and we obtain the final differential equation for the induced field [6], Ab = -
4"rro" c-----T-[WX[u x a 0 ] ].
(2)
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
257
It can also be shown [6], that the same equation (2) is valid for a physical situation, where the so-called magnetic Prandtl number is small: Pm = 4~rtrk/c z << 1 (k denotes the thermal diffusivity of the melt). For melts with a relatively low diffusivity (so that k << ul) this condition may be satisfied for a wider range of applications than the inequality Re m << 1. The solution of the Poisson-like equation (2) is given by the Biot-Savart formula, so that for non-magnetic materials we obtain
~r t [VX [u(r') XBo(r')] ] Ir'-rl
b(r)=~-Jv
dV'.
(3)
To reconstruct the velocity distribution in the conducting liquid from the induced magnetic field we have to solve the integral equation (IE) (3). It is easy to show that if we can perform measurements only for one single configuration of the external field and if the induced field can be measured only outside the fluid volume, this reconstruction can in general not be performed uniquely. Indeed, the solution of (3) is equivalent to the solution of the inverse problem of potential theory [8]: find the charge (mass) distribution density inside a body from the electrostatic (gravitational) potential measured outside it. It is well known, that this problem does not have a unique solution in the general case. For example, a sphere with a radially symmetric charge distribution produces an outside potential fb = Q/r, depending only on its total charge Q, so that no other information about the charge distribution density is available in this case. The mathematical background of this problem can be understood as follows [9]: outside a body the potential ~b is a harmonic function (A~b = 0, in our case Ab = 0). In the 3D case a harmonic function (which tends to zero at the infinity) outside some closed surface is completely defined by its values on this surface (Dirichlet problem). So if we measure the potential on some closed surface outside the charged body, we have all the information about the potential in the whole space, and no other additional information is available. But the values of the potential on the closed surface provide only an amount of 2D information from which a truly 3D charge distribution obviously cannot be recovered uniquely. In our case, however, we can change the configuration of the external magnetic field B0(r), and measuring the induced field b(r) for different B0-configurations we can (at least in principle) obtain the necessary amount of 3D information. For example, if it would be possible to apply the external field only on the given point Bo(r)= B o ~ ( r - ro), we could reconstruct the velocity at this point measuring the induced field somewhere outside. Scanning with this "point" external field through the fluid volume, we could obtain the complete flow picture.
3. Flow reconstruction for a 1D problem
It is impossible, of course, to create a 3-functional magnetic field or even a field which exists only in a small region of free space. To find out how an actual reconstruction can be performed, we proceed with a simplified version of the problem. We consider a pipe which is assumed to be infinitely long in the z-direction with the liquid flowing along it, so that the velocity does not depend on the z-coordinate and can be written as u(r) = uz(x, y) ez. Furthermore we assume that the external magnetic field lies in the x-y plane and does not depend on the z-coordinate also: Bo(r) = Bx(X, y) ex + By(x, y) ey. It is easy to show that in this situation only a z-component of the induced field exists, so that b(r) = bz(x, y) ez. Rewriting the rotation of the vector product in (3) as
[vx [u Xno]] =
(no
. v),,
- (,, . V)no
258
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
(here we have used the Maxwell equation VB = 0 and the incompressibility of the liquid Vu = 0), we obtain for the geometry explained above
" [B Ou~ where l o g ( r - r') 2 represents the potential for the 2D problems and the integration is performed over the cross section of the pipe. The task is to reconstruct the function uz(x, y) from the external field values bz(x, y). One should realize, of course, that in a real physical experiment such a flow would not induce any external field at all, because all electric currents in the liquid in this case would flow perpendicular to the pipe axis. These currents would produce corresponding closing currents in the pipe walls to fulfil the discontinuity condition for the current density. Hence there would be only closed current loops in the x - y plane, each forming a solenoid infinitely extended in the z-direction, which produces no external field. Here we ignore the magnetic field creating by the closing wall currents mentioned above. So the main purpose of these simplest liquid flow examples discussed below is to demonstrate how the missing information necessary for the unique flow recovering can in principle be extracted and optimally used. For comparison with any actual experimental results a complete 3D model should be solved. For our example, analogously to the 3D case, the induced field is completely determined by its values on some closed curve around the pipe, and hence we have at our disposal a 1D information amount. This means, in turn, that with a single external field configuration only 1D information about the flow can be extracted from the measured data. If we know, for example, in advance that the velocity depends on the radial coordinate only, u z = uz(r) (for the circular pipe), then such a flow can be reconstructed uniquely from the induced field measured for a single external field configuration. Below we consider this simplest but important case in detail to demonstrate the method solving integral equations like (3) and (4). For this problem it is convenient to rewrite the IE (4) in cylindrical coordinates,
bro(~p) =
f:dr'
j0f2~rd~b'(B~ cos ~b' + By sin ¢ ' ) l o g [ r 2 + r 2 - 2 r o r '
r' 0--u-U0r'
cos(tb - 4¢)],
(5)
where R denotes the pipe radius and the notation br0 means that measurements of the induced field are performed on a circle with radius r0(> R) around the pipe. After introducing the notation Ou/Or' = g(r') and
Kro(r', q~) = r ' f : ~ d ¢ ' (B x cos ~b'+ B r sin qS')log[r 2 + r ' 2 - 2ror' cos(~b - q~')],
(6)
Eq. (5) takes the form of a Fredholm integral equation of the first kind for the function g(r') with the kernel Kro(r', ~b): br0(~b ) = / : d r '
K~o(r' , 6) g(r').
(7)
It is well known that the solution of this IE is a so-called ill-conditioned problem in the Hadamard sense [10]. This means that small errors in the input data (unavoidable measurement errors in bro(q~)) can lead to arbitrarily large changes in the solution. The mathematical nature of this instability can be explained using eigenvalues Ai and eigenfunction h i of the integral operator/~ appearing in IE (7), which can be written as b =/~g. The eigenfunctions of K, which can be found as the solutions of the IE /?,hi = Aihi, form a complete orthogonal system of basis functions on the given interval (in our case (0, R)). Expanding b(~b) and g(r) in this system as b(q5) = 2CBihi(cb) and g(r) = EGihi(r) and substituting these expansions into the original IE, we obtain
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
259
the relation between the expansion coefficients as B i = AiG i. Hence the coefficients G i required to calculate g ( r ) are G i = Bg/Ai. The eigenvalues of the integral o p e r a t o r / { with a "smooth" kernel (see [11] for more details) decay usually very rapidly so that all Ai except the first few eigenvalues are very small. For this reason even quite small errors in the coefficients B~ due to measurement errors in b(6), being divided by very small A~, lead to large errors in coefficients G i especially for large i, for which the corresponding eigenfunctions are rapidly oscillating. For this reason strong unphysical oscillations occur in the recovered function g ( r ) , which does not allow to extract any meaningful information from the input data. Obviously some smoothing procedure, or regularization of the solution, is required. Among various techniques described in the literature [10-13] we chose the so-called singular value decomposition (SVD) because this method does not contain adjustable parameters in contrast to the methods using some regularization functionals or to the maximum entropy algorithm. Another important advantage of the SVD is the possibility to optimize the experimental setup in advance so that maximum information can be extracted from the measurement data for the given measurement accuracy. The first step in the numerical solution of (7) is the discretization of the problem: applying some numerical quadrature formula to the integral in (7) we obtain the system of linear algebraic equations b =/
(8)
Here b is the M-vector consisting of M induced field values b i ==-b(6g) measured with the experimental errors tr~ (in the following we assume that all values are measured with the same accuracy, i.e., tri = or, but the procedure described below can be extended also on the case of different ~). The N-vector g contains N values g i - g ( r j ) of the function to be recovered in the grid points rj used by the approximation of the integral (7), and K is an M x N matrix obtained from the same integral by numerical quadrature. It is reasonable to have more data values than grid points where the function is to be recovered, so that M > N and (8) can be solved in the least-square sense only: find a vector g which minimizes the deviation X 2 = IIb - Kg II z To obtain this solution we perform the SVD of the system matrix decomposing it as K = U S V T, where U and V are M X M and N x N orthogonal matrices (i.e., U T = U -~) and S is an M x N rectangular matrix which has non-zero elements only on its diagonal beginning in the upper left corner. These elements s i (i = 1, N ) are called singular values of the matrix K and are very useful by analyzing least-square fitting problems [12,13]. This decomposition allows us to rewrite (8) as b = u S V T g . Introducing new vectors /~ = U T b and = V X g , we obtain a diagonal system /1 = S~ with the solution gi = [~i/si • The function to be recovered can be found as g = V~. For an ill-conditioned problem some singular values of K are very small. In fact, they are closely related to the eigenvalues of the integral operator/(" described above: for the symmetric operator, converted to the symmetric matrix K these singular values are simply the corresponding eigenvalues of the matrix. To eliminate the instabilities connected with the small singular values the following procedure is applied [12,13]: we compute g = V~ taking gi = [~i/si for s i > Smin and gi = 0, if the corresponding s i < Smin. The choice of the threshold Smin depends on the concrete problem and is extensively discussed in the literature [8,11-13]. The only reduced quantity we have at our disposal to determine the Smin-value is the relative measurement error tr, so the natural choice would be Smi, ~ ~r. Some additional factors depending on the number of measurement data are discussed in a literature [12,13], but we obtained a good reconstruction quality simply choosing Smin = tr. We also found, that this threshold can be lowered if the number of measurement data is increased; further study of this problem is in progress. Now it is clear that we have to choose the experimental setup so that the system matrix K obtained for our experimental conditions has as large relative singular values (values normalized to Sm,x = max{s~}) as possible. The more sufficiently large singular values are available, the more oscillations in g ( r ) can be
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
260
Table 1 Normalized singular values of the matrix K for various experimental methods (i)
(ii)
(iii)
(iv)
External field
Gaussian
Gaussian
Dipolar
Dipolar
Measurement method
radial
circular
radial
circular
1 2 3 4 5
1.00 7.79 x 10- 3 1.35 x 10- 4 8.18X 10 -6 6.93 X 10 8
1.00 8.48 x 9.89 x 1.94x 3.09 x
1.00 4.55 x 1.96 x 5.53x 3.93 x
10 -1 10- 3 10 -3 10-4
10 -2 10 - 2 10 -4 10-s
1.00 3.01 x 9.34 x 2.94x 9.09 x
10-1 10- 2 10 -2 10-3
reconstructed. In our case the system matrix K depends both on the measurement geometry and on the external field configuration, which allows us to search for the optimal experimental strategy for the flow reconstruction. First of all, for our concrete problem (7) it should be pointed out that no reconstruction with a homogeneous external field is possible. Substituting in the kernel (6), for example, B x = B 0 = const., By = 0, we immediately obtain the result 2~'BoC 2 bz(r ) =
- - , r
where C2 denotes the second moment of the velocity derivative. This means that b z ( r ) does not depend on the polar angle ~b and decays as r-1 no matter which flow distribution exists. Performing the SVD analysis of this experimental situation, we would find out that all singular values of the K-matrix except the first one are zeros. Hence an inhomogeneous field is required. To find out which experimental setup is optimal for the flow reconstruction, we performed an analysis of four experimental situations choosing between Gaussian or dipolar external fields and radial or circular measurements. A Gaussian field was chosen as B~ = - 0 A / a y , By = a A / ~ x , A = A o e x p { - [(x - x 0 ) 2 + y 2 ] / 2 A 2 } with x 0 = 0.5 and A = 0.1; dipolar field means that the field of a 2D magnetic dipole which was located near the pipe wall with the orientation perpendicular to the pipe radius was inserted. Radial measurements stand for the induced field measurements performed along the continuation of the pipe radius (for r > R); circular measurements are performed on a circle around the pipe. Calculations were performed with M = 64 measurement points, and N = 32 grid points for the flow reconstruction. The first five singular values for all strategies are listed in Table 1. Clearly circular measurements for the dipolar external field are superior to the other configurations. And indeed, comparing the flow reconstruction for the last two methods for various flow configurations, we found that the setup (iv) is qualitatively better - see Fig. 1. In both cases the induced field was calculated from the given flow profile via the integral (7); then a random Gaussian noise with the amplitudes listed in the figure caption was added to the computed bz-values, and the flow reconstruction using the SVD algorithm described above was performed.
4. The 2D case and possible generalization of the reconstruction method for 3D problems If the flow distribution in the pipe is not one-dimensional (i.e., u z = u z ( r , ~b)), then such a flow cannot be recovered uniquely from the measurements performed with a single external field configuration for the reasons explained above. Hence we have to measure the induced field for different B0-configura-
D. V. Berkov, 1V.L. Gorn / Computer Physics Communications 86 (1995) 255-263
261
010
II' r
(I')
i
000
-0.10
-0.20
-0.30
r/R -0.40 0.0
0.2
0.4
o3o010020 "'(r) r 0.00
0.6
0.8
1.0
;i( • m~mm
,0.I0
r/R -0.20 00
0.2
04
0.6
0.8
10
Fig. 1. Two examples of the flow reconstruction for different measurement methods. Solid line: velocity derivative u'r to be reconstructed. Line with filled squares: reconstruction from radial measurements. Line with circles: reconstruction from circular measurements. The amplitude of the Gaussian random noise added to the computed values of the induced field was ~r = 1% for (a) and ~r = 0.5% for (b).
tions. The simplest possibility for the dipolar external field is to rotate the magnet creating this field around the pipe and to measure the induced field on a circle around the pipe for each magnet location. For 2D and 3D problems it is convenient to expand the function to be recovered over some suitable complete orthogonal basis is the corresponding functional space. In our problem of a flow in a cylindrical pipe the most natural possibility is to represent the fluid velocity as
uz(r, ~b) = E
E [a~l) sin(k~b) + a ~ ) cos(k~b)] Jk
,
(9)
k=0l=l
where Jk(x) is the Bessel function of order k and fl~t is its lth root. This particular choice of the radial-dependent functions ensures that the boundary condition uz(R, ~b) = 0 is satisfied automatically. Substituting this expansion into the IE (7), we obtain a system of linear algebraic equations connecting the measured values of the induced field with the expansion coefficients a,t, which also should be solved in the least-square sense. The maximal orders of basis functions kmax and lmax should be chosen in agreement with the experimental precision and with the desired resolution. Usually the SVD analysis of the obtained system gives a good guide: the number of basis functions for each spatial dimension should not be much larger than the number of "good" singular values of the system matrix.
D.V. Berkov, N.L. Gorn/ ComputerPhysicsCommunications 86 (1995) 255-263
262
0.4 0.2 o
.~
~'-o.2 -0.4 -0.6 0.5
o
o Y
-1
-1
x
(b)
0'] 0.2
"~'0.2 -0.4
-0.~ 0.5
...._~.i~
1 0.5
0 Y
-1
-1
x
Fig. 2. 2D flow reconstruction. (a) The velocity distribution to be reconstructed. (b) Reconstruction for a Gaussian random noise amplitude o- = 1%.
Representation (9) has several advantages, which are important for a numerical implementation of the algorithm. First of all, in the IE (3) for the general case
r ( B o " V)u - ( u . V ) B o dV b(r) = ~Jv ]r' r]
(10)
both the velocity components and their derivatives are present. Using (9) these derivatives can be expressed analytically via the same basis functions and expansion coefficients a, which makes it possible to treat terms with and without derivatives in the same manner. The analytical representation (9) also makes it easy to adjust the number of variables according to the available amount of information. And in addition, it makes it possible to fulfil the boundary condition or other constraints imposed on the solution in advance. An example of the velocity reconstruction for the 2D case is given in Fig. 2. The induced field was calculated from the given velocity profile
Uz(r, 4)) =
-
~
[sin 4) + cos(2~b)]
(11)
(which is represented in Fig. 2a) for Nm~a~= 12 angular positions of the magnet creating the external dipolar field. If no random noise was added to the induced field, a good reconstruction was achieved proving that a unique flow reconstruction from the measurements of the induced magnetic field can in principle be performed. When a random Gaussian noise with tr = 1% is added to the input data, the
D.V. Berkov, N.L. Gorn / Computer Physics Communications 86 (1995) 255-263
263
quantitative agreement becomes poorer, but the main features of the flow are still reproduced quite well (Fig. 2b). The amplitude of the random noise was chosen according to experimental results published recently [6], where it was shown that the m e a s u r e m e n t accuracy within several percents can be achieved. The reconstruction of a flow in the 3D case can be performed in principle in the same way. A possible experimental procedure include m e a s u r e m e n t s of the induced field on the (closed) surface surrounding the vessel with the fluid for several configurations of the external magnetic field to provide the necessary amount of information. However, the expansion (9) in 3D should be performed over a function system which fulfils not only the no-slip boundary conditions on the surface u ( r ~ S) = 0, but also the constraint Vu -- 0 due to the incompressibility of the liquid (for our example this was always true, because the only existing velocity component u z did not depend on the z-coordinate). For this purpose one has to choose a suitable function basis for each concrete geometry, as described, for example, in [14].
5. Conclusion The possibility of a complete reconstruction of the flow velocity distribution from the measurements of the induced magnetic field is shown. The validity of the suggested method is demonstrated on 1D and 2D problems. The optimal experimental procedure among the tested ones involves measurements of the induced field on the circle (2D case) or on the surface (general 3D problems) surrounding the vessel. The magnet creating the external field should be rotated around the vessel, and the complete m e a s u r e m e n t set has to be performed for each magnet position. We found that for a satisfactory flow reconstruction the n u m b e r of measured induced field values for each magnet position should be a few times larger than the n u m b e r of grid points (1D) or basis functions chosen for the velocity field representation (in 2D and 3D) if the relative m e a s u r e m e n t s error is or ~ 1%.
Acknowledgements The authors greatly acknowledge many valuable discussions with Prof. Dr. A. H u b e r t and Dr. S. Meshkov. They thank Dr. J. Baumgartl for clarifying some magnetohydrodynamical questions and Prof. Dr. G. Miiller for supporting this work. Financial support of the Deutsche Forschungsgemeinschaft is gratefully acknowledged.
References [1] G. Miiller, Convection and inhomogenities in crystal growth from the melt, in: Crystals, Vol. 12 (Springer, Berlin-Heidelberg, 1988). 12] K. Kakimoto, M. Eguchi, H. Watanabe and T. Hibija, J. Cryst. Growth 88 (1988) 365. [3] Y. Takeda, Nucl. Techn. 79 (1987) 120. [4] C.B. Reed, B.F. Picologlou, P.V. Dauzvardis and J.L. Bailey, Fusion Techn. 10-3 (1986) 813. [5] S. Horani and L. Krebs, in: Experimental Heat Transfer, Fluid Mechanics and Thermodynamics, eds. R. Shah, E. Ganic and K. Yang (Elsevier, Amsterdam, 1988) p. 279. [6] J. Baumgartl, A. Hubert and G. Miiller, Phys. Fluids A 5 (1993) 3280. [7] P.H. Roberts, An Introduction to Magnetohydrodynamics(Longman, London, 1967). [8] V.G. Romanov, Inverse Problems of Mathematical Physics (VNU, Utrecht, 1987). [9] O.D. Kellogg, Foundation of Potential Theory (Springer, Berlin, 1967). [10] A.N. Tichonov, Solutions of Ill-posed Problems (Winston, Washington, 1977). [11] M. Bertero, C. De Mol and E.R. Pike, Inverse Problems (GB) 1 (1985) 301; 4 (1988) 573. [12] W.H. Press, S.A. Teukolsky, V.H. Vetterling and B.P. Flannery, Numerical Recipes in FORTRAN, 2nd Ed. (Cambridge Univ. Press, Cambridge, 1992). [13] G.E. Forsythe, M.A. Malcolm and C.B. Moler, Computer Methods for Mathematical Computations (Prentice-Hall, Englewood Cliffs, NJ, 1977). [14] P.M. Morse and H. Feshbach, Methods of Theoretical Physics, Part II (McGraw-Hill, New York, 1961).