Numerical investigation of drop deformation in shear Y. Renardy, S. Afkhami. Mathematics, Virginia Tech, Jie Li, Cambridge. D. Khismatullin, Tulane U. http://www.math.vt.edu/people/renardyy Re=15,Ca=0.2 equal visc Influence of viscoelasticity? K. Verhulst, P. Moldenaers, R. Cardinaels, KU-Leuven NSF, NCSA, Teragrid, VT-ARC,Michael Renardy,T.Chinyoka, M.A.Clarke. A useful way to recycle plastics is to melt and mix them to form incompatible polymer blends microstructure of HIPS. An incompatible polymer blend is a new material, combining the desired properties of the two original plastics. A cutaway diagram shows the geometry of the device in which a drop travels in the experiments of J. Waldmeyer (PhD 2008) Waldmeyer,Mackley,Renardy,Renardy 2008 ICR The cross-section on the right shows two sample stream lines of the base flow. If the drop is small, then we can set up the boundary conditions in our simple shear flow from the strain rates of the baseflow. The plan for this presentation ● We begin with a description of the Newtonian VOF-CSF method in the context of drop deformation. Nichols, Hirt, Hotchkiss, Los Alamos Sci. Lab. Rep. LA-8355, 1980. (SOLA-VOF) Kothe, Mjolsness, Torrey, Los Alamos Nat. Lab. Rep.,1991.(RIPPLE) Brackbill, Kothe, Zemach, J. Comput. Phys.1992(CSF) Zaleski et a lJ. C Phys. 1994 SURFER,(Institut d’Alembert) Li, Renardy,Renardy Phys.Fluids 2000 ● For some simulations, a higher-order method like VOF-PROST is needed. e.g. drop retraction when shearing stops. Renardy,Renardy, J C P (2002) VOF PROST Newtonian Khismatullin,Renardy,Renardy,JNNFM 2006 Giesekus Renardy,Renardy,Benyahia,Assighaou,PRE2009 ● Implementation of 3D Oldroyd B or Giesekus constitutive model. 2D: Hooper, de Almeida, Macosko, Derby, JNNFM, 2001 (finite element) Yue, Feng, Liu, Shen, JFM, 2005(phase field) 3D: Pillapakkam, Singh, J. Comp. Phys. 2001 (level set) Khismatullin,Renardy,Renardy JNNFM (2006) (3D) Aggarwal, Sarkar, JFM 2007. Front-tracking fin diff scheme. Hulsen et al… The dimensionless parameters for the Newtonian case are ●Viscosity ratio mh drop hmatrix ●Density ratio rd / rm ●Capillary number Cahm’a/ viscous force causing deformation / capillary force which keeps the drop together. Shear rate ’ (Utop-Ubottom plate )/ L ●Reynolds number Re=rm’ a2 / hm drop tensiona z Interfacial radius Initial z=Lz x=Lx y=Ly Governing equations for a volume-of-fluid method for Newtonian liquids u u 0, u u p S F t 1 1 1 S u u , Fs ns s 2 T R1 R2 C Fluid 1 ns C Fluid 2 •A color function C(x,y,z,t) is advected by the flowfield •We reconstruct the interface from the surface where C jumps •Body force includes the interfacial tension force •Periodic boundary conditions in the x and y directions, and no slip at the walls z=0, z=1 The spatial discretization is a regular Cartesian grid. Wijk C(i,j,k)=0 Δzk Δyj Pijk.Cijk Vijk C(i,j,k)=1 Δxi Uijk e.g. C(i,j,k)=1/3 Finite differences are used to calculate derivatives. Staggered grid. In a cell that is cut by the interface, fluid property is averaged. We began with SURFER, which treats high Re flows. (open source, S. Zaleski 1997) 1. The interface is reconstructed from C(i,j,k,t) as a plane in each 3D cell, or a line in a 2D cell: C=0 C=0.7 Ui-½ Ui+½ C=1 V1i,j,k V2i,j,k V3i,j,k In a cell that is cut by the interface, C = volume fraction of fluid 1. •PLIC – piecewise linear interface reconstruction •Lagrangian advection yields C at each timestep 2. VOF codes typically use the projection method on the momentum equations, and an explicit temporal discretization. Chorin 1967 • We know the solution at the nth timestep. Next, solve for an intermediate velocity field u* without p. u* un (un )un 1( (S) F)n t • Compute a correction p, using un 0, un1 0: p ( ) u * Solve with multigrid method t •Solve for un+1: un1 u* p t The explicit time integration scheme is not feasible for low Re drop-breakup simulations. The explicit scheme is stable when Dt << time scale of viscous diffusion: mesh2 / viscosity Implicit scheme would be slow because variables are coupled large full matrix Semi-implicit scheme with decoupling of u, v,w SURFER++ Li, Renardy, Renardy 1998 The Stokes operator is associated with viscous diffusion. • We know the solution at the nth timestep. Next, solve for an intermediate velocity field u* without p. u u * n 1 (u )u ( ( S ) F ) n n n t The Stokes operator causes the need for small t for low Re, so treat some of this expression implicitly * We developed and implemented a semi-implicit time integration scheme Li,Renardy,Renardy 1998 Let us take the X-component of the intermediate velocity field: u* un 1 n 1 * (u )u n n F (2 u ) t x x 1 * n 1 * ( u v ) ( u wn) y y x z z x The Stokes operator terms for u* are treated implicitly. t ) ( ( ) ( ( )u * =explicit terms I ( x (2 x y y z z We factorize this (Zang,Street,Koseff 1994). t ) I t ( ) I t ( )u * I x (2 x y y z z Finally, we invert tridiagonal matrices. Analogously for v* and w*. 3. How VOF-CSF-PLIC discretizes the interfacial tension force Fs . u u u p S F t Fs ns s , s |C|, ns C . |C | k = - div ns 0 in fluid 2, At the continuum level, C 1 in fluid 1. At the discrete level, C=volume fraction of fluid 1 per cell. Finite differences of the discontinuous function C give ns , more finite diffs & nonlinear combinations give k . The Continuous Surface Force Formulation (CSF) works well in an average sense for flows. • CSF Brackbill et.al 1992, Kothe,Williams,Puckett 1998,… Introduce a mollified C in Fs over a distance much larger than the mesh: ~ C(x ) C(x')(x' x,) dx' where f(x,e) is a kernel. Diffusion of surface tension force? • CSS Continuous Surface Stress Formulation Lafaurie et.al 1994, Zaleski,Li,Gueyffier,… T (I ns ns ) s , Fs T. Application: The Cambridge Shear System was used to obtain experimental data on drop deformation for oscillatory, step, & steady shear Wannaborworn, Mackley, Y. Renardy, 2002 PDMS drop in PIB = 4 mN/m h= 80 Pa.s for both Equal density VOF-CSF reproduces the initial transient for oscillatory shear experiment at 0.3Hz 250% strain 30.175 mm diameter drop. Numerical Expt Microchannel application: 3D Newtonian Stokes flow Ca=0.4, R0/H=0.34 VOF-CSF simulation (Re=0.1) Experimental data of Sibillo, Pasquariello, Simeone, Guido, Soc. of Rheol. Meeting, 2005. Inertia is destabilizing, so add a small amount to break the drop: H Re=2, Ca=0.4 Monodisperse droplets YRenardy, Rheol Acta 2007 VOF-CSF-PLIC does not converge with spatial refinement for capillary breakup of filament. Mesh View from top a/8 a/12 a/16 a/20 The simulation of surface-tension-dominated regimes produce small SPURIOUS CURRENTS Re=12, Ca=1.14Cac The simulation of a drop in another liquid with zero initial velocity. EXACT SOLUTION for all time: zero velocity. Prescribed surface Zero velocity boundary tension condition p Fsurface tension C, const. for sphere. p C const. Discretized in the same manner as l.h.s. SPURIOUS CURRENTS : VOF-Continuous Surface Force Formulation. Velocity vector plots across centerline in x-z plane for different mesh size at t=200Dt. Dx=a/20 Dx=a/12 Magnitudes of v increase in Lmax and L2, remain same in L1 Norms of v at 200th timestep Dt=10-5 should approach 0 as mesh-size decreases. Dx L = max∣ v∣ L1= ∣ v∣ dxdydz L2 = ∣ v∣ 2 dxdydz method 1/96 0.0017998 0.0000840 0.0000147 CSF 1/128 0.0018409 0.0000854 0.0000154 1/160 0.0018905 0.0000860 0.0000157 1/192 0.0019688 0.0000863 0.0000157 1/96 0.0037704 0.0001418 0.0000192 CSS 1/128 0.0038588 0.0001245 0.0000162 1/160 0.0036042 0.0001123 0.0000144 1/192 0.0039840 0.0001045 0.0000135 1/96 0.0000224 0.0000009 PROST 0.00000014 1/128 0.0000131 0.0000005 0.00000009 1/160 0.0000095 0.0000004 0.00000007 O(Dx)2 1/192 0.0000057 0.0000002 0.00000004 Moral of the story you can’t win the game by finite differencing C Our VOF-PROST algorithm implements 1. A sharp interface reconstruction and Lagrangian advection of the VOF function. JCP,Renardy,Renardy 2002 k n (x x 0 ) (x x 0 )A(x x 0 ) 0 Δzk x0 Δyj Optimization scheme yields the mean Δxi curvature = 2 tr(A) at cell centers. 2.A modified projection method with semi-implicit time integration is used for the momentum and constitutive equations. Li,Renardy,Renardy 1998 PROST achieves convergence of fragment volumes with mesh refinement Dx=a/12 t=22.5 Dx=a/16 t=24 Re=12 The next slides show the implementation of PROST for viscoelastic liquids Our governing equations use the Giesekus model: Model parameter relaxation time (shear-thinning) hp Initial condition for a drop in shear: zero viscoelastic stress and impulsive startup for velocities. Dimensionless parameters ●Reynolds number small Re=rm’ a2 / hm ● Density ratio 1= rd / rm ●Capillary number Cahm’a/ =viscous force deforms drop / capillary force retracts drop. ●Viscosity ratio mh drop hmatrix ● Weissenberg numbers We = ‘ l or Deborah numbers Dem= Wem(1-bm)/Ca Ded= Wed(1-bd)/(mCa) measures viscoelastic vs capillary effects. ● Retardation parameter b= hsolvent / htotal ● Giesekus model parameter 0<<0.5 Positive definite property of the extra stress tensor is retained in our algorithm. The time-dependent UCM eqns have an instability if an eigenvalue of the extra stress tensor is < -G(0). (Rutkevich 1967) This does not happen if the initial data are ‘physical’, but can happen numerically. Newtonian Viscoelastic Interface cell -- Fluid properties are interpolated. -- This 'partly elastic' fluid changes properties as the interface moves, and need not satisfy the stability constraint. We correct this numerical instability by adding a multiple of I to T over interface cells if eval < -G(0). We use our semi-implicit time integration scheme and operator factorization for the constitutive equation -λκ(T(n))2 Our choice of implicit terms allows for the decoupling of variables followed by inversion of tridiagonal matrices T(n+1) -λκ(T(n))2 VOF-PROST runs on shared-memory machines. Mesh Dx=Dy=Dz=a/12 typically use 16 cpus, SGI Altix, 10 days. 10 million unknowns per timestep, 200,000 timesteps 3D simulations for experimental results of Moldenaers et al are shown next. Viscosity System Drop Matrix ratio 1 VE NE 1.5 3 NE VE 1.5 4 NE NE 1.5 5 NE VE(BF2) 0.75 BF2 (Verhulst thesis) A Boger fluid drop in a Newtonian matrix. De_d = 1.54. Viscosity ratio 1.5, is more retracted with increased shear rate because the viscoelastic stresses at the drop tips act pull the drop in. - - - Newtonian CSF simulation ____ Oldroyd B CSF Ca=0.32 simulation o experiment Ca=0.14 Ca=0.14 Ca=0.32 Rotational flow inside the drop does not generate as much viscoelastic stress as when the fluids are reversed. A Newtonian drop in Boger fluid matrix. De_m = 1.89, viscosity ratio 1.5, b_m=0.68,Ca=0.35, increases the initial overshoot with increase in shear rate, and retracts over a long time scale. NE-NE steady state is here. o experiment (D decreases over longer time scale) Small def. theory : D is same as NE-NE. Greco JNNFM 2002 ____ Oldroyd B _._. Giesekus Shear-thinning… smaller stresses in the VE ‘ stress wake’ A new breakup scenario for a viscoelastic drop in a Newtonian matrix was found experimentally at visc ratio 1.5 Verhulst thesis 2008 Elongation and necking to t’=100 is followed by interfacial tension forming dumbbells joined by a filament. The viscoelastic filament thins but instead of breaking, pulls the ends to A second end-pinching elongates the coalescle. drop more than in the first. Beads form on the filament. Filament breaks. 3D Numerical simulations at a higher Ca=0.65, same Ded=0.92 , forms dumbbells and a uniform filament. The dumbbells end-pinch numerically when the filament is under-resolved. Dx=R0/12 Domain 16R0x16R0x8R0 Dt=.0005/g’ 12days,16cpus,SGI Altix Experimental Large stresses build up at the neck by t’=35 and grows on the filament side. Interfacial tension forms dumbbell shape. If the filament were constrained not to break then high stresses there would pull the ends together. 1D surface tension driven breakup of an Oldroyd- B filament in vacuum never breaks. (M.Renardy 1994,1995) Does the Boger fluid remain Oldroyd-B when the filament thins a long time? Effect of shear flow history A solution in Stokes flow does not depend on the initial condition. This uniqueness is lost when additional effects such as viscoelasticity are added (or for instance, inertia. IJMF 2008). NE-VE at visc ratio 0.75, Ca=0.5, Dem=1.54, breaks up VOF-PROST simulation with Giesekus parameter 0.1,mesh a/12. … but not when the shear rate is ramped up in small steps. Numerical simulations with the one-mode Giesekus model with model parameter 0.1. The same level of viscoelastic stress is associated with breakup or settling The same level of viscoelastic stress occurs with breakup or settling The End Questions?
Pages to are hidden for
"IMA sola"Please download to view full document