Surface a small segment of the front. For two-dimensional flows we - PDF document


 DNS of Multiphase Flows — Simple Front Tracking DNS of Multiphase Flows — Simple Front Tracking Direct Numerical Simulations of In this lecture we will complete our code by: Multiphase • Adding surface tension Flows-6 
 • Allowing the viscosities in the different fluids to be different • Make the time integration higher order Surface Tension, Unequal Viscosities, and Higher The resulting code is a fully functional one, and allows us to simulate simple multi fluid problems Order Time Integration Gretar Tryggvason DNS of Multiphase Flows DNS of Multiphase Flows The surface tension at a point is given by t e f σ = σκ n We need, however, the total force on − t s n l Surface a small segment of the front. For two-dimensional flows we can use Tension the definition of the curvature s is the coordinate along the interface κ n = ∂ t ∂ s The force on a small element l is Z Z Z ∂ t δ f σ f σ ds = l = σκ n ds = σ ∂ sds = σ � � t e − t s ∆ s l ∆ s l ∆ s l Z Thus, we do not need to find the curvature, just the tangent vectors, which is generally much simpler DNS of Multiphase Flows DNS of Multiphase Flows � Here we take the The surface tension is distributed to the fixed grid. On the − t l − 1 / 2 t l +1 / 2 interface segment front the force is per unit area and on the grid the force is per around each interface � � unit volume. The total force is conserved, so that: − x l x l +1 point to consist of half x l − 1 Z Z f σ f σ s ds = v dv the distance to the ∆ s l δ V Z Z points on either side Z The total force on an interface segment is The tangents at the end of the segment are found using a � Z X centered approximation F σ f σ δ f σ l = s ds ≈ l ∆ s l l � � / ∆ s t l +1 / 2 = x l +1 − x l where the integral is over the part of the interface Z contributing to a given grid point and the sum is over the where p ( x l +1 � x l ) 2 + ( y l +1 � y l ) 2 ∆ s = interface segment that do. X The surface force acting on the interface segment around The total force at a given grid point is point l is: δ f σ � � l = σ t l +1 / 2 − t l − 1 / 2 Z F σ f σ v dv ≈ f σ i,j = i,j ∆ x ∆ y δ V

DNS of Multiphase Flows DNS of Multiphase Flows F σ i,j = F σ Using that the force at each grid point, At solid walls we prescribe symmetry boundary conditions so l transferred from the front, is: the net force at the wall would be tangent to the wall l w l δ f σ For weights based on bilinear interpolation, the force is Sum over elements i,j X f σ i,j = “close” to grid point i,j distributed to the four nearest grid points so no adjustment is ∆ x ∆ y l needed for the normal velocity. For the tangential component we need to add the contribution from the “ghost” part f σ l w l The weights are generally taken to be the same as i,j those used to interpolate the velocities to the front For the bottom boundary j = 2 On a staggered grid, the x -component of the surface force is ( f x ) i,j =2 = ( f x ) i,j =2 + ( f x ) i,j =1 distributed to the u -nodes and the y -component to the v - nodes. Similar adjustments are done j = 1 for the other boundaries The surface forces on the fixed grid are then added to the discrete Navier-Stokes equations DNS of Multiphase Flows DNS of Multiphase Flows %------------------ Find surface tension ----------------------- Code to find fx=zeros(nx+2,ny+2);fy=zeros(nx+2,ny+2); % Set fx & fy to zero for l=1:Nf+1, surface tension ds=sqrt((xf(l+1)-xf(l))^2+(yf(l+1)-yf(l))^2); tx(l)=(xf(l+1)-xf(l))/ds; and distribute it ty(l)=(yf(l+1)-yf(l))/ds; % Tangent vectors end to the fixed grid tx(Nf+2)=tx(2);ty(Nf+2)=ty(2); Unequal for l=2:Nf+1 % Distribute to the fixed grid nfx=sigma*(tx(l)-tx(l-1));nfy=sigma*(ty(l)-ty(l-1)); ip=floor(xf(l)/dx)+1; jp=floor((yf(l)+0.5*dy)/dy)+1; Viscosities ax=xf(l)/dx-ip+1; ay=(yf(l)+0.5*dy)/dy-jp+1; fx(ip,jp) =fx(ip,jp)+(1.0-ax)*(1.0-ay)*nfx/dx/dy; fx(ip+1,jp) =fx(ip+1,jp)+ax*(1.0-ay)*nfx/dx/dy; fx(ip,jp+1) =fx(ip,jp+1)+(1.0-ax)*ay*nfx/dx/dy; fx(ip+1,jp+1)=fx(ip+1,jp+1)+ax*ay*nfx/dx/dy; ip=floor((xf(l)+0.5*dx)/dx)+1; jp=floor(yf(l)/dy)+1; ax=(xf(l)+0.5*dx)/dx-ip+1; ay=yf(l)/dy-jp+1; fy(ip,jp) =fy(ip,jp)+(1.0-ax)*(1.0-ay)*nfy/dx/dy; fy(ip+1,jp) =fy(ip+1,jp)+ax*(1.0-ay)*nfy/dx/dy; fy(ip,jp+1) =fy(ip,jp+1)+(1.0-ax)*ay*nfy/dx/dy; fy(ip+1,jp+1)=fy(ip+1,jp+1)+ax*ay*nfy/dx/dy; end fx(1:nx+2,2)=fx(1:nx+2,2)+fx(1:nx+2,1); % Correct boundary fx(1:nx+2,nx+1)=fx(1:nx+2,nx+1)+fx(1:nx+2,nx+2); % values for the fy(2,1:ny+2)=fy(2,1:ny+2)+fy(1,1:ny+2); % surface force fy(nx+1,1:ny+2)=fy(nx+1,1:ny+2)+fy(nx+2,1:ny+2); % on the grid DNS of Multiphase Flows DNS of Multiphase Flows For unequal viscosities we need to use the full deformation Integrating around the control volume gives tensor. The viscous term is discretized by integrating over the boundaries of the velocity control volume ( D x ) n i +1 / 2 ,j = v i,j +1/2 ! v i +1 ,j +1/2 ! ( ! D h = 1 1 I µ ∂ u µ ∂ u ⇣ ⌘ ⇣ ⌘ r u + ( r u ) T � 2 i +1 ,j � 2 ∆ y � · n ds = µ ∆ x ∆ y ∂ x ∂ x i,j ,j ! p i,j ! u i +1/2 ,j ! V p i +1 ,j ! S ! ) ⇣ ∂ u ∂ y + ∂ v ⇣ ∂ u ∂ y + ∂ v ⌘ ⌘ " # + µ i +1 / 2 ,j +1 / 2 � µ ∆ x 2 ∂ u ∂ u ∂ y + ∂ v ds = 1 I ∂ x ∂ x i +1 / 2 ,j − 1 / 2 ∂ x ∂ x µ · n ds v i,j -1/2 ! v i +1 ,j -1/2 ! ∂ y + ∂ v ∂ u 2 ∂ v V S ∂ x ∂ y ! ! ! u i +1/2 ,j +1 ! u i -1/2 ,j +1 p i,j +1 +1 Using that on vertical boundaries the normal is ( D y ) n i,j +1 / 2 = ! v i,j +1/2 ! ( ! 1 ⇣ ∂ v ∂ x + ∂ u ⇣ ∂ v ∂ x + ∂ u +1/2 ⌘ ⌘ µ i +1 / 2 ,j +1 / 2 − µ ∆ y n = (0 , 1) n = (0 , � 1) Horizontal boundaries and ∆ x ∆ y ∂ y ∂ y i − 1 / 2 ,j +1 / 2 ! ) ! u i -1/2 ,j ! p i,j ! u i +1/2 ,j ! µ ∂ v µ ∂ v Vertical boundaries and n = (1 , 0) n = ( � 1 , 0) ⇣ ⌘ ⇣ ⌘ + 2 i,j +1 − 2 ∆ x ∂ y ∂ y i,j

DNS of Multiphase Flows DNS of Multiphase Flows Diffusion term: discrete term for the u-velocity Diffusion term: discrete term for the v-velocity ( D y ) n ( D x ) n i,j +1 / 2 = i +1 / 2 ,j = (⇣ ⇣ u n i +1 / 2 ,j +1 − u n v n i +1 ,j +1 / 2 − v n ( ⇣ u n i +3 / 2 ,j − u n ⇣ u n i +1 / 2 ,j − u n ⌘) 1 1 ⌘ i +1 / 2 ,j i,j +1 / 2 ⌘ i +1 / 2 ,j i − 1 / 2 ,j µ n 2 µ n − 2 µ n + i +1 / 2 ,j +1 / 2 i +1 ,j i,j ∆ x ∆ x ∆ x ∆ x ∆ y ∆ x ⇣ u n i +1 / 2 ,j +1 − u n v n i +1 ,j +1 / 2 − v n ⇣ u n i − 1 / 2 ,j +1 − u n v n i,j +1 / 2 − v n ⌘) ( + 1 i +1 / 2 ,j i,j +1 / 2 ⌘ i − 1 / 2 ,j i − 1 ,j +1 / 2 − µ n µ n + + i +1 / 2 ,j +1 / 2 i − 1 / 2 ,j +1 / 2 ∆ y ∆ x ∆ y ∆ y ∆ x ⇣ v n i,j +3 / 2 − v n ⇣ v n i,j +1 / 2 − v n ⇣ u n i +1 / 2 ,j − u n v n i +1 ,j − 1 / 2 − v n ⌘) ( ⌘) + 1 i,j +1 / 2 ⌘ i,j − 1 / 2 i +1 / 2 ,j − 1 i,j − 1 / 2 2 µ n − 2 µ n − µ n + i,j +1 i,j i +1 / 2 ,j − 1 / 2 ∆ y ∆ y ∆ y ∆ y ∆ x ! ! ! u i +1/2 ,j +1 ! v i,j +1/2 ! v i +1 ,j +1/2 ! u i -1/2 ,j +1 p i,j +1 +1 i +1 / 2 ,j +1 / 2 = 1 µ n 4( µ n i +1 ,j + µ n i +1 ,j +1 + µ n i,j +1 + µ n i,j ) i +1 / 2 ,j +1 / 2 = 1 µ n 4( µ n i +1 ,j + µ n i +1 ,j +1 + µ n i,j +1 + µ n i,j ) ! v i,j +1/2 ! ! p i,j ! u i +1/2 ,j ! p i +1 ,j ! +1/2 i +1 / 2 ,j − 1 / 2 = 1 i − 1 / 2 ,j +1 / 2 = 1 µ n 4( µ n i,j + µ n i,j +1 + µ n i − 1 ,j +1 + µ n µ n 4( µ n i +1 ,j − 1 + µ n i +1 ,j + µ n i,j + µ n i,j − 1 ) i − 1 ,j ) ! u i -1/2 ,j ! p i,j ! u i +1/2 ,j ! v i,j -1/2 ! v i +1 ,j -1/2 ! DNS of Multiphase Flows DNS of Multiphase Flows New Code for the diffusion terms, allowing unequal viscosities After the density is updated, we need to set the viscosity for i=2:nx,for j=2:ny+1 % Temporary u-velocity-viscosity as a function of the marker function ut(i,j)=ut(i,j)+(2.0/(r(i+1,j)+r(i,j)))*dt*(... +(1./dx)*2.*(m(i+1,j)*(1./dx)*(u(i+1,j)-u(i,j)) - ... m(i,j) *(1./dx)*(u(i,j)-u(i-1,j)) ) ... for i=1:nx+2,for j=1:ny+2 % Update the viscosity +(1./dy)*( 0.25*(m(i,j)+m(i+1,j)+m(i+1,j+1)+m(i,j+1))* ... m(i,j)=m1+(m2-m1)*chi(i,j); ((1./dy)*(u(i,j+1)-u(i,j)) + (1./dx)*(v(i+1,j)-v(i,j)) ) - ... 0.25*(m(i,j)+m(i+1,j)+m(i+1,j-1)+m(i,j-1))* ... end,end ((1./dy)*(u(i,j)-u(i,j-1))+ (1./dx)*(v(i+1,j-1)- v(i,j-1))) ) ) ; end,end For a simple code as the one developed here we could for i=2:nx+1,for j=2:ny % Temporary v-velocity-viscosity vt(i,j)=vt(i,j)+(2.0/(r(i,j+1)+r(i,j)))*dt*(... use the density as a marker function. However, using a +(1./dx)*( 0.25*(m(i,j)+m(i+1,j)+m(i+1,j+1)+m(i,j+1))* ... separate marker function and then set both density and ((1./dy)*(u(i,j+1)-u(i,j)) + (1./dx)*(v(i+1,j)-v(i,j)) ) - ... viscosity as a function of the marker function is more 0.25*(m(i,j)+m(i,j+1)+m(i-1,j+1)+m(i-1,j))* ... ((1./dy)*(u(i-1,j+1)-u(i-1,j))+ (1./dx)*(v(i,j)- v(i-1,j))) )... general and allows us, for example, to examine flows of +(1./dy)*2.*(m(i,j+1)*(1./dy)*(v(i,j+1)-v(i,j)) - ... two fluid with the same density but different viscosities m(i,j) *(1./dy)*(v(i,j)-v(i,j-1)) ) ) ; end,end DNS of Multiphase Flows DNS of Multiphase Flows A widely used third order Runge-Kutta (RK-3) method. To integrate du dt = L ( u ) Higher Order in in time, RK-3 involves taking three steps u 1 = u n + ∆ tL ( u n ) Time u 2 = 3 4 u n + 1 4 u 1 + 1 4 ∆ tL ( u 1 ) u n +1 = 1 3 u n + 2 3 u 2 + 2 3 ∆ tL ( u 2 ) We start by rewriting the steps as: u 1 = u n + ∆ tL ( u n ) u 2 = 3 4 u n + 1 h u n + ∆ tL ( u n ) i 4 u n +1 = 1 3 u n + 2 h u 2 + ∆ tL ( u 2 ) i 3