Computer-readable recording medium, simulation method, and simulation device
A simulation program, simulation method, or simulation device simulates the temporal shape change in a fluid. The simulation program, simulation method, and simulation device are characterized in that with respect to a fluid model representing the fluid as a collection of particles, a surface tension of an interface against another phase which is different from the phase of the fluid is calculated using a function to calculate a surface energy of the interface, and the shape change in the fluid is calculated on the basis of the calculated surface tension.
1. A non-transitory computer-readable recording medium having stored therein a simulation program for causing a computer to execute a process for suppressing behaviors including generating an unnatural velocity distribution when simulating a temporal shape change in a fluid, the process comprising:
obtaining a position vector and a velocity vector of a particle forming the fluid from a storage unit;
calculating, with respect to a fluid model representing the fluid as a collection of particles, a force applied to each particle and a surface tension applied to each particle using the obtained position vector and velocity vector;
adding the calculated force applied to each particle and surface tension applied to each particle to the velocity vector so as to update a velocity of each particle;
calculating a shape change in the fluid using the updated velocity of each particle; and
outputting the calculated shape change to an output unit, wherein
[
Formula
12
]
v
i
*
,
n
+
1
=
v
_
i
n
-
2
ⅆ
t
[
∑
j
m
j
(
p
ij
n
+
1
/
2
ρ
i
n
ρ
j
n
)
(
L
(
r
_
i
*
)
+
L
(
r
_
i
*
)
2
)
∂
W
(
r
_
i
*
-
r
j
*
,
h
)
∂
r
i
*
)
-
1
2
∑
j
m
j
(
4
μξ
ρ
i
ρ
j
(
v
_
i
n
-
v
_
j
n
)
·
(
r
_
i
*
-
r
_
j
*
)
r
_
i
*
-
r
_
j
*
2
+
η
2
)
v
ij
i
∂
W
(
r
_
i
*
-
r
_
j
*
,
h
)
∂
r
_
i
*
-
g
y
^
+
f
ij
m
i
+
f
i
,
s
1
m
i
,
is used when updating the velocity of each particle, [Formula 12] including a second term on the right side expressing a pressure gradient term, a third term on the right side expressing a viscous stress term, a fourth term on the right side expressing a gravity term, and a fifth term on the right side which is a force for preventing the distance between particles from becoming too close, to calculate the force applied to each particle, as well as a sixth term on the right side which is a surface tension, and followed by the following calculation:
r
i
*
,
n
+
1
=
r
_
i
*
+
ⅆ
t
2
v
i
*
,
n
+
1
[
Formula
13
]
r
i
n
+
1
=
{
r
i
*
,
n
+
1
y
i
*
,
n
+
1
>
0
r
i
*
,
n
+
1
-
y
i
*
,
n
+
1
v
i
,
y
*
,
n
+
1
v
i
*
,
n
+
1
y
i
*
,
n
+
1
≤
0
,
v
i
n
+
1
=
{
v
i
*
,
n
+
1
y
i
*
,
n
+
1
>
0
(
1
-
y
i
*
,
n
+
1
dtv
i
,
y
*
,
n
+
1
)
v
i
*
,
n
+
1
y
i
*
,
n
+
1
≤
0
[
Formula
14
]
ρ
i
n
+
1
=
ρ
i
n
+
2
dt
∑
j
m
j
ρ
i
n
ρ
j
n
(
v
i
s
,
n
+
1
+
v
_
i
s
,
n
2
-
v
ij
s
,
n
+
1
/
2
)
∂
∂
(
r
_
i
*
-
r
_
j
*
)
W
(
r
_
i
*
-
r
_
j
*
,
h
)
,
[
Formula
15
]
wherein
a subscript expresses the number of the particle, that is,
r i , v i , ρ i , p i represent the position vector, velocity vector, density, and pressure of the i-th particle respectively,
m j is the mass of the j-th particle, and
ξ, η are constants introduced to calculate a viscous term, wherein
f
ij
=
8
r
min
2
(
1
-
(
r
_
i
*
-
r
_
j
*
r
min
)
2
)
r
_
i
*
-
r
_
j
*
r
_
i
*
-
r
_
j
*
(r min is a constant), wherein
[
Formula
6
]
f
k
,
s
1
=
-
∂
∂
r
k
E
s
,
p
1
(
r
1
,
r
2
,
…
,
r
n
)
≈
Q
∑
j
[
(
Γ
g
-
Γ
s
)
(
m
k
ρ
k
(
m
j
ρ
j
)
1
-
1
d
f
k
χ
j
+
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
χ
k
)
+
Γ
s
(
m
k
ρ
k
(
m
j
ρ
j
)
1
-
1
d
f
k
+
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
)
]
∂
W
(
r
k
-
r
j
,
h
)
∂
r
k
-
Q
(
Γ
g
-
Γ
s
)
∑
j
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
∂
χ
k
∂
r
k
W
(
r
k
-
r
j
)
is applied, and
χ i =χ i ( y i *, ε)
if a plain surface of y=0 is a wall surface of a solid, wherein
g is gravity acceleration,
y is a unit vector in the y direction,
h expresses an influence radius which is a radius of a spherical area in which a value of a kernel function is non-zero, and
W is an overlap of kernel functions used in the SPH method, wherein
L( r i *)
is a renormalization matrix of
L
(
r
i
)
=
[
∑
j
m
j
ρ
j
(
r
j
-
r
i
)
⊗
∂
∂
r
i
W
(
r
i
-
r
j
,
h
)
]
-
1
,
and if
n
_
ij
*
=
(
L
(
r
_
i
*
)
+
L
(
r
_
i
*
)
2
)
3
r
_
i
*
-
r
_
j
*
r
_
i
*
-
r
_
j
*
v i s,n = v i n · n ij *
v j s,n = v j n · n ij *,
wherein
p ij n+1/2 ,v ij s,n+1/2
are intermediate values of space-time determined by solving the one-dimensional Riemann problem between the particles i and j, the values being determined by a process including:
first setting the following characteristic quantities with respect to the particles i and j:
q
i
n
,
+
=
log
(
ρ
i
n
)
+
v
i
s
,
n
c
[
Formula
16
]
q
i
n
,
-
=
log
(
ρ
i
n
)
-
v
i
s
,
n
c
[
Formula
17
]
q
j
n
,
+
=
log
(
ρ
j
n
)
+
v
j
s
,
n
c
[
Formula
18
]
q
j
n
,
-
=
log
(
ρ
j
n
)
-
v
j
s
,
n
c
;
[
Formula
19
]
and
then calculating each gradient as follows:
∇
log
(
ρ
)
|
i
=
∑
k
m
k
ρ
i
(
log
(
ρ
k
n
)
-
log
(
ρ
i
n
)
)
L
(
r
_
i
*
)
∂
W
(
r
_
i
*
-
r
_
k
*
,
h
)
∂
r
_
i
*
[
Formula
20
]
∇
v
|
i
=
∑
k
m
k
ρ
k
L
(
r
_
i
*
)
∂
W
(
r
_
i
*
-
r
_
k
*
,
h
)
∂
r
_
i
*
⊗
(
v
_
k
n
-
v
_
i
n
)
[
Formula
21
]
∇
q
i
n
,
+
=
∇
log
(
ρ
)
i
+
∇
v
|
i
n
_
ij
*
c
[
Formula
22
]
∇
q
j
n
,
+
=
∇
log
(
ρ
)
j
+
∇
v
|
j
n
_
ij
*
c
[
Formula
23
]
∇
q
i
n
,
-
=
∇
log
(
ρ
)
i
+
∇
v
|
i
n
_
ij
*
c
[
Formula
24
]
∇
q
j
n
,
-
=
∇
log
(
ρ
)
j
-
∇
v
|
j
n
_
ij
*
c
,
[
Formula
25
]
and wherein
r ij *= r i *− r j *,
r ij *=| r ij *|
and using the above quantities,
p ij n+1/2 ,v ij n+1/2
are determined as follows:
q
ij
n
+
1
/
2
,
+
=
q
j
n
,
+
+
(
r
_
ij
*
2
-
cdt
2
)
(
r
_
ij
*
·
∇
q
|
j
n
,
+
)
[
Formula
26
]
q
ij
n
+
1
/
2
,
-
=
q
i
n
,
+
-
(
r
_
ij
*
2
+
cdt
2
)
(
r
_
ij
*
·
∇
q
|
j
n
,
-
)
[
Formula
27
]
ρ
ij
n
+
1
/
2
=
exp
(
q
ij
n
+
1
/
2
,
+
+
q
ij
n
+
1
/
2
,
-
2
)
[
Formula
28
]
v
ij
n
+
1
/
2
=
c
(
q
ij
n
+
1
/
2
,
+
-
q
ij
n
+
1
/
2
,
-
2
)
[
Formula
29
]
p
ij
n
+
1
/
2
=
c
2
(
ρ
ij
n
+
1
/
2
+
ρ
0
)
.
[
Formula
30
]
2. The non-transitory computer-readable recording medium according to claim 1 , wherein the function is for calculating the sum of the result of the calculation of the surface energy of the interface with a phase of the air, and the result of the calculation of the surface energy of the interface with a phase of other than the phase of the air.
3. The non-transitory computer-readable recording medium according to claim 1 , wherein in the calculation of a shape change, the particle forming the fluid is stopped when the fluid collides against a wall surface.
4. A simulation method for suppressing behaviors including generating an unnatural velocity distribution when simulating a temporal shape change in a fluid by using a processor, the simulation method comprising:
obtaining a position vector and a velocity vector of a particle forming the fluid from a storage unit by using the processor;
calculating, with respect to a fluid model representing the fluid as a collection of particles, a force applied to each particle and a surface tension applied to each particle using the obtained position vector and velocity vector by using the processor;
adding the calculated force applied to each particle and surface tension applied to each particle to the velocity vector so as to update a velocity of each particle by using the processor;
calculating a shape change in the fluid using the updated velocity of each particle by using the processor; and
outputting the calculated shape change to an output unit by using the processor, wherein
[
Formula
12
]
v
i
*
,
n
+
1
=
v
_
i
n
-
2
dt
[
∑
j
m
j
(
p
ij
n
+
1
/
2
ρ
i
n
ρ
j
n
)
(
L
(
r
_
i
*
)
+
L
(
r
_
i
*
)
2
)
∂
W
(
r
_
i
*
-
r
j
*
,
h
)
∂
r
i
*
]
-
1
2
∑
j
m
j
(
4
μξ
ρ
i
ρ
j
(
v
_
i
n
-
v
_
j
n
)
·
(
r
_
i
*
-
r
_
j
*
)
r
_
i
*
-
r
_
j
*
2
+
η
2
)
v
ij
i
∂
W
(
r
_
i
*
-
r
_
j
*
,
h
)
∂
r
_
i
*
-
g
y
^
+
f
ij
m
i
+
f
i
,
s
1
m
i
,
is used when updating the velocity of each particle, [Formula 12] including a second term on the right side expressing a pressure gradient term, a third term on the right side expressing a viscous stress term, a fourth term on the right side expressing a gravity term, and a fifth term on the right side which is a force for preventing the distance between particles from becoming too close, to calculate the force applied to each particle, as well as a sixth term on the right side which is a surface tension, and followed by the following calculation:
r
i
*
,
n
+
1
=
r
_
i
*
+
dt
2
v
i
*
,
n
+
1
[
Formula
13
]
r
i
n
+
1
=
{
r
i
*
,
n
+
1
y
i
*
,
n
+
1
>
0
r
i
*
,
n
+
1
-
y
i
*
,
n
+
1
v
i
,
y
*
,
n
+
1
v
i
*
,
n
+
1
y
i
*
,
n
+
1
≤
0
,
v
i
n
+
1
=
{
v
i
*
,
n
+
1
y
i
*
,
n
+
1
>
0
(
1
-
y
i
*
,
n
+
1
dtv
i
,
y
*
,
n
+
1
)
v
i
*
,
n
+
1
y
i
*
,
n
+
1
≤
0
[
Formula
14
]
ρ
i
n
+
1
=
ρ
i
n
+
2
dt
∑
j
m
j
ρ
i
n
ρ
j
n
(
v
i
s
,
n
+
1
+
v
_
i
s
,
n
2
-
v
ij
s
,
n
+
1
/
2
)
∂
∂
(
r
_
i
*
-
r
_
j
*
)
W
(
r
_
i
*
-
r
_
j
*
,
h
)
,
[
Formula
15
]
wherein
a subscript expresses the number of the particle, that is,
r i , v i , ρ i , p i represent the position vector, velocity vector, density, and pressure of the i-th particle respectively,
m j is the mass of the j-th particle, and
ξ, η are constants introduced to calculate a viscous term, wherein
f
ij
=
8
r
min
2
(
1
-
(
r
_
i
*
-
r
_
j
*
r
min
)
2
)
3
r
_
i
*
-
r
_
j
*
r
_
i
*
-
r
_
j
*
(r min is a constant), wherein
f
k
,
s
1
=
-
∂
∂
r
k
E
s
,
p
1
(
r
1
,
r
2
,
…
,
r
n
)
≈
Q
∑
j
[
(
Γ
g
-
Γ
s
)
(
m
k
ρ
k
(
m
j
ρ
j
)
1
-
1
d
f
k
χ
j
+
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
χ
k
)
+
Γ
s
(
m
k
ρ
k
(
m
j
ρ
k
)
1
-
1
d
f
k
+
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
)
]
∂
W
(
r
k
-
r
j
,
h
)
∂
r
k
-
Q
(
Γ
g
-
Γ
s
)
∑
j
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
∂
χ
k
∂
r
k
W
(
r
k
-
r
j
)
[
Formula
6
]
is applied, and
χ i =χ i ( y i *:ε)
if a plain surface of v=0 is a wall surface of a solid, wherein
g is gravity acceleration,
y is a unit vector in the y direction,
h expresses an influence radius which is a radius of a spherical area in which a value of a kernel function is non-zero, and
W is an overlap of kernel functions used in the SPH method, wherein
L( r i *)
is a renormalization matrix of
L
(
r
i
)
=
[
∑
j
m
j
ρ
j
(
r
j
-
r
i
)
⊗
∂
∂
r
i
W
(
r
i
-
r
j
,
h
)
]
-
1
,
and if
n
_
ij
*
=
(
L
(
r
_
i
*
)
+
L
(
r
_
i
*
)
2
)
r
_
i
*
-
r
_
j
*
r
_
i
*
-
r
_
j
*
v i s,n = v i n · n ij *
v j s,n = v j n · n ij *
, wherein
p ij n+1/2 ,v ij s,n+1/2
are intermediate values of space-time determined by solving the one-dimensional Riemann problem between the particles i and i, the values being determined by a process including:
first setting the following characteristic quantities with respect to the particles i and j:
q
i
n
,
+
=
log
(
ρ
i
n
)
+
v
i
s
,
n
c
[
Formula
16
]
q
i
n
,
-
=
log
(
ρ
i
n
)
-
v
i
s
,
n
c
[
Formula
17
]
q
j
n
,
+
=
log
(
ρ
j
n
)
+
v
j
s
,
n
c
[
Formula
18
]
q
j
n
,
-
=
log
(
ρ
j
n
)
-
v
j
s
,
n
c
;
[
Formula
19
]
and
then calculating each gradient as follows:
∇
log
(
ρ
)
|
i
=
∑
k
m
k
ρ
i
(
log
(
ρ
k
n
)
-
log
(
ρ
i
n
)
)
L
(
r
_
i
*
)
∂
W
(
r
_
i
*
-
r
_
k
*
,
h
)
∂
r
_
i
*
[
Formula
20
]
∇
v
|
i
=
∑
k
m
k
ρ
k
L
(
r
_
i
*
)
∂
W
(
r
_
i
*
-
r
_
k
*
,
h
)
∂
r
_
i
*
⊗
(
v
_
k
n
-
v
_
i
n
)
[
Formula
21
]
∇
q
i
n
,
+
=
∇
log
(
ρ
)
i
+
∇
v
|
i
n
_
ij
*
c
[
Formula
22
]
∇
q
j
n
,
+
=
∇
log
(
ρ
)
j
+
∇
v
|
j
n
_
ij
*
c
[
Formula
23
]
∇
q
i
n
,
-
=
∇
log
(
ρ
)
i
-
∇
v
|
i
n
_
ij
*
c
[
Formula
24
]
∇
q
j
n
,
-
=
∇
log
(
ρ
)
j
-
∇
v
|
j
n
_
ij
*
c
,
[
Formula
25
]
and wherein
r ij *= r i *− r j *,
r ij *=| r ij *|
and using the above quantities,
p ij n+1/2 ,v ij s,n+1/2
are determined as follows:
q
ij
n
+
1
/
2
,
+
=
q
j
n
,
+
+
(
r
_
ij
*
2
-
cdt
2
)
(
r
_
ij
*
·
∇
q
|
j
n
,
+
)
[
Formula
26
]
q
ij
n
+
1
/
2
,
-
=
q
i
n
,
+
-
(
r
_
ij
*
2
+
cdt
2
)
(
r
_
ij
*
·
∇
q
|
i
n
,
-
)
[
Formula
27
]
ρ
ij
n
+
1
/
2
=
exp
(
q
ij
n
+
1
/
2
+
q
ij
n
+
1
/
2
,
-
2
)
[
Formula
28
]
v
ij
n
+
1
/
2
=
c
(
q
ij
n
+
1
/
2
,
+
-
q
ij
n
+
1
/
2
,
-
2
)
[
Formula
29
]
p
ij
n
+
1
/
2
=
c
2
(
ρ
ij
n
+
1
/
2
+
ρ
0
)
.
[
Formula
30
]
5. A simulation device for suppressing behaviors including generating an unnatural velocity distribution when simulating a temporal shape change in a fluid, the simulation device comprising:
a processor which performs a process including:
obtaining a position vector and a velocity vector of a particle forming the fluid from a storage unit,
calculating, with respect to a fluid model representing the fluid as a collection of particles, a force applied to each particle and a surface tension applied to each particle using the obtained position vector and velocity vector;
adding the calculated force applied to each particle and surface tension applied to each particle to the velocity vector so as to update a velocity of each particle;
calculating a shape change in the fluid using the updated velocity of each particle; and
an output unit configured to output the calculated shape change, wherein
v
i
*
,
n
+
1
=
v
_
i
n
-
2
ⅆ
t
[
∑
j
m
j
(
p
ij
n
+
1
/
2
ρ
i
n
ρ
j
n
)
(
L
(
r
_
i
*
)
+
L
(
r
_
i
*
)
2
)
∂
W
(
r
_
i
*
-
r
_
j
*
,
h
)
∂
r
i
*
]
-
1
2
∑
j
m
j
(
4
μξ
ρ
i
ρ
j
(
v
_
i
n
-
v
_
j
n
)
·
(
r
_
i
*
-
r
_
j
*
)
r
_
i
*
-
r
_
j
*
2
+
η
2
)
v
ij
i
∂
W
(
r
_
i
*
-
r
_
j
*
,
h
)
∂
r
i
*
-
g
y
^
+
f
ij
m
i
+
f
i
,
s
1
m
i
,
[
Formula
12
]
is used when updating the velocity of each particle, [Formula 12] including a second term on the right side expressing a pressure gradient term, a third term on the right side expressing a viscous stress term, a fourth term on the right side expressing a gravity term, and a fifth term on the right side which is a force for preventing the distance between particles from becoming too close, to calculate the force applied to each particle, as well as a sixth term on the right side which is a surface tension, and followed by the following calculation:
r
i
*
,
n
+
1
=
r
_
i
*
+
ⅆ
t
2
v
i
*
,
n
+
1
[
Formula
13
]
r
i
n
+
1
=
{
r
i
*
,
n
+
1
y
i
*
,
n
+
1
>
0
r
i
*
,
n
+
1
-
y
i
*
,
n
+
1
v
i
,
y
*
,
n
+
1
y
i
*
,
n
+
1
≤
0
,
v
i
n
+
1
=
{
v
i
*
,
n
+
1
y
i
*
,
n
+
1
>
0
(
1
-
y
i
*
,
n
+
1
dtv
i
,
y
*
,
n
+
1
)
v
i
*
,
n
+
1
y
i
*
,
n
+
1
≤
0
[
Formula
14
]
ρ
i
n
+
1
=
ρ
i
n
+
2
dt
∑
j
m
j
ρ
i
n
ρ
j
n
(
v
i
s
,
n
+
1
+
v
_
i
s
,
n
2
-
v
ij
s
,
n
+
1
/
2
)
∂
∂
(
r
_
i
*
-
r
_
j
*
)
W
(
r
_
i
*
-
r
_
j
*
,
h
)
,
[
Formula
15
]
wherein
a subscript expresses the number of the particle, that is,
r i ,v i , ρ i , p i represent the position vector, velocity vector, density, and pressure of the i-th particle respectively,
m j is the mass of the j-th particle, and
ξ, η are constants introduced to calculate a viscous term, wherein
f
ij
=
8
r
min
2
(
1
-
(
r
_
i
*
-
r
_
j
*
r
min
)
2
)
3
r
_
i
*
-
r
_
j
*
r
_
i
*
-
r
_
j
*
(r min is a constant), wherein
f
k
,
s
1
=
-
∂
∂
r
k
E
s
,
p
1
(
r
1
,
r
2
,
…
,
r
n
)
≈
Q
∑
j
[
(
Γ
g
-
Γ
s
)
(
m
k
ρ
k
(
m
j
ρ
j
)
1
-
1
d
f
k
χ
j
+
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
χ
k
)
+
Γ
s
(
m
k
ρ
k
(
m
j
ρ
k
)
1
-
1
d
f
k
+
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
)
]
∂
W
(
r
k
-
r
j
,
h
)
∂
r
k
-
Q
(
Γ
g
-
Γ
s
)
∑
j
m
j
ρ
j
(
m
k
ρ
k
)
1
-
1
d
f
j
∂
χ
k
∂
r
k
W
(
r
k
-
r
j
)
[
Formula
6
]
is applied, and
χ i =χ i ( y i *,ε)
if a plain surface of y=0 is a wall surface of a solid, wherein
g is gravity acceleration,
y is a unit vector in the y direction,
h expresses an influence radius which is a radius of a spherical area in which a value of a kernel function is non-zero, and
W is an overlap of kernel functions used in the SPH method, wherein
L( r i *)
is a renormalization matrix of
L
(
r
i
)
=
[
∑
j
m
j
ρ
j
(
r
j
-
r
i
)
⊗
∂
∂
r
i
W
(
r
i
-
r
j
,
h
)
]
-
1
,
and if
n
_
ij
*
=
(
L
(
r
_
i
*
)
+
L
(
r
_
i
*
)
2
)
r
_
i
*
-
r
_
j
*
r
_
i
*
-
r
_
j
*
v i s,n = v i n · n ij *
v j s,n = v j n · n ij *
, wherein
p ij n+1/2 ,v ij s,n+1/2
are intermediate values of space-time determined by solving the one-dimensional Riemann problem between the particles i and j, the values being determined by a process including:
first setting the following characteristic quantities with respect to the particles i and j:
q
i
n
,
+
=
log
(
ρ
i
n
)
+
v
i
s
,
n
c
[
Formula
16
]
q
i
n
,
-
=
log
(
ρ
i
n
)
-
v
i
s
,
n
c
[
Formula
17
]
q
j
n
,
+
=
log
(
ρ
j
n
)
+
v
j
s
,
n
c
[
Formula
18
]
q
j
n
,
-
=
log
(
ρ
j
n
)
-
v
j
s
,
n
c
;
[
Formula
19
]
and
then calculating each gradient as follows:
∇
log
(
ρ
)
|
i
=
∑
k
m
k
ρ
i
(
log
(
ρ
k
n
)
-
log
(
ρ
i
n
)
)
L
(
r
_
i
*
)
∂
W
(
r
_
i
*
-
r
_
k
*
,
h
)
∂
r
_
i
*
[
Formula
20
]
∇
v
|
i
=
∑
k
m
k
ρ
k
L
(
r
_
i
*
)
∂
W
(
r
_
i
*
-
r
_
k
*
,
h
)
∂
r
_
i
*
⊗
(
v
_
k
n
-
v
_
i
n
)
[
Formula
21
]
∇
q
i
n
,
+
=
∇
log
(
ρ
)
i
+
∇
v
|
i
n
_
ij
*
c
[
Formula
22
]
∇
q
j
n
,
+
=
∇
log
(
ρ
)
j
+
∇
v
|
j
n
_
ij
*
c
[
Formula
23
]
∇
q
i
n
,
-
=
∇
log
(
ρ
)
i
-
∇
v
|
i
n
_
ij
*
c
[
Formula
24
]
∇
q
j
n
,
-
=
∇
log
(
ρ
)
j
-
∇
v
|
j
n
_
ij
*
c
,
[
Formula
25
]
and wherein
r ij *= r i *− r j *,
r ij *=| r ij *|
and using the above quantities,
p ij n+1/2 ,v ij s,n+1/2
are determined as follows:
q
ij
n
+
1
/
2
,
+
=
q
j
n
,
+
+
(
r
_
ij
*
2
-
cdt
2
)
(
r
_
ij
*
·
∇
q
|
j
n
,
+
)
[
Formula
26
]
q
ij
n
+
1
/
2
,
-
=
q
i
n
,
+
-
(
r
_
ij
*
2
+
cdt
2
)
(
r
_
ij
*
·
∇
q
|
i
n
,
-
)
[
Formula
27
]
ρ
ij
n
+
1
/
2
=
exp
(
q
ij
n
+
1
/
2
+
q
ij
n
+
1
/
2
,
-
2
)
[
Formula
28
]
v
ij
n
+
1
/
2
=
c
(
q
ij
n
+
1
/
2
,
+
-
q
ij
n
+
1
/
2
,
-
2
)
[
Formula
29
]
p
ij
n
+
1
/
2
=
c
2
(
ρ
ij
n
+
1
/
2
+
ρ
0
)
.
[
Formula
30
]