IMA sola

Document Sample
IMA sola Powered By Docstoc
					                         Numerical investigation of
                         drop deformation in shear
                 Y. Renardy, S. Afkhami. Mathematics, Virginia Tech,
                     Jie Li, Cambridge. D. Khismatullin, Tulane U.

                                                                         equal visc

                              Influence of viscoelasticity?
                            K. Verhulst, P. Moldenaers, R. Cardinaels,

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

●   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
 mh drop hmatrix
 ●Density ratio rd / rm
 ●Capillary number Cahm’a/
 viscous force causing deformation / capillary force which keeps
 the drop together.                     Shear rate ’ (Utop-Ubottom plate )/
 ●Reynolds      number Re=rm’ a2 / hm drop tensiona
                                      Interfacial radius


Governing equations for a volume-of-fluid method for
Newtonian liquids
                   u          
       u  0,       u  u   p     S  F
                   t          
                              1   1 
                        
     S  u   u  , Fs      ns s

                               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.


                          Δyj                     Pijk.Cijk

                       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:

                     Ui-½        Ui+½

                                        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,  un1  0:

           ( )    u *                 Solve with multigrid method
                     t

   •Solve for   un+1:   un1  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

     • We know the solution at the nth timestep.

u u
 *       n
      (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          
        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 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 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
   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.

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





                         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.

VOF-Continuous Surface Force Formulation.
Velocity vector plots across centerline in x-z plane for different
mesh size at t=200Dt.

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∣
                                                               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
1/128     0.0000131      0.0000005          0.00000009
1/160     0.0000095      0.0000004          0.00000007
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
                        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

The next slides show the implementation of PROST for
                  viscoelastic liquids
Our governing equations use the Giesekus model:

                                               Model parameter
relaxation time                                (shear-thinning)


 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 Cahm’a/
                   =viscous force deforms drop / capillary force retracts drop.

 ●Viscosity      ratio mh     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

                 Interface cell
                -- Fluid properties are interpolated.
                -- This 'partly elastic' fluid changes properties
as the interface moves, and need not satisfy the stability

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

Our choice of implicit terms allows for the decoupling of variables
followed by inversion of tridiagonal matrices


 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.

System   Drop   Matrix
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




                       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

                                                                      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.

                                                              Domain 16R0x16R0x8R0


                                                              12days,16cpus,SGI Altix


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


Shared By: