- Email: [email protected]
A numerical method for a mathematical model in gravitational hydrodynamics
A Numerical Method for a Mathematical Model in Gravitational Hyclrodynamics* U. Bulgarelli Z.N.S.E.A.N. Via di Vallemno 139 00128 Ram, Italy and Z.A.C. CN.R Vi& de1 Policlinico 137 00161 Romq ztazy M. M. Cerimele Z.A.C. CNR Vkzik de1Policlinico 137 00161 Ram, Italy and A. Zaretti
ZXpartimento di Matematka de1 Politecnico di Milano Piazza Z.ema& du Vinci, 32 20133 Milano, Ztaly Transmitted
by S. Riddi
ABSTRACT A mathematical model describing the interaction of two classical gravitational fields, produced respectively by concentrated and distributed masses, is considered. A numerical methodology for the treatment of that model is developed, and a prototype problem is solved.
*Supported by Gruppi Nazionali di Ricerca MPI.
APPLIED
MATHEMATICS
AND COMPUTATION
(D Elsevier Science Publishing Co., Inc., 1988 52 Vanderbilt Ave., New York, NY 10017
27351-264
251
oo983oo3/88/$03.50
U.BUIG4BELLI,M.M.CEFUMELE,ANDA.ZABETTI
252 INTRODUCTION
Gravitational hydrodynamics or fluid dynamics with long range forces is of special interest for astronomers, physicists, and geologists in studying the dynamical behavior of the components of a galaxy: the stars and the interstellar gas. The models are usually but not always [4] formulated in terms of classical continuum mechanics and are described by highly nonlinear equations (e.g. [7, 81); in fact this kind of problem deals with a self-consistent gravitational field generated by the fled added to the pressure forces. Their solution therefore may be more easily obtained by modem computer methods. In comparison with classical hydrodynamics, the equations of fhrid dynamics with long range forces are less well defined, and the physical phenomena described are less well understood. In this paper we are going to develop a numerical method for describing the behavior of a viscous compressible fluid, with long range forces, when another component of the force field exists and is produced by a single moving mass. For a theoretical study of this problem see [9]. The method is obtained by combining the MAC method for the Navier-Stokes equations [lo] with the Courant-Isaacson-Rees method for hyperbolic systems and with several special numerical procedures. The study of the stability of the methods we are going to use, can be found in [l, 111.
GOVERNING
EQUATIONS
The governing equations for the model described above can be derived from the classical principles of conservation of mass and momentum of a fluid, from a state equation, and from the law of motion of the particle, by considering two additional body forces in the momentum equation. The first is the self gravitational field of the fluid, and the second is a gravitational field produced by the moving particle. For the sake of simplicity the gas is supposed to be moving in a bounded fixed region Q c R 2 (whose boundary is denoted by ag), and the ball S to have mass m and center at the point Y(t) = (x(t), y(t)). Its radius 6 > 0 is to be small enough that one may consider each ball nearly pointlike. These equations constitute a system of six differential equations with six unknowns, as follows: i
‘In what folIows we shall always asmme that the equations are written in dimensionless form.
253
Crauitational Hydmdymdcs
The conservation of momentum:
where u = U(X, y, t), u = U(X,y, t) [(x, y) E Q c R’, 0 Q t Q Tl are the components of the velocity q of the fluid, p = p(x, y, t) [(r, y) E D c R2, 0 < t Q T] is the density of the fluid, ~=p(r,y,t)[(r,y)~ncR~,Ogt~T]isthepressweofthefluid, p, ji (assumed constant) are the viscosity coefficients, @=qx,y,t) [(r,y)~R’, Oat
The conservation of mass: Pt +
(PUL+ WV = 0.
(2)
The equations of motion of the moving particle:
$x(t)=-
\k,(x(t),y(t),
t),
(3) de ,eY(O
= - \k,bW¶YW, t>*
The state equation: P = CP7,
(4
where c and y are constants (c > 0, y > l), with y the specific heat ratio (see e.g. PI).
The potential \k(r, y, t), which will be called the continuous potential (see e.g. [6]), is given by
U. BULGABELLI, M. M. CEBIMELE, AND A. ZABEXTI
254
where P is the domain of R2 in which the equations (1) and (2) hold (that is, the region filled by the continuous medium), and (Y a constant greater than zero. The potential (P(r, Y, t), which will be called the particle potential (see e.g. [6]), is given by
Q(x9y.t)
=P/
log\l(r - 1)2+(Y - 9J2dlds,
w(t))
(6)
where Y(t) =(x(t), y(t)) E R2 is the center of the particle S at the time t, B is a constant greater than zero, and r, y are the coordinates of any point of R 2 outside of S. We want here to emphasize that the ftmctions \k(x, Y, t), @(x, y, t) are defined in the whole space R 2, up to an additive constant; hence the functions x(t), y(t), given by the system (3), can be the coordinates of any point of R 2, while the equations (l), (2), and (4) hold in 0 only, which has been assumed bounded, but arbitrarily large. In solving (l), (2), and (3), the following initial conditions are employed:
and
&o) =x0 x’(t,) = Xl
YOO) = Yo,
(8) YVO) =Yl-
The boundary conditions for (1) are given by
q(r, y, t)-ntx, Y)=07
where n is the outward unit vector normal to a Q, and I, is the unit tangential vector. The equation (2) does not need boundary conditions, because of its hyperbolic behavior (see e.g. [3, 91). We notice that, because of (4), the remaining unknowns of the problem are u(r, II, t), 4x, Y, t), P(X, Y, t), 4th
y(t).
255
Grauitational Hydmdynamics
THE COMPUTATIONAL, METHOD Our finite approximations for the equations (l), (2), (3), (4) have been obtained as in [l], taking into account (5), (6), (7), (8), (9). For computing purposes only, R’ is replaced by a bounded domain U as large as possible, taking into account the available computer resources. This domain, strictly containing Q U iJQ, has been discretized by a staggered grid whose steps are Ax and Ay. We denote by 8, the set of the interiors cells of 51, and by a h, the set of the boundary cells [l]. The components u, u of the velocity and the density p are computed at the following points for each interior and boundary cell: (a) the horizontal component u of the velocity is defined at the center of each vertical bounding segment; (b) the vertical component u is defined at the center of each horizontal bounding segment; (c) the density p is defined at the center of each cell. Finally the functions \k(x, y, t) and @(x, y, t) are computed at the vertices of each cell of the computational mesh. Assuming the time step At > 0, and the vahzs uIfr, utI, & known at the TIS~ t, = nntf + to (n integer with 1 d n Q [T/At]), the values ?I+1 ui ‘“Lf ,P~,~ at the time L+r, in the cell (i, j) are computed by the foiiowing finite difference equations, which are consistent with the equations (l), (2):
-
iw,,+~““>”+(5”)“+1+(T”)n+1 w”
1 f*I
,
+
a”
+
II
w
I
n
c.5
-Id
-
--
n
d
v
0”
14- Idn
II
eel (3
1 @)x,+ ‘“y,)” - --4”(0;.,)+ (AZ)8 %
p’“+‘=;(l+
E;')(p;;')
The field of values ‘p;;‘, P+’ will be obtained as indicated in the following. The finite difference ‘*&uation for the conservation of mass is simply given by
where
258
U. BULCABElLI, M. M. CEBIMELE, AND A. ZABETX’X
We notice that the continuous potential *(x, y, t) given by (5) is used, both in the equations (1) which are defined in 8, and in the equations (3) which hold in the whole space R 2; therefore we have to compute the right hand side of (5) in the whole U. In the numerical integration which has to be performed in order to compute the function \k(x, y, t), the points at which this function is defined are always different from the points at which the function p(&v, t) is evaluated; hence the following formula can be used:
* “+‘=k((i-l)Ax,(j-l)Ay,(n+l)At) ‘*j
h,k,i,
jEZ,
(12)
Since the particle potential @(x, y, t ) is used only in the equations (l), which hold in B only, we confine ourselves to computing this function for (x, y) E Q, wherever the center Y(t) of S, whose coordinates are (x(t), y(t)), is located in U. To this purpose we introduce a system of polar coordinates in U of the form
~(t)=x(t)+rcos8,
0<8<27r,
O
?J(t)=y(t)+rsid,
Og7g6,
where 6 +Z min(Ax, Ay); the last assumption means that the radius of the ball S has to be much smaller than the discretization steps in order to have a reasonable approximation of the physical problem. Hence (6) can be written as follows:
+ [y -y(t)
- rsid]“}
&de.
(14
CravitationaI Hy&odytultnics
259
When we perform the discretization of (14), we find that its right hand side becomes singular when S is in h, U afit,, and the ball contains one of the vertices of a computational cell. To handle this singularity we proceed in the following way. If r(S), y(i) are the coordinates of the center of S at the tune S, the vertices of the cell to which the center belongs are
Q,, = [ik Q2=
jAyI,
[iAr,(j+l)Ay],
Ql=
[(i+l)kjAy],
Q3= [(i+l)Ax,(j+l)Ay],
where i = [x(i)/Ar], j = [y(t)/Ay]. The distance between the center of S and each of these points Q, (p = 0,1,2,3), is
d,=
[(4u- xp.)2+ (Y(i) - YQp)2]1'2*
Hence: if 0 < d, Q S, the point QP belongs to S; if d, > S, the point Q,, belongs to (k?, U ah,) - S; if d p = 0, the center of S falls in the point Q,. To avoid this singularity, it is enough to select the step AT, for the numerical integration in (14), in such a way that d p is not a multiple of Ar when 0 c d, Q 8; finally, when d, = 0 we assume as the value of potential in Q, the mean of the potential computed at the closest four points. By these remarks, a discrete approximation of (14) is given as follows: 4) ““=~((i-1)Ax,(j-1)Ay,(n+1)At) i.I
=$
-
; (Z-1)Ar m-l
+[(k-1)Ay-y(t,,+,)-(Z-l)Arsin(m-~)At3)]2)ArAf3
(15)
with nr = [S/Ar] and n, = [29r/A8]. To describe the particle motion, that is, to compute the two terms x( t,+ l),y( t,+ 1), we have to integrate the system (3) with the initial condi-
U. BUXABELLI, M. M. CEBIMELE, AND A. ZABE’ITI
260
tions (8). First of all we observe that the system (3) can be written as
with
dto) = x0*
40) = x1,
4to) = Yo,
4to) = Yl.
To integrate the system (16), since the radius S of the ball is much smaller than the discretization steps, we prefer to use a method [ll] giving a higher accuracy than the one chosen to integrate the NavierStokes equations. More precisely, we have
(2p+l=
zp+ At’(Z - BAt’J) -lfp)n+l,
07)
where:
(q+‘, zg+l, g+l,
zp+‘=
p = (zzp,-
\k,p, zf;, -
e+:,1), tP
=
z,p+ ‘y, qy,
O
pAt’+ t,,,
0 Q p < [At/At’],
and, because of the linearization used, ‘k,p= q+l(zp,
ZH),
‘k;l’&;+l(zp,
21).
Moreover we have
I
(1 - 0At’J) -’ =
‘1
BAt’
0
0
0 0
1 0
0 1
0 Bbt’
_o
0
0
1
261
Crauituttonul Hydmdynwn~s Therefore (17) becomes z1p+l= zp + At+&‘-
BAt’\k,P),
=sp+l = zc - At’*:, @+I = z3p+ At’( z: - BAt’\k,P),
(18)
,$+I = z; - At”@;. Note that the right hand side of (18) is computed at the point where the ball center is located at the time t,, which usuallyis not a computational point; hence we use a bilinear approximation using the four closest grid points.
REMARK. We recall that the stability condition for the equations (lo), (ll), which can be obtained only in the uncoupled linearized case [l], requires that the time step satisfy the more restrictive condition between the following two:
luig + WA -Ax
+2p
Ay
(19)
while for the system (18), which is unconditionally stable, we do not need any restriction.
SUMMARY OF COMPUTATIONAL
STEPS
Assuming the values &, utj, oIfr tobeknownateachpointof atthetimet,= n At + t,,, the computation proceeds in this way: 1.
pt; ’ is computed
&UC?&
in the whole b, by using (11) and the initial condi-
tiOKL%
2. 3.
3”+ ’ is computed in the whole U by using (12) and the values of pt; ’ ccZputed at step 1. The coordinates of the center of S in U, at the time t,+ 1, are computed by using (18).
282
U.BULCAFiELLI,M.M.CEBIMELE,ANDA.ZABElTI
4. aJ;;1iscomputed in the whole h, u ~36, by using (15) in which the 5.
coordinates of the center of S just computed at step 3 are used. ulf;i, $1 are computed in the whole 8, U ah, by using (lo), the initial boundary conditions, and the values py:l, 9$yf1, Qfff i given in the computational steps 1,2,4.
Then because ah variables are computed at the time t,,+ i, the method starts again from step 1. The FCHWRWprogram used is available in Appendix A of the technical report [2]. DISCUSSION OF NUMERICAL
RESULTS
In this section we are going to integrate the equations (l), (2), (3), (4), (5), (6) as a suitable example, by using the numerical procedure discussed above. More precisely, we suppose that the domains Q and U respectively are defined by Q= {10.~r<21.,10.dy<21.}, U= {O.
\\
/a
FIG. 1. Density contours, t = 4.0. Density (a) 0.1700, (b) 0.1500, (c) O.OWO,(d) 0.0500.
//
I
rf
FIG.2. Density contours, t = 8.0. Density (a) 1.8oo0, (b) 0.1700, (c) 0.0700, (d) 0.0300, (e) 0.0100, (f) 0.@080.
Gmuitational Hydrudms
FIG.3
Particle trajectoxies, t = 4.0.
FIG.4
Particle trajectories, t = 8.0.
located at the point {x(&J = 18, y( to) = 14.) with initial velocity { x’(ts) = 0.1935, y’(t,) = O.}. The fluid filling Q has initial constant density p0 = 0.1, and it is supposed to be initiahy at rest. Because of the initial conditions of the flmd, only the first step of the computation procedure is obtained by the numerical integration of the system (3) by using (5). Figures 1 and 2 show, at the times t = 4.0 and t = 8.0 respectively, the contour lines of the density p of the fluid for the following values: p = {0.170,0.150,0.090,0.050}
(t = 4.0),
p = { 1.800,0.170,0.070,0.030,0.010,0.008}
(t =S.O).
Figures 3 and 4 respectively show, at the same times t = 4.0 and t = 8.0, the trajectory of the particle S starting from the initial point and the fluid velocity field in which S is moving. We notice that in Figures 1 and 2 the density is increasing from the boundary to the center of 9, as can be easily deduced from the equations; moreover, in Figures 3 and 4 we can observe that the behavior of the velocity field and the path of S are consistent with the proposed model. Finally we would like to underline that the described method gives a promising tool for the study of the interactions between two classical gravitational fields. The extension of this method to a model taking into account also the energy equation and the interaction of several particles seems to be quite straightforward.
We wish to erpress our gmtitude to Giuuanni Pause of the Dipartimento di Mat&mat&a &l Politemico di Milano jin- his encoumgement and a useful
U. BULGABEILI,
264 suggestion jbr
hf. M. CEFUMEXE, AND A. ZABETTI
the development of this mod+d, and to Marcella Prosped of
Z.A.C. Rome, who supjdid
the jyaphic mutine.~.
REFERENCES 1
U. Bulgamlh, V. Casulli, and D. Greenspan, hssure Metti for the Num.ericuZ Solution of Free !&j&cc Fhdd Flows, Pineridge Press, Swansea, U.K., 1984.
2
U. BulgamIi, M. M. Cerimele, and A. Zaretti, A Numerical Method for a Mathematical Model Describing Interaction Between Two Gravitational Fields, Quad. IAC No. 16,1985. R. Courant and K. 0. Friedrh, SupersonicFh and Shock Waves, Interscience, New York, 1948. D. Greenspan, Disc&e Numerical Methodi in PhtrJics and Engineering, Academic, 1974. F. B. Hildebrandt, Finite Diffbmnce Equutiom and Simulutions,Prenti~Hall, 1988. 0. D. Kelkqp, Foundation of Potential Theog, Springer, Berlin, 1929. C. C. Lin, Stellar dynamical theory of normal spirals,in Lecturea in Appl. Math., Vol. 9 (G. Ehlers, Ed.), Amer. Math.!%c.,Providence, 1987. D. Potter, ChmputaiionulPh@cs, Wiley, New York, 1973. G. Prouse, F. Bohii, and A. Zaretti, On a mathematical model describing interaction between two classicalgravitatiod fields, Rend. Accad. Nuz. Sci. dei XL Mem. Mat. 11, to appear. J. E. Welch, F. H. H&w, J. P. Shannon,and B. J. Daly, The MAC Method: a computing technique for solving viscous, incompressible, transient fluid-flow problems involving frea surfaces, LA 3425, Los Alama Scientific Labs., Los Alamas, N.M., 1988. D. Trigiante, Asymptotic stability and discretization on an iufinite interval, Computing l&117-119 (1977).
10
11
ZXpartimento di Matematka de1 Politecnico di Milano Piazza Z.ema& du Vinci, 32 20133 Milano, Ztaly Transmitted
by S. Riddi
ABSTRACT A mathematical model describing the interaction of two classical gravitational fields, produced respectively by concentrated and distributed masses, is considered. A numerical methodology for the treatment of that model is developed, and a prototype problem is solved.
*Supported by Gruppi Nazionali di Ricerca MPI.
APPLIED
MATHEMATICS
AND COMPUTATION
(D Elsevier Science Publishing Co., Inc., 1988 52 Vanderbilt Ave., New York, NY 10017
27351-264
251
oo983oo3/88/$03.50
U.BUIG4BELLI,M.M.CEFUMELE,ANDA.ZABETTI
252 INTRODUCTION
Gravitational hydrodynamics or fluid dynamics with long range forces is of special interest for astronomers, physicists, and geologists in studying the dynamical behavior of the components of a galaxy: the stars and the interstellar gas. The models are usually but not always [4] formulated in terms of classical continuum mechanics and are described by highly nonlinear equations (e.g. [7, 81); in fact this kind of problem deals with a self-consistent gravitational field generated by the fled added to the pressure forces. Their solution therefore may be more easily obtained by modem computer methods. In comparison with classical hydrodynamics, the equations of fhrid dynamics with long range forces are less well defined, and the physical phenomena described are less well understood. In this paper we are going to develop a numerical method for describing the behavior of a viscous compressible fluid, with long range forces, when another component of the force field exists and is produced by a single moving mass. For a theoretical study of this problem see [9]. The method is obtained by combining the MAC method for the Navier-Stokes equations [lo] with the Courant-Isaacson-Rees method for hyperbolic systems and with several special numerical procedures. The study of the stability of the methods we are going to use, can be found in [l, 111.
GOVERNING
EQUATIONS
The governing equations for the model described above can be derived from the classical principles of conservation of mass and momentum of a fluid, from a state equation, and from the law of motion of the particle, by considering two additional body forces in the momentum equation. The first is the self gravitational field of the fluid, and the second is a gravitational field produced by the moving particle. For the sake of simplicity the gas is supposed to be moving in a bounded fixed region Q c R 2 (whose boundary is denoted by ag), and the ball S to have mass m and center at the point Y(t) = (x(t), y(t)). Its radius 6 > 0 is to be small enough that one may consider each ball nearly pointlike. These equations constitute a system of six differential equations with six unknowns, as follows: i
‘In what folIows we shall always asmme that the equations are written in dimensionless form.
253
Crauitational Hydmdymdcs
The conservation of momentum:
where u = U(X, y, t), u = U(X,y, t) [(x, y) E Q c R’, 0 Q t Q Tl are the components of the velocity q of the fluid, p = p(x, y, t) [(r, y) E D c R2, 0 < t Q T] is the density of the fluid, ~=p(r,y,t)[(r,y)~ncR~,Ogt~T]isthepressweofthefluid, p, ji (assumed constant) are the viscosity coefficients, @=qx,y,t) [(r,y)~R’, Oat
The conservation of mass: Pt +
(PUL+ WV = 0.
(2)
The equations of motion of the moving particle:
$x(t)=-
\k,(x(t),y(t),
t),
(3) de ,eY(O
= - \k,bW¶YW, t>*
The state equation: P = CP7,
(4
where c and y are constants (c > 0, y > l), with y the specific heat ratio (see e.g. PI).
The potential \k(r, y, t), which will be called the continuous potential (see e.g. [6]), is given by
U. BULGABELLI, M. M. CEBIMELE, AND A. ZABEXTI
254
where P is the domain of R2 in which the equations (1) and (2) hold (that is, the region filled by the continuous medium), and (Y a constant greater than zero. The potential (P(r, Y, t), which will be called the particle potential (see e.g. [6]), is given by
Q(x9y.t)
=P/
log\l(r - 1)2+(Y - 9J2dlds,
w(t))
(6)
where Y(t) =(x(t), y(t)) E R2 is the center of the particle S at the time t, B is a constant greater than zero, and r, y are the coordinates of any point of R 2 outside of S. We want here to emphasize that the ftmctions \k(x, Y, t), @(x, y, t) are defined in the whole space R 2, up to an additive constant; hence the functions x(t), y(t), given by the system (3), can be the coordinates of any point of R 2, while the equations (l), (2), and (4) hold in 0 only, which has been assumed bounded, but arbitrarily large. In solving (l), (2), and (3), the following initial conditions are employed:
and
&o) =x0 x’(t,) = Xl
YOO) = Yo,
(8) YVO) =Yl-
The boundary conditions for (1) are given by
q(r, y, t)-ntx, Y)=07
where n is the outward unit vector normal to a Q, and I, is the unit tangential vector. The equation (2) does not need boundary conditions, because of its hyperbolic behavior (see e.g. [3, 91). We notice that, because of (4), the remaining unknowns of the problem are u(r, II, t), 4x, Y, t), P(X, Y, t), 4th
y(t).
255
Grauitational Hydmdynamics
THE COMPUTATIONAL, METHOD Our finite approximations for the equations (l), (2), (3), (4) have been obtained as in [l], taking into account (5), (6), (7), (8), (9). For computing purposes only, R’ is replaced by a bounded domain U as large as possible, taking into account the available computer resources. This domain, strictly containing Q U iJQ, has been discretized by a staggered grid whose steps are Ax and Ay. We denote by 8, the set of the interiors cells of 51, and by a h, the set of the boundary cells [l]. The components u, u of the velocity and the density p are computed at the following points for each interior and boundary cell: (a) the horizontal component u of the velocity is defined at the center of each vertical bounding segment; (b) the vertical component u is defined at the center of each horizontal bounding segment; (c) the density p is defined at the center of each cell. Finally the functions \k(x, y, t) and @(x, y, t) are computed at the vertices of each cell of the computational mesh. Assuming the time step At > 0, and the vahzs uIfr, utI, & known at the TIS~ t, = nntf + to (n integer with 1 d n Q [T/At]), the values ?I+1 ui ‘“Lf ,P~,~ at the time L+r, in the cell (i, j) are computed by the foiiowing finite difference equations, which are consistent with the equations (l), (2):
-
iw,,+~““>”+(5”)“+1+(T”)n+1 w”
1 f*I
,
+
a”
+
II
w
I
n
c.5
-Id
-
--
n
d
v
0”
14- Idn
II
eel (3
1 @)x,+ ‘“y,)” - --4”(0;.,)+ (AZ)8 %
p’“+‘=;(l+
E;')(p;;')
The field of values ‘p;;‘, P+’ will be obtained as indicated in the following. The finite difference ‘*&uation for the conservation of mass is simply given by
where
258
U. BULCABElLI, M. M. CEBIMELE, AND A. ZABETX’X
We notice that the continuous potential *(x, y, t) given by (5) is used, both in the equations (1) which are defined in 8, and in the equations (3) which hold in the whole space R 2; therefore we have to compute the right hand side of (5) in the whole U. In the numerical integration which has to be performed in order to compute the function \k(x, y, t), the points at which this function is defined are always different from the points at which the function p(&v, t) is evaluated; hence the following formula can be used:
* “+‘=k((i-l)Ax,(j-l)Ay,(n+l)At) ‘*j
h,k,i,
jEZ,
(12)
Since the particle potential @(x, y, t ) is used only in the equations (l), which hold in B only, we confine ourselves to computing this function for (x, y) E Q, wherever the center Y(t) of S, whose coordinates are (x(t), y(t)), is located in U. To this purpose we introduce a system of polar coordinates in U of the form
~(t)=x(t)+rcos8,
0<8<27r,
O
?J(t)=y(t)+rsid,
Og7g6,
where 6 +Z min(Ax, Ay); the last assumption means that the radius of the ball S has to be much smaller than the discretization steps in order to have a reasonable approximation of the physical problem. Hence (6) can be written as follows:
+ [y -y(t)
- rsid]“}
&de.
(14
CravitationaI Hy&odytultnics
259
When we perform the discretization of (14), we find that its right hand side becomes singular when S is in h, U afit,, and the ball contains one of the vertices of a computational cell. To handle this singularity we proceed in the following way. If r(S), y(i) are the coordinates of the center of S at the tune S, the vertices of the cell to which the center belongs are
Q,, = [ik Q2=
jAyI,
[iAr,(j+l)Ay],
Ql=
[(i+l)kjAy],
Q3= [(i+l)Ax,(j+l)Ay],
where i = [x(i)/Ar], j = [y(t)/Ay]. The distance between the center of S and each of these points Q, (p = 0,1,2,3), is
d,=
[(4u- xp.)2+ (Y(i) - YQp)2]1'2*
Hence: if 0 < d, Q S, the point QP belongs to S; if d, > S, the point Q,, belongs to (k?, U ah,) - S; if d p = 0, the center of S falls in the point Q,. To avoid this singularity, it is enough to select the step AT, for the numerical integration in (14), in such a way that d p is not a multiple of Ar when 0 c d, Q 8; finally, when d, = 0 we assume as the value of potential in Q, the mean of the potential computed at the closest four points. By these remarks, a discrete approximation of (14) is given as follows: 4) ““=~((i-1)Ax,(j-1)Ay,(n+1)At) i.I
=$
-
; (Z-1)Ar m-l
+[(k-1)Ay-y(t,,+,)-(Z-l)Arsin(m-~)At3)]2)ArAf3
(15)
with nr = [S/Ar] and n, = [29r/A8]. To describe the particle motion, that is, to compute the two terms x( t,+ l),y( t,+ 1), we have to integrate the system (3) with the initial condi-
U. BUXABELLI, M. M. CEBIMELE, AND A. ZABE’ITI
260
tions (8). First of all we observe that the system (3) can be written as
with
dto) = x0*
40) = x1,
4to) = Yo,
4to) = Yl.
To integrate the system (16), since the radius S of the ball is much smaller than the discretization steps, we prefer to use a method [ll] giving a higher accuracy than the one chosen to integrate the NavierStokes equations. More precisely, we have
(2p+l=
zp+ At’(Z - BAt’J) -lfp)n+l,
07)
where:
(q+‘, zg+l, g+l,
zp+‘=
p = (zzp,-
\k,p, zf;, -
e+:,1), tP
=
z,p+ ‘y, qy,
O
pAt’+ t,,,
0 Q p < [At/At’],
and, because of the linearization used, ‘k,p= q+l(zp,
ZH),
‘k;l’&;+l(zp,
21).
Moreover we have
I
(1 - 0At’J) -’ =
‘1
BAt’
0
0
0 0
1 0
0 1
0 Bbt’
_o
0
0
1
261
Crauituttonul Hydmdynwn~s Therefore (17) becomes z1p+l= zp + At+&‘-
BAt’\k,P),
=sp+l = zc - At’*:, @+I = z3p+ At’( z: - BAt’\k,P),
(18)
,$+I = z; - At”@;. Note that the right hand side of (18) is computed at the point where the ball center is located at the time t,, which usuallyis not a computational point; hence we use a bilinear approximation using the four closest grid points.
REMARK. We recall that the stability condition for the equations (lo), (ll), which can be obtained only in the uncoupled linearized case [l], requires that the time step satisfy the more restrictive condition between the following two:
luig + WA -Ax
+2p
Ay
(19)
while for the system (18), which is unconditionally stable, we do not need any restriction.
SUMMARY OF COMPUTATIONAL
STEPS
Assuming the values &, utj, oIfr tobeknownateachpointof atthetimet,= n At + t,,, the computation proceeds in this way: 1.
pt; ’ is computed
&UC?&
in the whole b, by using (11) and the initial condi-
tiOKL%
2. 3.
3”+ ’ is computed in the whole U by using (12) and the values of pt; ’ ccZputed at step 1. The coordinates of the center of S in U, at the time t,+ 1, are computed by using (18).
282
U.BULCAFiELLI,M.M.CEBIMELE,ANDA.ZABElTI
4. aJ;;1iscomputed in the whole h, u ~36, by using (15) in which the 5.
coordinates of the center of S just computed at step 3 are used. ulf;i, $1 are computed in the whole 8, U ah, by using (lo), the initial boundary conditions, and the values py:l, 9$yf1, Qfff i given in the computational steps 1,2,4.
Then because ah variables are computed at the time t,,+ i, the method starts again from step 1. The FCHWRWprogram used is available in Appendix A of the technical report [2]. DISCUSSION OF NUMERICAL
RESULTS
In this section we are going to integrate the equations (l), (2), (3), (4), (5), (6) as a suitable example, by using the numerical procedure discussed above. More precisely, we suppose that the domains Q and U respectively are defined by Q= {10.~r<21.,10.dy<21.}, U= {O.
\\
/a
FIG. 1. Density contours, t = 4.0. Density (a) 0.1700, (b) 0.1500, (c) O.OWO,(d) 0.0500.
//
I
rf
FIG.2. Density contours, t = 8.0. Density (a) 1.8oo0, (b) 0.1700, (c) 0.0700, (d) 0.0300, (e) 0.0100, (f) 0.@080.
Gmuitational Hydrudms
FIG.3
Particle trajectoxies, t = 4.0.
FIG.4
Particle trajectories, t = 8.0.
located at the point {x(&J = 18, y( to) = 14.) with initial velocity { x’(ts) = 0.1935, y’(t,) = O.}. The fluid filling Q has initial constant density p0 = 0.1, and it is supposed to be initiahy at rest. Because of the initial conditions of the flmd, only the first step of the computation procedure is obtained by the numerical integration of the system (3) by using (5). Figures 1 and 2 show, at the times t = 4.0 and t = 8.0 respectively, the contour lines of the density p of the fluid for the following values: p = {0.170,0.150,0.090,0.050}
(t = 4.0),
p = { 1.800,0.170,0.070,0.030,0.010,0.008}
(t =S.O).
Figures 3 and 4 respectively show, at the same times t = 4.0 and t = 8.0, the trajectory of the particle S starting from the initial point and the fluid velocity field in which S is moving. We notice that in Figures 1 and 2 the density is increasing from the boundary to the center of 9, as can be easily deduced from the equations; moreover, in Figures 3 and 4 we can observe that the behavior of the velocity field and the path of S are consistent with the proposed model. Finally we would like to underline that the described method gives a promising tool for the study of the interactions between two classical gravitational fields. The extension of this method to a model taking into account also the energy equation and the interaction of several particles seems to be quite straightforward.
We wish to erpress our gmtitude to Giuuanni Pause of the Dipartimento di Mat&mat&a &l Politemico di Milano jin- his encoumgement and a useful
U. BULGABEILI,
264 suggestion jbr
hf. M. CEFUMEXE, AND A. ZABETTI
the development of this mod+d, and to Marcella Prosped of
Z.A.C. Rome, who supjdid
the jyaphic mutine.~.
REFERENCES 1
U. Bulgamlh, V. Casulli, and D. Greenspan, hssure Metti for the Num.ericuZ Solution of Free !&j&cc Fhdd Flows, Pineridge Press, Swansea, U.K., 1984.
2
U. BulgamIi, M. M. Cerimele, and A. Zaretti, A Numerical Method for a Mathematical Model Describing Interaction Between Two Gravitational Fields, Quad. IAC No. 16,1985. R. Courant and K. 0. Friedrh, SupersonicFh and Shock Waves, Interscience, New York, 1948. D. Greenspan, Disc&e Numerical Methodi in PhtrJics and Engineering, Academic, 1974. F. B. Hildebrandt, Finite Diffbmnce Equutiom and Simulutions,Prenti~Hall, 1988. 0. D. Kelkqp, Foundation of Potential Theog, Springer, Berlin, 1929. C. C. Lin, Stellar dynamical theory of normal spirals,in Lecturea in Appl. Math., Vol. 9 (G. Ehlers, Ed.), Amer. Math.!%c.,Providence, 1987. D. Potter, ChmputaiionulPh@cs, Wiley, New York, 1973. G. Prouse, F. Bohii, and A. Zaretti, On a mathematical model describing interaction between two classicalgravitatiod fields, Rend. Accad. Nuz. Sci. dei XL Mem. Mat. 11, to appear. J. E. Welch, F. H. H&w, J. P. Shannon,and B. J. Daly, The MAC Method: a computing technique for solving viscous, incompressible, transient fluid-flow problems involving frea surfaces, LA 3425, Los Alama Scientific Labs., Los Alamas, N.M., 1988. D. Trigiante, Asymptotic stability and discretization on an iufinite interval, Computing l&117-119 (1977).
10
11