DNS of Multiphase Flows DNS of Multiphase Flows Direct Numerical Simulations of Multiphase The boundary Flows-3 conditions A Simple Solver for Variable Density Flow (3 of 3) Gretar Tryggvason DNS of Multiphase Flows DNS of Multiphase Flows Tangent velocity Boundary Conditions v i,2 ! v i+1,2 ! Boundary Conditions for the velocity: Velocity of wall is given, u wall (no-slip) u i,2 ! p i+1,2 ! Interpolate linearly p i,2 Here we take the domain to be a rectangle with ! prescribed velocities at all boundaries. u wall = u i, 1 + u i, 2 v i,1 ! v i+1,1 ! 2 The boundaries coincide with the edges of the Solve for the “ghost” velocity pressure control volumes so the normal velocity is p i,1 ! u i,1 ! p i+1,1 ! given where it is needed. u i, 1 = 2 u wall − u i, 2 If the wall velocity is zero: The implementation of the tangent velocity is slightly u wall = 0 then the ghost more complex and we will use “ghost” points outside % tangential velocity at boundaries velocity is a reflection of u(1:nx+1,1)=2*usouth-u(1:nx+1,2); the domain to impose the correct tangential velocity u(1:nx+1,ny+2)=2*unorth-u(1:nx+1,ny+1); the interior velocity v(1,1:ny+1)=2*vwest-v(2,1:ny+1); v(nx+2,1:ny+1)=2*veast-v(nx+1,1:ny+1); u i, 1 = � u i, 2 DNS of Multiphase Flows DNS of Multiphase Flows Boundary Conditions for the A simple way to implement the pressure boundary v i,j +1/2 ! v i +1 ,j +1/2 ! pressure. Use the normal conditions is to set the value of the density in the ghost cell velocity given at the to a large value and the velocity to the normal velocity boundary when we p i, j ! p i +1, j ! u i +1/2 ,j ! U b,j = u i -1/2 ,j ! substitute the correction p α +1 = i,j velocity into the mass ⌘# − 1 " 1 1 1 1 1 1 ⇣ ⌘ ⇣ v i,j -1/2 ! v i +1 ,j -1/2 ! β + + + conservation equation: ρ n +1 i +1 ,j + ρ n +1 ρ n +1 + ρ n +1 ρ n +1 i,j +1 + ρ n +1 ρ n +1 + ρ n +1 ∆ x 2 ∆ y 2 i,j i,j i − 1 ,j i,j i,j i,j − 1 " p α +1 p α +1 p α p α 1 1 ⇣ ⌘ ⇣ ⌘ u n +1 v n +1 i,j +1 / 2 � v n +1 i +1 ,j i − 1 ,j i,j +1 i,j − 1 i +1 / 2 ,j � U b,j + + + i,j − 1 / 2 ρ n +1 i +1 ,j + ρ n +1 ρ n +1 + ρ n +1 ρ n +1 i,j +1 + ρ n +1 ρ n +1 + ρ n +1 ∆ x 2 ∆ y 2 + = 0 i,j i,j i − 1 ,j i,j i,j i,j − 1 ∆ x ∆ y ⌘# ⇣ u ∗ i +1 / 2 ,j − u ∗ v ∗ i,j +1 / 2 − v ∗ − 1 i − 1 / 2 ,j i,j − 1 / 2 + + (1 − β ) p α The pressure equation at the boundary is therefore: i,j 2 ∆ t ∆ x ∆ y ! ! 1 1 p i +1 ,j � p i,j p i,j +1 � p i,j � p i,j � p i,j − 1 + ρ n +1 i +1 ,j + ρ n +1 ρ n +1 i,j +1 + ρ n +1 ρ n +1 + ρ n +1 ∆ x 2 ∆ y 2 i,j i,j i,j i,j − 1 rt=r; lrg=1000; ⇣ u ∗ i +1 / 2 ,j � U b,j v ∗ i,j +1 / 2 � v ∗ 1 ⌘ i,j − 1 / 2 = + rt(1:nx+2,1)=lrg; rt(1:nx+2,ny+2)=lrg; 2 ∆ t ∆ x ∆ y rt(1,1:ny+2)=lrg; rt(nx+2,1:ny+2)=lrg;
DNS of Multiphase Flows DNS of Multiphase Flows Advecting the Density: Centered Differences with Diffusion ∂ρ ∂ρ Write as: ∂ t = � u · r ρ ∂ t = �r · ( ρ u ) Advecting the Add a diffusion term ∂ρ ∂ t = �r · ( ρ u ) + µ 0 r 2 ρ Density Discrete version i,j � ∆ t ρ i +1 ,j + ρ i,j ρ i − 1 ,j + ρ i,j ⇣ ⌘ ρ n +1 = ρ n u i +1 / 2 ,j � u i − 1 / 2 ,j i,j ∆ x 2 2 � ∆ t ρ i,j +1 + ρ i,j ρ i,j − 1 + ρ i,j ⇣ ⌘ v i,j +1 / 2 � u i,j − 1 / 2 ∆ y 2 2 + µ 0 ∆ t ∆ x 2 ( ρ i +1 ,j � 2 ρ i,j + ρ i − 1 ,j ) + µ 0 ∆ t ∆ y 2 ( ρ i,j +1 � 2 ρ i,j + ρ i,j − 1 ) DNS of Multiphase Flows DNS of Multiphase Flows Discrete version i,j � ∆ t ρ i +1 ,j + ρ i,j ρ i − 1 ,j + ρ i,j ⇣ ⌘ ρ n +1 = ρ n u i +1 / 2 ,j � u i − 1 / 2 ,j i,j ∆ x 2 2 � ∆ t ρ i,j +1 + ρ i,j ρ i,j − 1 + ρ i,j ⇣ ⌘ v i,j +1 / 2 � u i,j − 1 / 2 ∆ y 2 2 The full method + µ 0 ∆ t ∆ x 2 ( ρ i +1 ,j � 2 ρ i,j + ρ i − 1 ,j ) + µ 0 ∆ t ∆ y 2 ( ρ i,j +1 � 2 ρ i,j + ρ i,j − 1 ) %=====ADVECT DENSITY using centered difference plus diffusion ======== ro=r; for i=2:nx+1,for j=2:ny+1 r(i,j)=ro(i,j)-(0.5*dt/dx)*(u(i,j)*(ro(i+1,j)+ro(i,j))-u(i-1,j)*(ro(i-1,j)+ro(i,j)) )... -(0.5*dt/dy)*(v(i,j)*(ro(i,j+1)+ro(i,j))-v(i,j-1)*(ro(i,j-1)+ro(i,j)) )... +(m0*dt/dx/dx)*(ro(i+1,j)-2.0*ro(i,j)+ro(i-1,j))... +(m0*dt/dy/dy)*(ro(i,j+1)-2.0*ro(i,j)+ro(i,j-1)); end,end DNS of Multiphase Flows DNS of Multiphase Flows Discretization in time Algorithm Initialize parameters and arrays, set Initial field given time step and initial density 1. Set velocity at ghost points and update the marker function to find new density and viscosity for is=1:nstep Determine u , v set boundary conditions for 2. Find a temporary velocity using the advection and the boundary conditions tangential velocity (ghost points) diffusion terms only: Advect density ! 1 Advect the density ρ n h u n � � A n h + g h + D n � u ∗ h = h + ∆ t h ρ n +1 ! h 1 ρ n ⇣ h u n � A n h + g h + D n � � u ∗ h = h + ∆ t Find predicted velocity 3. Find the pressure needed to make the velocity h ρ n +1 ⇣ ! h h h field incompressible ! ! h h 1 = 1 Solve for pressure using SOR r h · r h p h ∆ t r h · u ∗ ! 1 = 1 h ρ n +1 r h p h ∆ t r h · u ∗ r h · h h ρ n +1 Find the projected velocity h h � ∆ t r h p h u n +1 by adding the pressure gradient = u ∗ 4. Correct the velocity by adding the pressure gradient: h ρ n +1 h PLOT h � ∆ t r h p h t = t + Δ t u n +1 = u ∗ end h ρ n +1 h
Recommend
More recommend