Interface Treating Method for Multi-Medium Flow on Adaptive - PowerPoint PPT Presentation

� Sino-German Symposium on Modern Numerical Methods for Compressible Fluid Flows and Related Problems, May 21-27, 2014, Beijing Interface Treating Method for Multi-Medium Flow on Adaptive Triangular Meshes Chunwu Wang Joint work with Professor Chi-Wang Shu College of Science Nanjing University of Aeronautics and Astronautics Nanjing 210016 , P.R.China Email:wangcw@nuaa.edu.cn • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Outline � Introduction � Equations � The RGFM on triangular meshes � Some numerical results � Conclusions • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

1. Introduction • The RKDG method is used to discretize of the Euler equation. Advantages: Easy to handle the complicated fluid boundary, compactness,...... Disadvantages: More degrees of freedom, not cost efficient,...... • Triangular mesh: Easy to decompose the fluid field with complicated boundary. Motivation • Present an interface treating method for multi-medium flow on triangular meshes, easy to handle complicated geometry and fluid boundary. • Couple the adaptive moving mesh method to the interface treating method to reduce the numerical and conservative errors and improve the solution resolution. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

The ghost fluid method(GFM) 1D Original GFM( R. Fedkiw et al, JCP,1999 ). 1D New GFM( R. Fedkiw et al, JCP,2002 ). • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

� � � � � The modified GFM( T.G.Liu et al, JCP, 2003 ). interface B Medium 1 Medium 2 i 1 i 2 i 3 p I A � real nodes u � ghost nodes I d IL Medium 1 Medium 2 i 2 i 1 i The 1D and 2D Real GFM( C.W.Wang et al, SISC,2006 ). • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

2. Equations Euler equation (Chi-Wang Shu et al, JCP, 1998) The two-dimensional system for a compressible fluid can be written as follows: ∂ � ∂ t + ∂ � ∂ x + ∂ � U F G ∂ y = 0 where � U =[ ρ , ρ u , ρ v , E ] , U ) =[ ρ u , ρ u 2 + p , ρ uv , ( E + p ) u ] T , � F ( � U ) =[ ρ v , ρ uv , ρ v 2 + p , ( E + p ) v ] T . � G ( � E = ρ e + 1 2 ρ ( u 2 + v 2 ) For closure of the system, the equation of state (EOS) is required. In the present work, our in- terest centered on the compressible gas and water media. The EOS for gases or water medium can be written uniformly as p = ( γ − 1 ) ρ e − γ B where γ and B are treated as fluid constants. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Level set equation (Chi-Wang Shu et al, JCP, 2007) To track the moving fluid interface, we employ the level set technique. The level set equation can be written as ∂φ ∂ t + � u · ∇ φ = 0 . • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

3. The RGFM on triangular meshes Case 1 There is one neighbor of triangle 1 in another medium. Case 2 There are two neighbors of triangle 1 in another medium. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

The unit outward normal The normal is needed in the construction of the Riemann problem. In this work, the normal direction of the level set function at the barycenter of triangle ∆ l can be approximated by using the Green formula as � � ∂ P ( x , y ) − ∂ Q ( x , y ) dxdy = Q ( x , y ) dx + P ( x , y ) dy ∂ x ∂ y △ l ∂ △ l If P ( x , y ) = φ ( x , y ) , Q ( x , y ) = 0, we have � � ∂φ ( x , y ) dxdy = φ ( x , y ) dy ∂ x △ l ∂ △ l So φ x ( x , y ) can be approximated at the barycenter of triangle △ l as � � ∂ △ l φ ( x , y ) dy 1 φ x ( x b , y b ) = � = φ ( x , y ) dy |△ l | △ l dxdy ∂ △ l Similar procedures can be carried out to obtain φ y ( x b , y b ) . So the unit normal can be written as φ x ( x b , y b ) φ y ( x b , y b ) � � N = ( , ) φ 2 x ( x b , y b )+ φ 2 φ 2 x ( x b , y b )+ φ 2 y ( x b , y b ) y ( x b , y b ) • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Redefinition of the fluid velocity When large gradient appears in the velocity field, the level set contours may be severely dis- torted or even the code may collapse. To overcome this difficulty, the extension velocity is applied by solving a convection equation ∂ I ∂ t + N · ∇ I = 0 . to extrapolate the interface velocity to the whole computational domain. By solving the level set function with the extension velocity, the uniformly distributed level set contour can be obtained. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

The moving triangular meshes (H. Z. Tang et al, JSC, 2008) A one-to-one coordinate transformation from the logical or computational domain to the phys- ical domain is denoted by x = x ( ξ ) , ξ ∈ Ω l The map is to find the minimizer of the following functional � 2 E ( x ) = 1 ( � ∇ x i ) T G i ( � � ∑ ∇ x i ) d ξ 2 Ω l i = 1 By using the choice of Winslow, we deduce the Euler-Lagrange equations of the functional ∇ x · ( ω � � ∇ x ) = 0 • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Discretize the equation de solve it with a relaxed iteration method as ℵ ( i ) ℵ ( i ) W i j x ν ∑ ∑ x i = � i j / W i j j = 1 j = 1 x ν + 1 x i +( 1 − µ i ) x ν = µ i � i i where W i j = ∆ τ | V i | ω i j | ˆ l i j | / | ξ i j ξ i | ℵ ( i ) ∑ µ i = max ( W i j , σ ) j = 1 The monitor function is taken as � ω 2 ω 2 ω 2 ω E i = 1 + α ρ ¯ E i ( β ρ , ρ )+ α p ¯ E i ( β p , p )+ α q ¯ E i ( β q , q ) where | ∇ ξ ρ | ω E i ( β ρ , ρ ) = : min ( 1 , ¯ β ρ max E i {| ∇ ξ ρ |} ) � u 2 + v 2 q = • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Remarks • It is necessary to avoid the mesh moving across the inter- face which may lead to the mixed fluid around the inter- face. • It is worth noting that the fluid interface does not move physically in the iterative mesh redistribution. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Update the solution After each mesh iterative step, we need to update the approximate solutions U and φ on the new mesh from the old mesh . When the high-resolution conservative interpolation formulas for single-medium flow are di- rectly applied to updating the solutions U of multi-medium or multi-phase compressible flows, numerical inaccuracies will occur at the material interfaces due to the mesh-redistribution. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

1 case in the RKDG method, the approximate solution u h ( x , y , t ) inside the triangle K For the P can be written as 3 ∑ u i ( t ) φ i ( x , y ) u h ( x , y , t ) = i = 1 where the degrees of freedom u i ( t ) are values of the numerical solution at the midpoints of edges. The basis function φ i ( x , y ) = a i + b i x + c i y , i = 1 , 2 , 3 and   1 , ifi = j φ i ( x j , y j ) =  0 , otherwise The L 2 projection method is used to remap the numerical solution to the new mesh from the old mesh � � h ( x , y ) φ i ( x , y ) dxdy = ∑ K u N u O h ( x , y ) φ i ( x , y ) dxdy K i i • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

4. Some numerical results Problem 1. Sod problem   ( 1 , 0 , 0 , 1 , 1 . 4 , 0 ) , x < 0 . 5 ( ρ , u , v , p , γ , B ) =  ( 0 . 125 , 0 , 0 , 0 . 1 , 1 . 2 , 0 ) , x ≥ 0 . 5 The positive parameters in the monitor function are respectively taken as α ρ = 20 , α p = 20 , α q = 0 , β ρ = 0 . 1 , β p = 0 . 1 , β q = 0 . 1 , • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Problem 2. gas bubble explosion problem The initial interface lies at the middle of the computational domain. The initial states on the left and right sides of the interface are taken as   ( 1 . 5 , 0 , 0 , 100 , 1 . 4 , 0 ) , inside the bubble , ( ρ , u , v , p , γ , B ) =  ( 5 , 0 , 0 , 1 , 1 . 67 , 0 ) , outside of the bubble . α ρ = 50 , α p = 0 , α q = 0 , β ρ = 0 . 1 , β p = 0 . 1 , β q = 0 . 1 , • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

Problem 4. Shock bubble interaction problem The gas shock with strength 2 . 35 impact on the cylindrical water bubble. The pre- and post- shock states of the incident shock wave are   ( 1 , 0 , 0 , 1 , 2 . 8 , 3306 ) , inside the water cylinder , ( ρ , u , v , p , γ , B ) =  ( 0 . 001 , 0 , 0 , 1 , 1 . 4 , 0 ) , outside of the water cylinder , Schematics of the computational domain α ρ = 10 , α p = 0 , α q = 30 , β ρ = 0 . 1 , β p = 0 . 1 , β q = 0 . 1 , • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit

10 10 10 8 8 8 6 6 6 4 4 4 2 2 2 0 0 0 0 0 2 2 4 4 6 6 8 8 10 10 0 2 4 6 8 10 10 10 10 8 8 8 6 6 6 4 4 4 2 2 2 0 0 0 0 2 2 4 4 6 6 8 8 10 10 0 0 2 4 6 8 10 Pressure contours at time t = 0 . 02 , t = 0 . 04,respectively. • First • Prev • Next • Last • Go Back • Full Screen • Close • Quit