Engineering Analysis with Boundary Elements 16 (1995) 29-33
ELSEVIER
Copyright © 1995 Elsevier Science Limited Printed in Great Britain. All rights reserved 0955-7997/95/$09.50
0955-7997(95)00056-9
Interrelated fundamental solutions for various heterogeneous potential, wave and advective-diffusive problems Richard Paul Shaw Department of Civil Engineering, State University of New York at Buffalo, 212 Ketter Hall, Box 604300, Buffalo, New York, USA
&
G. Steven Gipson School of Civil Engineering, Oklahoma State University, Oklahoma, USA The Green's fundamental solutions for several second order elliptic partial differential equations governing various physical processes are found to be interrelated. These solutions are required for boundary integral equation formulations. Thus a solution for one such process is essentially a solution for several others. This is particularly useful for heterogeneous media where the catalog of existing fundamental solutions is meager. After a discussion of these relationships, a new class of two dimensional fundamental solutions involving axisymmetric material variations is given. Key words: Heterogeneous media, fundamental solutions, boundary element method. in the radial coordinate. These Green's functions in turn lead to related fundamental solutions for the other equations.
INTRODUCTION The use of fundamental ,Green's functions is vital to the formation of boundary integral equation (BIE) representations which have become an extremely useful tool in modern comput~ttional methods, e.g. boundary element methods (BEM). Unfortunately, the number of available fundamental solutions is limited, especially for problems involving heterogeneous media, thereby restricting the usefulness of the BIE/M approach by requiring various approximate techniques such as iteration based on a known homogeneous medium solution. However, fundamental Green's functions for a number of different physical problems for heterogeneous media are themselves related. "['his is not surprising since the standard technique for iheterogeneous media problems is to reduce them to some equivalent homogeneous media problem with a known fundamental solution. Consider the BIE formulation for physical processes described by potential, advective-diffusive and time harmonic wave equations. Fundamental solutions derived for one equation lead to corresponding but not identical fundamental solutions for the others. As a particular example, a Green's fundamental solution for 2D axisymmetric problems can be developed for a physical parameter whiclh is a homogeneous polynomial
GOVERNING EQUATIONS Steady diffusion through a medium with position dependent conductivity or diffusivity, K(F) (case 1), is governed by V . [K(F)V U 1 (r')] = K ( F ) V 2 U 1 (r') + VK(F). V U 1(~) = -Ql(7)
(1)
while the steady advective-diffusive equation (case 2), with advective velocity 17(F) and diffusivity D(F) (using a different symbol for diffusivity here for clarity) is, V . [Z)(7)VV2(~')] - ¢(7). VU2(7) = D(7)V2U2(7) + [VD(7) - ¢(7)]. VU2(r) = -Q2(7)
(2)
with Q(7) as a volume source in both cases. The corresponding time harmonic, e x p ( - i w t ) , wave equation (case 3), would be V . [K(V)VU3(F)] + N(F)U3(F) = -Q3(7) 29
(3)
R. P. Shaw, G. S. Gipson
30
with K(7) and N(7) representing different material properties for different physical processes (including the time harmonic frequency dependence in the term N). There are many physical processes which are governed by these three equation forms; details of these applications will be left to a later, applied paper. If eqn (1) is divided through by K(7), the resulting form is similar to that for eqn (2) divided through by D(7), i.e.
V 2U(7) + {VK(V)/K(7)}. VU(7) = -Q(7)/K(7) (4a)
v 2v(7) + {[VD(7) -
(4b)
The equations are equivalent if
is made, z eqn (8c) may be further reduced to
V2 W3(7) -- {K-I/2[V2{K1/2} q- A(7)] q-- B(7)} W3(7 ) =-Q3(7)K
I7(7) = -V¢(7) = -V[Do In (Keq(7)) ]
(6)
introducing an 'equivalent' Keq. A constant diffusivity and an irrotational velocity will be used for the advective-diffusive case for the remainder of this discussion. A simple dependent variable transformation will lead these three equations to corresponding Helmholtz forms. Take (7)
which eliminates the term V W(7) from the resulting equations in W(7) for all three equations, using the equivalent form K = exp [C/Do] for eqn (2), i.e.
(8C')
These are all in the form of Helmholtz or modified Helmholtz equations with a position dependent index of refraction, (11)
where fl(7) > 0 (< 0) leads to the standard (modified) Helmholtz form and
(12a)
/32 =-Keq'/2V2{KeU2} or V2(Ke~/2) +/32Kelq/2= 0 (12b)
/33
=
- K Uz[Vz{KU2} + A(7)] + B(F) or V2(KI/2) + [/33+ B(7)] KU2 = -A(7)
BOUNDARY INTEGRAL EQUATION FORMULATION There are several different forms for the governing BIE, e.g. if a fundamental solution, G(F, 7o), is based on eqn (1) with a unit point source, the BIE would be
cV(&) = [ Q(Y)G(Y, Yo) dV(r) Jv
+ Is K(7)[G(7, 7o)VV(7)
(8a)
- U(V)VG(F, 70)] • h(V) dS(F)
V 2 W2(7) - {V0. V0/[4D2o] + V2¢/[2Do]} W2(7) = V2 W2(7) -- {K~ql/2v2{Kel/2}} W2(7)
=-Q2(V)KI/2(V)/Do =-Q(F)eqKeql/2(7) (8b)
(8c)
using the identity
[VK. VK]/[4K 2] + [V2K]/[ZK] = K-U2V2{K 1/2}
(9)
(13)
On the other hand, eqn (1 1) with a unit point source replacing Q(V)KU2(7) for a fundamental solution G*(F, Yo), leads to a BIE cW(Fo)
V2W3(7) - {K I/2V2{KI/2} + N/K}W3(7 )
(12c)
The interrelationship between these equations implies that a fundamental solution for one case may be used to generate a fundamental solution for the others, although this generally implies different material variations in the different physical processes.
V2Wl(7)- {K 1/2V2{K1/2}}Wl(7)
= -Q3(V)K-I/2(7)
1/2(7)
(5)
Thus fundamental solutions for the heterogeneous potential problem for specific K(7) will generate corresponding fundamental solutions for the advectivediffusive problem and vice-versa. For the many advective-diffusive problems with constant diffusivity, i.e. D(7) equals Do, the velocity field in eqn (5) would have to be irrotational, i.e. I7 would be derivable as the negative gradient of a potential,
=-QI(V)K 1/2(7)
(10)
fll = -K-1/Zv2{ K1/2} or VZ(K 1/2) -[- fll K1/2 = 0
VK(V)/K(7) = [VD(7)- 17(r)]/D(V)
U(7) = r(7)W(7) ----K(7) -]/2 W(7)
U(Y) = A(V)K U2(7)+ B(F)K(7)
V 2 Wi(r) + fli(7) Wi(7) = -Qi(7)K-U2(7)
• VV(7)
= -Q(V)/D(7)
If the restrictive, but useful, assumption that N and K are related in eqn (8c),
=
Jv G*(V,Fo)K-Ue(V)Q(V) dV(7)
+ Is{a*(7, 7o)Vw(7) - W(V)VG*(F, Fo)} • fi(7) as(F)
(14)
If U(V) is used as the proper dependent variable here, eqn (13) is recovered with G*(r,ro) as
Potential, wave and advective-diffusive problems KI/2(F)KI/2(Fo)G(F, Fo). As usual, h(F) is the outward (from V(F)) normal to S(F) and c is (0, 1/2, 1) for F exterior to V(F), on S(F) at a smooth point, or interior to V(F) respectively. Clearly a solution for W is equivalent to one for U. SOLUTION PROCEDURE Rather than finding W(F), or equivalently U(F), for given K(7) and Q(F) (and N(7) if required), or the converse, it is simpler here to find both K(7) and U(F), after determininl~; G(7, Fo), from a given /3(?'). Although eqns (11) and (12a, b) appear superficially the same, they are not, in the sense that the equation in K(F) could for example be one-dimensional, e.g. a layered medium, while that in G(F, Fo) could be twoor three-dimensional. This approach is fairly general in that the only restrictions that apply are constant diffusivity and irrotational velocity in case 2, and the relationship between N(F) and K(F) in case 3, but in practice only lead to useful fundamental solutions for limited forms of/3(F), two of which are given next. CONSTANT/3 For example, the solution for a constant/3, i.e./30, has been studied for some time, e.g. Gheorghitza 3 and Cheng, 4 for the heterogeneous medium potential and advective-diffusive problems where it requires V 2 { K 1/2} q-/30 KI/2 =: 0
(15)
and leads to a simple equation on G*(F,70)
= 1/(47rR);
/3o > 0
(17a) (17b)
/30 = 0
= exp [-AIo/2R]/(aTrR);
/30
=
-Ao < 0 (17c)
with similar results in 213 with R = I F - 701, G*(F, 70) = (i/4)H(ol)(/31o/2R); = - In (R)/(27r);
= Ko(Alo/2R)/(2rc);
/30 > 0 /30 = 0
K1/2(x) = c 1 cos (/~x) --~ c 2 sin (Ax)
(19a)
while if/3 = - a 2
KU2(x) = el exp (ax) + c2 exp ( - a x )
(19b)
and the solution for/3 = 0 is
Kl/2(x) = ClX + c 2
(19c)
If the second form, eqn (1 9b), involves an infinite domain, it may be appropriate to set el equal to zero, in which case K(x) would simply be c2 exp ( - 2 a x ) . If K(x) were zero at x = 0 in the first form, eqn (19a), Cl would again be set to zero leaving K(x) = c2 sinZ(Ax) = {1 - cos (2Ax)}/2, which models a periodically layered material, although with an unlikely starting value since zero conductivity is unrealistic physically. However, neither x = 0 nor x = 7r/A need not be in the original domain. Case 2 would have a velocity potential related to these equivalent K(F) values through eqn (6) such that (20)
(16)
whose fundamental solution is well known, i.e. in 3D with R = 17- 7ol, G*(V, 70) = exp [i/31o/2R]/(4zrR);
has been solved, i.e. guessing fli has led to a Green's fundamental solution, but for what physical circumstance. The solution for constant/3 requires eqn (15) to be satisfied for K 1/2 which still leaves a large number of possible physical applications. For case 1, K(7) is an actual material parameter, such as diffusivity or conductivity. In that instance, one of the many forms for K(F) allowable for the same Green's function could be a layered medium where K is a function only of a single cartesian coordinate, x; if/3 is -+-A2, then
(b(x) = cDo In [Keq(X)]
vZo* (F, Fo) +/3oG*(F, 70) = - 6 ( 7 - 70)
31
(18a) (18b)
/3o = - A o < 0 (18c)
Case 3 with both/3 and B(F) constant corresponds to the above, but still has an additional degree of freedom in the term A(F).
for either of the forms for K(x) given above. This would lead to either an exponentially varying velocity potential or a periodic velocity potential (and thus velocity in either case since there is only one component of velocity, Vx, defined by -d~(x)/dx). In fact, if Keq(X) were a simple exponential, the corresponding velocity field would be a constant flow in the x direction. Case 3 for the Helmholtz equation directly is discussed at length in the paper of Shaw and Makris, 2 where a solution is given for a layered material using A(x) as a linear function of x plus an infinite sum of decaying exponentials in x. Clearly other solutions of this form are possible. But more general forms are allowed as well; one apparently new family of solutions is described below.
APPLICATIONS FOR C O N S T A N T / 3
SOLUTIONS FOR AXISYMMETRIC MATERIALS WITH A POWER VARIATION
One of the major difficulties with a semi-inverse approach, as used here, is to determine which 'problem'
Consider fl(F) to be - a , 2r 2n-2 where r is the 2D radial coordinate for cases 1 and 2 and n is an integer.
R. P. Shaw, G. S. Gipson
32 Equation (11) in G* becomes
same form for this derivative for either z o = 0 or
Zo #0,
?o26*(7, 70)/0? + (1/r)OG*(7,7o)/0~
OG*/OR = - c 2 / R
+ (1/rZ)O2a*(F, 70)/002 - a]rZna*(F, Fo) = -r26(F - 7o) = - r g 6(F- Fo)
(21)
A change of variable, p = r" and ~b = nO, transforms this equation to O2G*(/5, ~o)/Op 2 + (1/p)OG*(~, ~o)/Op
+ (1/p2)azG*(~,/~)/Ogb 2 - aZG*(~, ~o)/n 2 = - 6 ( F - Fo)/(n2r2o"-2 )
(22)
The coordinate system in this equation may now be translated to have an origin at Po such that the differential operator, which is a Laplacian, may be expressed in terms of the variable P = 1 ~ - Po[ without changing its form. While the right hand side of this equation is no longer a simple delta function, it is zero when P ~ 0 and singular when P : 0. The solution, for P ¢ 0 may be found by separation of variables. The separation constant should be m2/n 2, where m is an integer to have cyclic behavior in 0, leading to solutions involving Im/n(o
P=ff--ffO = ~- nzo-lh;
zn --zon =
(23)
while for z o = 0, the form is simply /5 = kn;
P = R"
(24)
In both cases, this leaves only one term having the appropriate behavior,
G* (P) = cKo (anP/n)
where the requirement of unit flux through a sphere of radius e around R = 0 gives c2 = 1/(270 in all cases. This gives the original Green's function as
G* (F, 70) = (1/ (27r))Ko(an[F n - Fn]/n)
(28)
which, as expected, is not axisymmetric, i.e. cannot be expressed as a simple function of R = IF-Fo[, except for the case, n = 1 where c2 is still 1/27r. A similar derivation exists for/3 > 0 leading to the standard Bessel function forms. For the potential problem (case 1), eqn (12a) leads to a material property
K(r) = [cllo(anrn/n) + czKo(olnrn/n)] 2
(29)
For the example of a perfect fluid in case 2 such that the velocity potential is harmonic, i.e. satisfied a Laplace equation leaving the form of B as
/3 = - V f b . Vdp/aD 2 = _C~nr2 2 n - 2
----
(1/4D2o) d2~b/dr 2
(30) this leads to a potential which is a homogeneous polynomial of degree n,
O(r) = -[402a2r2"/{ (2n - 1)(2n)}]
(31)
Case 3 for the Helmholtz equation has analogous Bessel function solutions, compounded by the presence of B(F) and A(F). If B(F) were zero and A(F) were merged into the source term, Q3(7), the Green's function would be derived as above, except that the condition at infinity would be an outgoing wave form rather than a decaying one. If on the other hand, B(F) were chosen as a2r 2n-2 to cancel the fl term, and K(F), A(F) were axisymmetric, Kl/2(F) would be given for arbitrary A(r),
K1/Z(r) = Ip=r (1/P) Jf==S P'A(p')dp'
[/~..~zo]n __zon
P = nr~-lR
(27)
+ cl In (r) + e2
(32)
which allows for a fairly general form for N(r), e.g. if A were a constant, Ao,
Kl/2(r) = Ao(r2/4) + c~ In (r) + c~
(33a)
U(r) = AoK1/2(r) + a2r2n-2K(r)
(33b)
(25)
such that
CONCLUSION
OG*/OR = -c2K 1(a.P/n)(an/n)nr~ -l for z o ¢ 0 - c2Kl (anP/n)(a./n)nR "-1 for Zo = 0 (26) An expansion of K l for small argument leads to the
While the usefulness of a boundary integral equation formulation in the solution of engineering problems is clearly established, the restriction on this approach when heterogeneous media are involved is also clear. To make the method completely general would require
Potential, wave and advective-diffusive problems knowledge of Green's fundamental solutions for any form of heterogeneity; this is beyond our current abilities. The approach developed here is still restricted, especially in that it is inw,'rse in nature, i.e. the parameter t3 is chosen so as to allow for the solution for the Green's function and then the appropriate heterogeneity is 'recovered' through solution of a similar equation. Such solutions are very limited, but do represent an improvement over no solution at all. One benefit is that any solution here will apply, in slightly modified form, to several different problems. Perhaps even better though, this indicates an approach to other solutions for other heterogeneities in a manner different from but analogous to that described by Clements. 5 As a fringe benefit, it also allows the testing of various approximate procedures, e.g. ray theory for time harmonic waves, against 'known' solutions. In this sense, this is an exploratory paper, not aimed at the specific solution of a specific problem but rather at
33
an 'approach' to such problems. Details of such applications will be forthcoming in subsequent papers. REFERENCES 1. Brebbia, C. A. & Dominguez, J. Boundary Elements: An Introductory Course, 2nd edn. CMP/McGraw-HiI1, New York, 1992. 2. Shaw, R. P. & Makris, N. Green's functions for Helmholtz and Laplace equations in heterogeneous media. Engineering Analysis, 1992, 10(3), 179-83. 3. Gheorghitza, St I. On the plane steady flow of water through inhomogeneous media. First Symp. on Fundamental of Transport Phenomena in Porous Media, IAHR, Haifa, Israel, 1969. 4. Cheng, A. H-D. Darcy flow with variable permeability. Water Resources Research, 1984, 20, 980-4. 5. Clements, D. L. A boundary integral equation method for the numerical solution of a second order elliptic equation with variable coefficients. Journal of Australian Mathematics Society, 1980, 22(B), 218-28.