1. We are almost done with the complete code. However, a few important steps remain. DNS of Multiphase Flows — Simple Front Tracking Direct Numerical Simulations of Multiphase Flows-6 Surface Tension, Unequal Viscosities, and Higher Order Time Integration Gretar Tryggvason 2. In this lecture we will do three things. We will add a way to include surface tension, we will modify the DNS of Multiphase Flows — Simple Front Tracking code so that the viscosities in the different fluids can be different, and we will increase the order of the time integration. In this lecture we will complete our code by: • Adding surface tension • Allowing the viscosities in the different fluids to be different • Make the time integration higher order The resulting code is a fully functional one, and allows us to simulate simple multi fluid problems 3. First the surface tension. DNS of Multiphase Flows Surface Tension
4-1. The key thing about surface tension is that it is a tension. The force is parallel to the interface and if DNS of Multiphase Flows the interface is flat, the pull on one end of a small interface segment is balanced by the pull on the other The surface tension at a point is given by side. If, however, the interface is curved, there is a net force in the normal direction that is proportional t e f σ = σκ n to the curvature. For our purpose it is important that we need the force on the part of an interface that We need, however, the total force on − t s lies within a control volume---not the force per unit area. While we may be tempted to write down the n l a small segment of the front. For force per unit area and then multiply it by the area of the interface segment, it is better to go back to the two-dimensional flows we can use the definition of the curvature s is the coordinate definition of the surface tension. We can get there from the standard definition of the interface force by along the interface κ n = ∂ t using that the curvature times the normal is, by definition, equal to the derivative of the tangent vector ∂ s The force on a small element l is with respect to the interface length s. The integral over the curvature times the normal over the interface Z Z Z ∂ t δ f σ f σ ds = ∂ sds = σ � t e − t s � segment delta s is, if sigma is constant, exactly equal to the difference in the pull on either edge times l = σκ n ds = σ ∆ s l ∆ s l ∆ s l the surface tension coefficient, or the surface tension coefficient times the difference between the Z Thus, we do not need to find the curvature, just the tangents at the edges. tangent vectors, which is generally much simpler 4-2. Thus, we only need to find the tangents to the interface and not the curvature, which simplifies the DNS of Multiphase Flows computations considerably. Furthermore, we end up with a more robust algorithm that is conservative in The surface tension at a point is given by the sense that the force on the end of one segment is exactly equal and opposite to the force on the t e f σ = σκ n segment attached to it. This is important since for closed surfaces the net surface force should be zero We need, however, the total force on − t s and when we sum over all the interface elements, the forces on adjacent elements cancel. n l a small segment of the front. For two-dimensional flows we can use 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 = ∂ sds = σ � t e − t s � l = σκ n ds = σ ∆ 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 5. The numerical implementation is straightforward. We compute the surface force on an interface DNS of Multiphase Flows segment that consists of half of the elements attached to an interface point. The tangent force that pulls � Here we take the on the segment to the right of point l, halfway between point l and l+1, is given by the difference in the t l +1 / 2 − t l − 1 / 2 interface segment coordinates of these points divided by delta s, the distance between them. Delta s is computed as the around each interface � � − x l x l +1 point to consist of half x l − 1 square root of the difference in the x coordinates squared plus the difference in the y coordinates the distance to the points on either side squared. Thus, the force on the interface segment containing point l, and consisting of half the element The tangents at the end of the segment are found using a � before and after the point is given by delta f equal to sigma times the tangent vector at l+1/2 minus the centered approximation tangent vector at l-1/2. � � / ∆ s t l +1 / 2 = x l +1 − x l ( x l +1 � x l ) 2 + ( y l +1 � y l ) 2 where p ∆ s = The surface force acting on the interface segment around point l is: δ f σ � � l = σ t l +1 / 2 − t l − 1 / 2
6. The main consideration when transferring the force from the interface to the fixed fluid grid is that the DNS of Multiphase Flows total force must be conserved. Thus, the surface integral over the force per unit area (or length in two- The surface tension is distributed to the fixed grid. On the dimensions) must be equal to the volume integral over the force per unit volume (or area in two front the force is per unit area and on the grid the force is per dimension) on the fixed grid. The total force on the interface segment that will be transferred to a given unit volume. The total force is conserved, so that: Z Z grid point is the integral over that segment. Each segment generally consists of several interface f σ f σ s ds = v dv ∆ s l δ V Z Z elements and we must sum over those elements. For each element we found the force, denoted by delta Z The total force on an interface segment is f_l, on the last slide by subtracting the tangents. On the fixed grid we need the force per unite area in Z X F σ f σ δ f σ l = s ds ≈ l ∆ s l l two-dimensions, and the total force is the integral over the control volume surrounding each grid point, where the integral is over the part of the interface Z or f subscript i,j times delta x and delta y. contributing to a given grid point and the sum is over the interface segment that do. X The total force at a given grid point is Z F σ f σ v dv ≈ f σ i,j = i,j ∆ x ∆ y δ V 7. The total force on a surface segment that is transferred to a given grid point must be equal to the DNS of Multiphase Flows total force on the grid. The force transferred from each element is already the integral over the element, Using that the force at each grid point, F σ i,j = F σ so we just need to sum those over all elements that are close to each grid points, and to recover the l transferred from the front, is: force per unit volume on the fixed grid we divide by delta x time delta y. In general each surface element l w l δ f σ Sum over elements i,j X can contribute to more than one point on the fixed grid and we divide the surface force between the f σ i,j = ∆ x ∆ y “close” to grid point i,j l different grid points based on weights w, where the superscript l denotes the element and the subscript l w l f σ The weights are generally taken to be the same as i,j denotes the fixed grid points. Obviously, the weights for each element must add up to one, so that the i,j those used to interpolate the velocities to the front total force is accounted for. When we use a staggered grid, where the different velocity components On a staggered grid, the x -component of the surface force is each have their own control volume, the different component of the surface force are distributed to distributed to the u -nodes and the y -component to the v - different locations on the fixed grid. Once we have the surface force per unite volume on the grid, those nodes. can be added to the discrete Navier-Stokes equations. The surface forces on the fixed grid are then added to the discrete Navier-Stokes equations 8. For the fluid flow we can assume that a flat boundary can be represented as a symmetry boundary. DNS of Multiphase Flows Although this is mostly used for inviscid flows, we can also use it for no-slip boundaries if we add the At solid walls we prescribe symmetry boundary conditions so constraint of a zero tangent velocity. We assume that the same is true if we add surface tension so the the net force at the wall would be tangent to the wall surface tension for the boundary nodes can be set by assuming a symmetric, or a ghost, interface on the For weights based on bilinear interpolation, the force is other side of the boundary. This results in a zero normal force by symmetry and for bilinear interpolation, distributed to the four nearest grid points so no adjustment is needed for the normal velocity. For the tangential component or area weighting, no adjustment is needed for the normal velocity since we are not solving for the we need to add the contribution from the “ghost” part velocity at the boundary. The tangent force, at a point half a grid spacing away from the boundary, has to For the bottom boundary be modified by adding the contribution from the ghost interface on the other side. j = 2 ( f x ) i,j =2 = ( f x ) i,j =2 + ( f x ) i,j =1 Similar adjustments are done j = 1 for the other boundaries
Recommend
More recommend
