Wave Propagation Methods for
Compressible Multicomponent Flow with
Moving Interfaces and Boundaries
Keh-Ming Shyue
Department of Mathematics National Taiwan University
Overview
Illustrative examples Mathematical model
Fluid-mixture type equations of motion for homogeneous two-phase flow
General pressure law for real materials Numerical techniques
Finite volume method based on wave propagation Surface tracking for moving boundaries
Volume tracking for moving interfaces Numerical results
Future work
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900 1000
air t = 0
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900 1000
air
t = 0.2s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 0.4s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 0.6s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 0.8s
Falling Liquid Drop Problem
100 200 300 400 500 600 700 800 900
air
t = 1s
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water t = 0
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water t = 12ms
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water t = 24ms
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water t = 36ms
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water t = 48ms
Flying Projectile & Ocean Surface
−10 0 10 20 30 40 50
−20
−15
−10
−5 0 5 10 15 20
air
water t = 60ms
Two Phase Flow Problem
Ignore physical effects such as gravity, viscosity, surface tension, mass diffusion, and so on
Each fluid component satisties
Eulerian form conservation laws
ρt + ∇ · (ρ~u) = 0 (ρ~u)t + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (ρE)t + ∇ · (ρE~u + p~u) = 0 General pressure law p(ρ, e)
ρ: density, ~u: vector of particle velocity, p: pressure E: specific total energy, e: specific internal energy
Two-Phase Flow Model
Model derivation based on averaging theory of Drew (Theory of Multicomponent Fluids, D.A. Drew & S. L.
Passman, Springer, 1999)
Namely, introduce indicator function χk as χk(M, t) =
(1 if M belongs to phase k
0 otherwise
Denote < ψ > as volume averaged for flow variable ψ, hψi = 1
V Z
V
ψ dV Gauss & Leibnitz rules
hχk∇ψi = h∇(χkψ)i−hψ∇χki & hχkψti = h(χkψ)ti−hψ(χk)ti
Two-Phase Flow Model (cont.)
Take product of each conservation law with χk & perform averaging process. In case of mass conservation equation, for example, we have
hχkρkit + ∇· < χkρk~uk >= hρk(χk)t + ρk~uk · ∇χki Since χk is governed by
(χk)t + ~u0 · ∇χk = 0 (~u0: interface velocity), this leads to mass averaged equation for phase k
hχkρkit + ∇· < χkρk~uk >= hρk (~uk − ~u0) · ∇χki Analogously, we may derive averaged equation for momentum, energy, & entropy (not shown here)
Two-Phase Flow Model (cont.)
In summary, averaged model system, we have, are
hχkρkit + ∇· < χkρk~uk >= hρk (~uk − ~u0) · ∇χki hχkρk~ukit + ∇· < χkρk~uk ⊗ ~uk > +∇ hχkpki = hpk∇χki +
hρk~uk (~uk − ~u0) · ∇χki hχkρkEkit + ∇· < χkρkEk~uk + χkpk~uk >= hpk~uk · ∇χki +
hρkE (~uk − ~u0) · ∇χki hχkit + h~uk · ∇χki = h(~uk − ~u0) · ∇χki
Note: existence of various interfacial source terms
Mathematical as well as numerical modelling of these terms are important (but difficult) for general multiphase flow
problems
Homogeneous 2-Phase Flow Model
Assume homogeneous (1-pressure & 1-velocity) flow;
i.e., across interfaces: pι = p & ~uι = ~u, ι = 0, 1, 2
Introduce volume fraction for phase k as αk = Vk/V Now, by dropping all interfacial terms, we may obtain a simplified model as
(αkρk)t + ∇ · (αkρk~u) = 0
(αkρk~u)t + ∇ · (αkρk~u ⊗ ~u) + ∇ (αkp) = p∇αk (αkρkEk)t + ∇ · (αkρkEk~u + αkp~u) = p~u · ∇αk
(α1)t + ~u · ∇α1 = 0
for k = 1, 2, & α1 + α2 = 1. Note this gives 2(2 + Nd) + 1
equations in total for a Nd-dimension 2-phase flow problem
Homogeneous Flow Model (cont.)
Note that, in practice, rather than using equations αkρk~u &
αkρkEk for each phase, we may write down a system of the form
(αkρk)t + ∇ · (αkρk~u) = 0 (ρ~u)t + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (ρE)t + ∇ · (ρE~u + p~u) = 0 (α1)t + ~u · ∇α1 = 0
ρ~u = P2
k=1 αkρk~u: total momentum ρE = P2
k=1 αkρkEk: total energy
This gives 4 + Nd equations in total, Nd + 1 less than previous model system
Homogeneous Flow Model (cont.)
Note that it is easy to include, for instance, gravity &
capillary effects to the model
(αkρk)t + ∇ · (αkρk~u) = 0 (k = 1, 2) (ρ~u)t + ∇ · (ρ~u ⊗ ~u) + ∇p = φ~
(ρE)t + ∇ · (ρE~u + p~u) = φ · ~u~ (α1)t + ~u · ∇α1 = 0
1. Gravity case: φ = ~g~
2. Capillary case: φ = σκ∇α~
~g: gravitational constant, σ: surface tension coef.
κ: curvature at interface
Homogeneous Flow Model (cont.)
Mixture equation of state: p = p(α2, α1ρ1, α2ρ2, ρe) Isobaric closure: p1 = p2 = p
For a class of EOS, explicit formula for p is available (examples are given next)
For some complex EOS, from (α2, ρ1, ρ2, ρe) in model equations we recover p by solving
p1(ρ1, ρ1e1) = p2(ρ2, ρ2e2) &
2
X
k=1
αkρkek = ρe This homogeneous two-phase model was called a five-equation model by Allaire, Clerc, & Kokh (JCP
2002) or a volume-fraction model by Shyue (JCP 1998)
Homogeneous Flow Model (cont.)
Polytropic ideal gas: pk = (γk − 1)ρkek
ρe =
2
X
k=1
αkρkek =
2
X
k=1
αk p
γk − 1 =⇒
p = ρe
2 X
k=1
αk γk − 1
Homogeneous Flow Model (cont.)
Polytropic ideal gas: pk = (γk − 1)ρkek
ρe =
2
X
k=1
αkρkek =
2
X
k=1
αk p
γk − 1 =⇒
p = ρe
2 X
k=1
αk γk − 1
Van der Waals gas: pk = (1−γkb−1kρk )(ρkek + akρ2k) − akρ2k
ρe =
2
X
k=1
αkρkek =
2
X
k=1
αk 1 − bkρk γk − 1
(p + akρ2k) − akρ2k
=⇒
p =
"
ρe −
2
X
k=1
αk 1 − bkρk
γk − 1 − 1
akρ2k
# 2 X
k=1
αk 1 − bkρk
γk − 1
Homogeneous Flow Model (cont.)
Two-molecular vibrating gas: pk = ρkRkT (ek), T satisfies
e = RT
γ − 1 + RTvib
exp Tvib/T − 1
As before, we now have
ρe =
2
X
k=1
αkρkek =
2
X
k=1
αk
ρkRkTk γk − 1
+ ρkRkTvib,k exp
Tvib,k/Tk
− 1
=
2
X
k=1
αk
p
γk − 1
+ pvib,k exp
pvib,k/p
− 1
(Nonlinear eq.)
Homogeneous Flow Model (cont.)
It is easy to show entropies, Sk, k = 1, 2, satisfy
∂p1
∂S1
ρ1
DS1
Dt − ∂p2
∂S2
ρ2
DS2
Dt = ρ1c21 − ρ2c22 ∇ · ~u
Homogeneous Flow Model (cont.)
It is easy to show entropies, Sk, k = 1, 2, satisfy
∂p1
∂S1
ρ1
DS1
Dt − ∂p2
∂S2
ρ2
DS2
Dt = ρ1c21 − ρ2c22 ∇ · ~u
Murrone & Guillard (JCP 2005) propsed a reduced two-phase flow model in which
(α1)t + ~u · ∇α1 = α1α2 ρ2c22 − ρ1c21
P2
k=1 αkρkc2k
!
and now entropy of each phase satisfy
DSk
Dt = ∂Sk
∂t + ~u · ∇Sk = 0, for k = 1, 2
Some Remarks
1. Model system is hyperbolic under suitable thermodynamic stability condition
Some Remarks
1. Model system is hyperbolic under suitable thermodynamic stability condition
2. When α2 = 0 (or = 1), ρ2 (or ρ1) can not be recovered from α2 & α2ρ2 (or α1 & α1ρ1), and so take αk ∈ [ε, 1 − ε]
Some Remarks
1. Model system is hyperbolic under suitable thermodynamic stability condition
2. When α2 = 0 (or = 1), ρ2 (or ρ1) can not be recovered from α2 & α2ρ2 (or α1 & α1ρ1), and so take αk ∈ [ε, 1 − ε]
3. In the model, it is not at all clear on how to compute nonlinear term ρι, ι > 1 from αk & αkρk
Some Remarks
1. Model system is hyperbolic under suitable thermodynamic stability condition
2. When α2 = 0 (or = 1), ρ2 (or ρ1) can not be recovered from α2 & α2ρ2 (or α1 & α1ρ1), and so take αk ∈ [ε, 1 − ε]
3. In the model, it is not at all clear on how to compute nonlinear term ρι, ι > 1 from αk & αkρk
4. In fact, there are other set of model systems proposed in the literature that are more robust for homogeneous flow & in other more complicated context (examples)
Some Remarks
1. Model system is hyperbolic under suitable thermodynamic stability condition
2. When α2 = 0 (or = 1), ρ2 (or ρ1) can not be recovered from α2 & α2ρ2 (or α1 & α1ρ1), and so take αk ∈ [ε, 1 − ε]
3. In the model, it is not at all clear on how to compute nonlinear term ρι, ι > 1 from αk & αkρk
4. In fact, there are other set of model systems proposed in the literature that are more robust for homogeneous flow & in other more complicated context (examples) 5. In cases when individual pressure law differs in form
(see below), new mixture pressure law should be
devised first & construct model equations based on that
Barotropic & Non-Barotropic Flow
Fluid component 1: Tait EOS
p(ρ) = (p0 + B) ρ ρ0
γ
− B
Fluid component 2: Noble-Abel EOS
p(ρ, e) = γ − 1 1 − bρ
ρe
Mixture pressure law (Shyue, Shock Waves 2006)
p =
(p0 + B) ρ ρ0
γ
− B if α = 1
γ − 1 1 − bρ
(ρe − B) − B if α 6= 1
Barotropic Two-Phase Flow
Fluid component ι: Tait EOS
p(ρ) = (p0ι + Bι)
ρ ρ0ι
γι
− Bι, ι = 1, 2
Mixture pressure law (Shyue, JCP 2004)
p =
(p0ι + Bι)
ρ ρ0ι
γι
− Bι if α = αι (0 or 1) (γ − 1) ρ
e + B ρ0
− γB if α ∈ (0, 1)
Wave Propagation Method
Finite volume formulation of wave propagation method, QnS gives approximate value of cell average of solution q over cell S at time tn
QnS ≈ 1 M(S)
Z
S
q(X, tn) dV
M(S): measure (area in 2D or volume in 3D) of cell S
C E
D F
G H
Wave Propagation Method (cont.)
First order version: Piecewise constant wave update Godunov-type method: Solve Riemann problem at each cell interface in normal direction & use resulting waves to update cell averages
Qn+1S := Qn+1S −M (Wp ∩ S)
M(S) Rp, Rp being jump from RP
↓
↓
Wp Wp
Wave Propagation Method (cont.)
First order version: Transverse-wave included
Use transverse portion of equation, solve Riemann problem in transverse direction, & use resulting
waves to update cell averages as usual
Stability of method is typically improved, while conservation of method is maintained
↓ ↓
Wpq
Wpq
Wave Propagation Method (cont.)
High resolution version: Piecewise linear wave update wave before propagation after propagation
a) b)
c) d)
αpr
p/2
αpr
p/2
λp∆ t
λp∆ t
Volume Tracking Algorithm
1. Volume moving procedure (a) Volume fraction update
Take a time step on current grid to update cell averages of volume fractions at next time step (b) Interface reconstruction
Find new interface location based on volume fractions obtained in (a) using an interface reconstruction scheme. Some cells will be subdivided & values in each subcell must be initialized.
2. Physical solution update
Take same time interval as in (a), but use a method to update cell averages of multicomponent model on new grid created in (b)
Interface Reconstruction Scheme
Given volume fractions on current grid, piecewise linear interface reconstruction (PLIC) method does:
1. Compute interface normal
Gradient method of Parker & Youngs Least squares method of Puckett
2. Determine interface location by iterative bisection
0 0
0 0
0 0
0 0
0 0
0
0 0
0 0
0 0
0 0
0
1 0.29 0.68 0.09
0.51 0.51
0.09 0.68 0.29
↓
↓ interface
interface
Data set Parker & Youngs Puckett
Volume Moving Procedure
(a) Volume fractions given in previous slide are updated with uniform (u, v) = (1, 1) over ∆t = 0.06
(b) New interface location is reconstructed
0
0 0 0 0
0 0 0 0
0 0
0 0
0 0
0 0
0
1 0.01
0.38 0.11
0.25 0.72
0.06 0.85
0.74 5(−3) 1(−3)
↑
↓
old interface
new interface
(a) (b)
Surface Moving Procedure
Solve Riemann problem at tracked interfaces & use
resulting wave speed of the tracked wave family over ∆t to find new location of interface at the next time step
o
o
o o
o
o
↑
↑ old front
old front
new front
Boundary Conditions
For tracked segments representing rigid (solid wall)
boundary (stationary or moving), reflection principle is used to assign states for fictitious subcells in each time step:
zC := zE (z = ρ, p, α)
~uC := ~uE − 2(~uE · ~n)~n + 2(~u0 · ~n)
~u0: moving boundary velocity
C
E F
G H
~n
Interface Conditions
For tracked segments representing material interfaces, pressure equilibrium as well as velocity continuity
conditions across interfaces are fulfilled by 1. Devise of the wave-propagation method
2. Choice of Riemann solver used in the method
Stability Issues
Choose time step ∆t based on uniform grid mesh size
∆x, ∆y as
∆t maxp,q (λp, µq)
min(∆x, ∆y) ≤ 1,
λp, µq: speed of p-wave, q-wave from Riemann
problem solution in normal-, transverse-directions Use large time step method of LeVeque (i.e., wave
interactions are assumed to behave in linear manner) to maintain stability of method even in the presence of
small Cartesian cut cells
Apply smoothing operator (such as, h-box approach of Berger et al. ) locally for cell averages in irregular cells
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction Problem
Shock-Bubble Interaction (cont.)
Approximate locations of interfaces
time=55µs
air R22
time=115µs time=135µs
time=187µs time=247µs time=200µs
time=342µs time=417µs time=1020µs
Shock-Bubble Interaction (cont.)
Quantitative assessment of prominent flow velocities:
Velocity (m/s) Vs VR VT Vui Vuf Vdi Vdf Haas & Sturtevant 415 240 540 73 90 78 78 Quirk & Karni 420 254 560 74 90 116 82 Our result (tracking) 411 243 538 64 87 82 60 Our result (capturing) 411 244 534 65 86 98 76
Vs (VR, VT) Incident (refracted, transmitted) shock speed t ∈ [0, 250]µs (t ∈ [0, 202]µs, t ∈ [202, 250]µs ) Vui (Vuf) Initial (final) upstream bubble wall speed t ∈ [0, 400]µs (t ∈ [400, 1000]µs)
Vdi (Vdf) Initial (final) downstream bubble wall speed t ∈ [200, 400]µs (t ∈ [400, 1000]µs)
Aluminum-Plate Impact Problem
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 0µs
Aluminum-Plate Impact Problem
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 0.5µs
Aluminum-Plate Impact Problem
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 1µs
Aluminum-Plate Impact Problem
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 2µs
Aluminum-Plate Impact Problem
0 1 2 3
−2
−1 0 1 2
Density
Al target Al flyer
vacuum
0 1 2 3
−2
−1 0 1 2
Pressure
0 1 2 3
−2
−1 0 1 2
Volume fraction
t = 4µs
Cylinder lift-off Problem
Moving speed of cylinder is governed by Newton’s law Pressure contours are shown with a 1000 × 200 grid
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.05 0.1 0.15 0.2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.05 0.1 0.15 0.2
0 0.05 0.1 0.15 0.2
t = 0
t = 0.1641s
t = 0.30085s
Cylinder lift-off Problem
A convergence study of center of cylinder & relative mass loss for at final stopping time t = 0.30085s
Mesh size Center of cylinder Relative mass loss 250 × 50 (0.618181, 0.134456) −0.257528
500 × 100 (0.620266, 0.136807) −0.131474 1000 × 200 (0.623075, 0.138929) −0.066984
Results are comparable with numerical appeared in literature
Falling Rigid Object in Water Tank
−2 0 2
−3
−2
−1 0 1 2 3
Density air
water
−2 0 2
−3
−2
−1 0 1 2 3
Pressure
−2 0 2
−3
−2
−1 0 1 2 3
Volume fraction
t = 0ms
Falling Rigid Object in Water Tank
−2 0 2
−3
−2
−1 0 1 2 3
Density air
water
−2 0 2
−3
−2
−1 0 1 2 3
Pressure
−2 0 2
−3
−2
−1 0 1 2 3
Volume fraction
t = 1ms
Falling Rigid Object in Water Tank
−2 0 2
−3
−2
−1 0 1 2 3
Density air
water
−2 0 2
−3
−2
−1 0 1 2 3
Pressure
−2 0 2
−3
−2
−1 0 1 2 3
Volume fraction
t = 2ms
Falling Rigid Object in Water Tank
−2 0 2
−3
−2
−1 0 1 2 3
Density air
water
−2 0 2
−3
−2
−1 0 1 2 3
Pressure
−2 0 2
−3
−2
−1 0 1 2 3
Volume fraction
t = 3ms
Future Work
3D volume tracking method
General curvilinear grid system
Body-fitted grid for complicated geometries Low Mach number flow
Remove sound-speed stiffness by preconditioning techniques or pressure-based method
Include more physics towards real applications
Diffusion, phase transition, or elastic-plastic effect Hybrid surface-volume tracking algorithm for balance laws with interfaces & boundaries
Liquid Drop Problem (Revisit)
Tracking Capturing
air air
t = 0 t = 0
Liquid Drop Problem (Revisit)
Tracking Capturing
air air
t = 0.2s t = 0.2s
Liquid Drop Problem (Revisit)
Tracking Capturing
air air
t = 0.4s t = 0.4s
Liquid Drop Problem (Revisit)
Tracking Capturing
air air
t = 0.6s t = 0.6s
Liquid Drop Problem (Revisit)
Tracking Capturing
air air
t = 0.8s t = 0.8s
Liquid Drop Problem (Revisit)
Tracking Capturing
air air
t = 1s t = 1s
Thank You
References
(JCP 1998) An efficient shock-capturing algorithm for compressible multicomponent problems
(JCP 1999, 2001) A fluid-mixture type algorithm for
compressible multicomponent flow with van der Waals (Mie-Grüneisen) equation of state
(JCP 2004) A fluid-mixture type algorithm for barotropic two-fluid flow Problems
(JCP 2006) A wave-propagation based volume tracking method for compressible multicomponent flow in two space dimensions
(Shock Waves 2006) A volume-fraction based algorithm for hybrid barotropic & non-barotropic two-fluid flow
problems
Thermodynamic Stability
Fundamental derivative of gas dynamics
G = −V 2
(∂2p/∂V 2)S
(∂p/∂V )S , S : specific entropy Assume fluid state satisfy G > 0 for thermodynamic stability, i.e.,
(∂2p/∂V 2)S > 0 & (∂p/∂V )S < 0 (∂2p/∂V 2)S > 0 means convex EOS
(∂p/∂V )S < 0 means real speed of sound, for c2 = ∂p
∂ρ
S
= −V 2 ∂p
∂V
S
> 0
Homogeneous Flow Model (cont.)
Mie-Grüneisen EOS: pk = pref(ρk) + ρkΓ(ρk)[ek − eref(ρk)]
ρe =
2
X
k=1
αkρkek =
2
X
k=1
αk
p − pref(ρk)
Γ(ρk) + ρkeref(ρk)
=⇒
p =
"
ρe −
2
X
k=1
αk
−pref(ρk)
Γ(ρk) + ρkeref(ρk)
#
Mie-Grüneisen Equations of State
(pref, eref) lies along an isentrope
1. Jones-Wilkins-Lee EOS for gaseous explosives
Γ(V ) = γ − 1, V = 1/ρ eref(V ) = e0 + A V0
R1 exp −R1V V0
+ B V0
R2 exp −R2V V0
pref(V ) = p0 + A exp −R1V V0
+ B exp −R2V V0
2. Cochran-Chan EOS for solid explosives
Γ(V ) = γ − 1
eref(V ) = e0 + −A V0 1 − E1
"
V V0
1−E1
− 1
#
+ B V0 1 − E2
"
V V0
1−E2
− 1
#
pref(V ) = p0 + A V −E1
− B V −E2
Mie-Grüneisen EOS (cont.)
(pref, eref) lies along a Hugoniot locus
Assume linear shock speed us & particle velocity up
us = c0 + s up We may derive the relations
Γ(V ) = Γ0 V V0
α
, Γ0 = γ − 1 pref(V ) = p0 + c02(V0 − V )
[V0 − s(V0 − V )]2 eref(V ) = e0 + 1
2 pref(V ) + p0 (V0 − V )
Material Quantities for Model EOS
JWL EOS ρ0(kg/m3) A(GPa) B(GPa) R1 R2 Γ
TNT1 1630 371.2 3.23 4.15 0.95 0.30
TNT2 1630 548.4 9.375 4.94 1.21 1.28
Water 1004 1582 −4.67 8.94 1.45 1.17
CC EOS ρ0(kg/m3) A(GPa) B(GPa) E1 E2 Γ
TNT 1840 12.87 13.42 4.1 3.1 0.93
Copper 8900 145.67 147.75 2.99 1.99 2
Shock EOS ρ0(kg/m3) c0(m/s) s Γ0 α
Aluminum 2785 5328 1.338 2.0 1
Copper 8924 3910 1.51 1.96 1
Molybdenum 9961 4770 1.43 2.56 1
MORB 2660 2100 1.68 1.18 1
Water 1000 1483 2.0 2.0 10−4