• 沒有找到結果。

# moving-mesh relaxation scheme for

N/A
N/A
Protected

Share "moving-mesh relaxation scheme for"

Copied!
37
0
0

(1)

### compressible two-phase barotropicflow with cavitation

Keh-Ming Shyue

Department of Mathematics National Taiwan University

Taiwan

(2)

### Outline

Introduce simple relaxation model for compressible 2-phase barotropic flow with & without cavitation

Describe finite-volume method for solving proposed model numerically on fixed mapped (body-fitted) grids Discuss adaptive moving mesh approach for efficient improvement of numerical resolution

Show sample results & discuss future extensions

(3)

### Conventional barotropic 2-phase model

Assume homogeneous equilibrium barotropic flow with single pressure p & velocity ~u in 2-phase mixture region

Equations of motion

t1ρ1) + ∇ · (α1ρ1~u) = 0 (Continuity α1ρ1)

t2ρ2) + ∇ · (α2ρ2~u) = 0 (Continuity α2ρ2)

t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0 (Momentum ρ~u) Saturation condition for volume fraction αk, k = 1, 2

α1 + α2 = 1 =⇒ α1ρ1

ρ1(p) + α2ρ2

ρ2(p) = 1

yielding scalar nonlinear equation for equilibrium p

(4)

### Equilibrium speed of sound

Non-monotonic sound speed ceq defined as 1

ρc2eq

= α1

ρ1c21 + α2

ρ2c22 (Wood’s formula)

exists in 2-phase region, air-water case shown below;

yielding numerical difficulty such as inaccurate wave transmission across diffused interface

0 0.2 0.4 0.6 0.8 1

0 500 1000 1500

α

c eq

(5)

### Equilibrium vs. frozen sound speed

Equilibrium (solid) & frozen (dashed) sound speeds, c2

f =

PYkc2k, in case of passive advection of air-water interface

(6)

### Equilibrium vs. frozen sound speed

When acoustic wave interacts with numerical diffusion zone, sound speeds difference leads to time delay τ of transmitted waves through interface

(7)

### Relaxation barotropic 2-phase model

To overcome these difficulties, we consider a relaxation

model proposed by Caro, Coquel, Jamet, Kokh, Saurel et al.

t1ρ1) + ∇ · (α1ρ1~u) = 0

t2ρ2) + ∇ · (α2ρ2~u) = 0

t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇ (α1p1 + α2p2) = 0

tα1 + ~u · ∇α1 = µ (p11) − p22)) (Transport α1) When taking infinite pressure relaxation µ → ∞, we have

p11) = p22) =⇒ p1

1ρ1 α1



− p2

 α2ρ2 1 − α1



= 0 yielding scalar nonlinear equation for volume fraction α1

(8)

### Relaxation barotropic model: Remarks

This model can be viewed as isentropic version of

relaxation model proposed by Saurel, Petitpas, Berry (JCP 2009, see below) for 2-phase flow

This model is hyperbolic & has monotonic sound speed c2f = P

Ykc2k

Cavitation is modeled as a simplified mechanical

relaxation process, occurring at infinite rate & not as a mass transfer process

i.e., cavitation pockets appear as volume fraction

increases for a small amount of gas present initially

(9)

### Relaxation 2-phase model

Saurel, Petitpas, Berry (JCP 2009)

t1ρ1) + ∇ · (α1ρ1~u) = 0

t2ρ2) + ∇ · (α2ρ2~u) = 0

t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇ (α1p1 + α2p2) = 0

t1ρ1e1) + ∇ · (α1ρ1e1~u) + α1p1∇ · ~u = −pµ (p¯ 1 − p2)

t2ρ2e2) + ∇ · (α2ρ2e2~u) + α2p2∇ · ~u = pµ (p¯ 1 − p2)

tα1 + ~u · ∇α1 = µ (p1 − p2)

¯

p = (p1Z2 + p2Z1)/(Z1 + Z2); Zk = ρkck

Model agrees with reduced 2-phase model of Kapila,

Menikoff, Bdzil, Son, Stewart (Phys. Fluid 2001) formally

(10)

### 1-phase barotropic cavitation models

Cutoff model

p =

(p(ρ) if ρ ≥ ρsat psat if ρ < ρsat Schmidt model

p =



p(ρ) if ρ ≥ ρsat

psat + pgl ln

 ρgc2gρlc2l l+α(ρgρl)) ρl(ρgc2gα(ρgc2gρlc2l ))



if ρ < ρsat

with pgl = ρgc2gρlc2lg − ρl)

ρ2gc2g − ρ2l c2l & ρ = αρg + (1 − α)ρl Modified Schmidt model & its variant

(11)

### Mapped grid method

We want to use finite-volume mapped grid approach to solve proposed relaxation model in complex geometry Assume mapped grids are logically rectangular & will

review method for hyperbolic system of conservation laws

tq + ∇ · f (q) = 0

i − 1 i − 1

i

i j

j

j + 1 j + 1

Cˆij Cij

ξ1

ξ2

mapping

∆ξ1

∆ξ2

logical grid physical grid

←−

x1 = x11, ξ2) x2 = x21, ξ2) x1

x2

(12)

### Mapped grid methods

On a curvilinear grid, a finite volume method takes Qn+1ij = Qnij− ∆t

κij∆ξ1

Fi+1 1

2,j − Fi−1 1 2,j



− ∆t κij∆ξ2

Fi,j+2 1 2

− Fi,j−2 1 2



∆ξ1, ∆ξ2 denote mesh sizes in ξ1- & ξ2-directions

κij = M(Cij)/∆ξ1∆ξ2 is area ratio between areas of grid cell in physical & logical spaces

Fi−1 1

2,j = γi−1

2,ji−1

2,j, Fi,j−2 1

2

= γi,j−1

2

i,j−1

2 are normal fluxes per unit length in logical space with γi−1

2,j = hi−1

2,j/∆ξ1 &

γi,j−1

2 = hi,j−1

2 /∆ξ2 representing length ratios

(13)

### Wave propagation method

First order wave propagation method devised by LeVeque on mapped grid is a Godunov-type finite volume method

Qn+1ij = Qnij− ∆t κij∆ξ1



A+1 ∆Qi−1

2,j + A1 ∆Qi+1

2,j



∆t κij∆ξ2



A+2 ∆Qi,j−1

2 + A2 ∆Qi,j+1

2



with right-, left-, up-, & down-moving fluctuations A+1 ∆Qi−1

2,j, A1 ∆Qi+1

2,j, A+2 ∆Qi,j−1

2, & A2 ∆Qi,j+1

2 that are entering into grid cell

To determine these fluctuations, one-dimensional Riemann problems in direction normal to cell edges are solved

(14)

### Wave propagation method

Speeds & limited versions of waves are used to calculate second order correction terms as

Qn+1ij := Qn+1ij 1 κij

∆t

∆ξ1

Fei+1 1

2,j Fe1

i−12,j



1 κij

∆t

∆ξ2

Fei,j+2 1

2

Fe2

i,j−12



For example, at cell edge (i − 12, j) correction flux takes e

F1

i−12,j = 1 2

Nw

X

k=1

λ1,ki−1

2,j

1 − κ ∆t

i−12,j∆ξ1

λ1,ki−1

2,j

! f W1,k

i−21,j

κi−1

2,j = (κi−1,j + κij)/2. To aviod oscillations near

discontinuities, a wave limiter is applied leading to limited waves Wf

(15)

### Wave propagation method

Transverse wave propagation is included to ensure second order accuracy & also improve stability that A±1 ∆Qi−1

2,j are each split into two transverse fluctuations: up- &

down-going A±2 A+1 ∆Qi−1

2,j & A±2 A1 ∆Qi−1

2,j, at each cell edge

This method can be shown to be conservative & stable

under a variant of CFL (Courant-Friedrichs-Lewy) condition of form

ν = ∆t max

i,j,k

λ1,ki−1

2,j

κip,j∆ξ1

,

λ2,ki,j−1

2

κi,jp∆ξ2

 ≤ 1,

ip = i if λ1,k

i−12,j > 0 & i − 1 if λ1,k

i−12,j < 0

(16)

### Extension to moving mesh

To extend mapped grid method to solution adaptive moving grid method one simple way is to take approach proposed by

H. Tang & T. Tang, Adaptive mesh methods for one- and two-dimensional hyperbolic conservation laws, SIAM J.

Numer. Anal., 2003

In each time step, this moving mesh method consists of three basic steps:

(1) Mesh redistribution

(2) Conservative interpolation of solution states (3) Solution update on a fixed mapped grid

(17)

### Mesh redistribution scheme

Winslow’s approach (1981)

Solve ∇ · (D∇ξj) = 0, j = 1, . . . , Nd

for ξ(x). Coefficient D is a positive definite matrix which may depend on solution gradient

Variational approach (Tang & many others)

Solve ξ · (D∇ξxj) = 0, j = 1, . . . , Nd for x(ξ) that minimizes “energy” functional

E(x(ξ)) = 1 2

Z

Nd

X

j=1

Tξ D∇xj

Lagrangian (ALE)-type approach (e.g., CAVEAT code)

(18)

### Conservative interpolation

Numerical solutions need to be updated conservatively, i.e.

XM Ck+1

Qk+1 = X

M Ck Qk

after each redistribution iterate k. This can be done by Finite-volume approach

M Ck+1

Qk+1 = M Ck

Qk

Ns

X

j=1

hjG˘j, G˘ = ( ˙x · n)Q

Geometric approach

"

X

S

M Cpk+1 ∩ Spk#

Qk+1C = X

S

M Cpk+1 ∩ Spk QkS

Cp, Sp are polygonal regions occupied by cells C & S

(19)

### Interpolation-free moving mesh

If we want to derive an interpolation-free moving mesh method, one may first consider coordinate change of equations via (x, t) 7→ (ξ, t), yielding transformed

conservation law as

tq˜ + ∇ξ · ˜f = G

˜

q = Jq, f˜j = J (q ∂tξj + ∇ξj · f) , J = det(∂ξ/∂x)1 G = q 

tJ + ∇ξ · (J∂tξj) +

XN j=1

fjξ · J∂xjξk

= 0 (if GCL & SCL are satisfied)

Numerical method can be devised easily to solve these equations

(20)

### Relaxation solver on moving meshes

In each time step, our numerical method for solving

barotropic 2-phase flows on a moving mesh consists of following steps:

(1) Moving mesh step

Determine cell-interface velocity & cell-interface location in physical space over a time step

(2) Frozen step µ → 0

Solve homogeneous part of relaxation model on a

moving mapped grid over same time step as in step 1 (3) Relaxation step µ → ∞

Solve model system with only source terms in infinite relaxation limit

(21)

### Water-vapor cavitation

Initially, in closed shock tube, flow is homogeneous (contains α = 106 gas in bulk liquid) at standard

atmospheric condition & exists interface separating flow with opposite motion (u = 100 m/s)

Result in pressure drop & formation of cavitation zone in middle; shocks form also from both ends

Eventually, shock-cavitation collision occurs

−u u

← Inteface

(22)

(23)

### Water-vapor cavitation

Physical grid in x-t plane

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1x 10−3

time

x

(24)

### Oscillating water column

Initially, in closed shock tube, water column moves at u = 1 from left to right, yielding air compression at right

& air expansion at left

Subsequently, pressure difference built up across water column resulting deceleration of column of water to

right, makes a stop, & then acceleration to left; a reverse pressure difference built up across water column redirecting flow from left to right again

Eventually, water column starts to oscillate

air

air water

u →

(25)

### Oscillating water column

Physical grid in x-t plane

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

time

x

(26)

(27)

### Oscillating water column

Time evolution of pressure at left & right boundaries

0 1 2 3 4 5 6 7 8 9 10

0.8 1 1.2 1.4 1.6 1.8

p L

0 1 2 3 4 5 6 7 8 9 10

0.8 1 1.2 1.4 1.6 1.8

p R

mmesh mesh0

mmesh mesh0

time

(28)

### Tait equation of state & parameters

Each fluid phase (liquid & gas) satisfies Tait equation of state

pk(ρ) = (p0k + Bk)

 ρ ρ0k

γk

− Bk for k = 1, 2. with parameters for liquid phase as

(γ, B, ρ0, p0)1 = 

7, 3000 bar, 103 kg/m3, 1 bar while parameters for gas phase as

(γ, B, ρ0, p0)2 = 

1.4, 0, 1 kg/m3, 1 bar

(29)

### Supersonic flow over a bluntbody

Formation of cavitation zone

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Volume fraction

0 0.5 1 1.5 2 2.5 x 108

Pressure

(30)

### Underwater explosion

Contours of density and pressure at selected times

Density t=0.24ms

air

water

Pressure

t=0.4ms

(31)

### Underwater explosion

Contours of density and pressure at selected times

t=0.8ms

t=1.2ms

(32)

### Underwater explosion

Physical grid at selected times

−1 0 1

−1

−0.5 0 0.5

−1 0 1

−1

−0.5 0 0.5

t = 0.2ms t = 0.4ms

−1

−0.5 0 0.5

−1

−0.5 0 0.5

t = 0.8ms t = 1.2ms

(33)

### Final remarks

Show preliminary results obtained using relaxation moving mesh methods for barotropic 2-phase flow problem

Cavitation problems are often occurred in low Mach

scenario & so suitable fixed up such as preconditioning

& others are necessary for solution accuracy improvement

Extension to non-barotropc cavitation with phase transition

. . .

(34)

## Thank you

(35)

### Reduced 2-phase flow model

Reduced 2-phase flow model of Kapila et al. is zero-order approximation of Baer-Nunziato equations with stiff

mechanical relaxation that takes

t1ρ1) + ∇ · (α1ρ1~u) = 0

t2ρ2) + ∇ · (α2ρ2~u) = 0

t (ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0

t (ρE) + ∇ · (ρE~u + p~u) = 0

tα1 + ~u · ∇α1 = ρ2c22 − ρ1c21 P2

k=1 ρkc2kk ∇ · ~u

(36)

### Baer-Nunziato Two-Phase Flow Model

Baer & Nunziato (J. Multiphase Flow 1986) (α1ρ1)t + ∇ · (α1ρ1~u1) = 0

1ρ1~u1)t + ∇ · (α1ρ1~u1 ⊗ ~u1) + ∇(α1p1) = p0∇α1 + λ(~u2 − ~u1) (α1ρ1E1)t + ∇ · (α1ρ1E1~u1 + α1p1~u1) = p02)t + λ~u0 · (~u2 − ~u1) (α2ρ2)t + ∇ · (α2ρ2~u2) = 0

2ρ2~u2)t + ∇ · (α2ρ2~u2 ⊗ ~u2) + ∇(α2p2) = p0∇α2 − λ(~u2 − ~u1)

2ρ2E2)t + ∇ · (α2ρ2E2~u2 + α2p2~u2) = −p02)t − λ~u0 · (~u2 − ~u1) (α2)t + ~u0 · ∇α2 = µ (p2 − p1)

αk = Vk/V : volume fraction (α1 + α2 = 1), ρk: density,

~uk: velocity, pk: pressure, Ek = ek + ~u2k/2: specific total energy, e : specific internal energy, k = 1, 2

(37)

### Baer-Nunziato Model (Cont.)

p0 & ~u0: interfacial pressure & velocity Baer & Nunziato (1986)

p0 = p2, ~u0 = ~u1 Saurel & Abgrall (1999)

p0 = P2

k=1 αkpk, ~u0 = P2

k=1 αkρk~uk

P2

k=1 αkρk λ & µ (> 0): relaxation parameters that determine rates at which velocities and pressures of two phases reach

equilibrium

With the proposed model equations, accurate results can be obtained on a mapped grid using a standard method, such as the high-resolution wave- propagation algorithm for a

Shock diffraction in 2-phase mixture Cavitating Richtmyer-Meshkov instability High pressure fuel injector.. Homogeneous relaxation

Review a high-resolution wave propagation method for solving hyperbolic problems on mapped grids (which is basic integration scheme implemented in CLAWPACK) Describe

Have shown results in 1 , 2 &amp; 3 D to demonstrate feasibility of method for inviscid compressible flow problems. Department of Applied Mathematics, Ta-Tung University, April 23,

The disadvantage of the inversion methods of that type, the encountered dependence of discretization and truncation error on the free parameters, is removed by

In this paper, we have studied a neural network approach for solving general nonlinear convex programs with second-order cone constraints.. The proposed neural network is based on

Since the generalized Fischer-Burmeister function ψ p is quasi-linear, the quadratic penalty for equilibrium constraints will make the convexity of the global smoothing function

Now we assume that the partial pivotings in Gaussian Elimination are already ar- ranged such that pivot element a (k) kk has the maximal absolute value... The growth factor measures