Dynamics of thin liquid films in zero gravity: stabilization and destabilization

Aneet DHARMAVARAM NARENDRANATH, Ph.D.
Lecturer, Mechanical Engineering-Engineering Mechanics, 
 Michigan Technological University,
 Houghton, MI, USA

Michigan Technological University

Houghton, MI 49931

Pattern Selection in evaporating liquid films: Short wavelength Instabilities

SystemOpen["/home/dnaneet/Research_4.0/WTC/2016/Movie2_ISTP_c.avi"]

Long Wave film dynamics and who is interested
(Subscript[h, t])+(Subscript[(f(h)(Subscript[h, x]) + g(h)(Subscript[h, xxx])), x])=0

L1 = 185.62; TMax = 100000; km = 0.0381; 
hSol1 = h /. NDSolve[{Complete[50, 1/3, 1, 6.5, 0., 0.25], 
       h[0, y, t] == h[L1, y, t], h[x, 0, t] == h[x, L1, t], 
       h[x, y, 0] == 1 + (-0.05*Cos[2*Pi*(x/L1)] - 0.05*Sin[2*Pi*(x/L1)])*
          Cos[2*Pi*(y/L1)]}, h, {x, 0, L1}, {y, 0, L1}, {t, 0, TMax}][[1]]; 
{TMin1, TRup1} = InterpolatingFunctionDomain[hSol1][[3]];

Animate[Plot3D[hSol1[x,y,facTRup1],{x,0,L1},{y,0,L1},PlotRange->{{0,L1},{0,L1},{0.8,3.}},AxesLabel->{"x","y","h"},PlotLabel->"Film evolution",PlotPoints->25,PlotStyle->Thick,Mesh->Full,ImageSize->Large],{fac,0,0.95,0.001},AnimationRate->0.1,AnimationRepetitions->2]

Pattern Selection

FData1=Abs[Fourier[Table[hSol1[x,y,0.9TRup1],{x,0,L1,L1/50},{y,0,L1,L1/50}]]];
FData1[[1,1]]=0;
FData=Threshold[FData1/Norm[FData1/.FData1[[1,1]]->0,"Frobenius"],{"Hard",0.1Max[FData1/Norm[FData1,"Frobenius"]]}];
MatrixPlot[FData[[1;;10,1;;10]],ImageSize->Medium]

Pendant drop: Effect of destabilizing gravity
Imagine liquid film on ceiling [isothermal non-porous ceiling],  (Subscript[h, T]) + (Subscript[((h^3)(Subscript[h, xxx])), x]) + (Subscript[((h^3)(Subscript[h, x])), x])

Clear[x,t,xMin, xMax,TMax,frac,m,time];
xMin=0; xMax=2.Pi(Sqrt[2]/1);TMax =20;
   hsol=u /. NDSolve[{D[u[t, x], t] == -D[u[t, x]^3*D[u[t, x], x, x, x], x] - 
        D[u[t, x]^3*D[u[t, x], x], x], u[0, x] == 1 - 0.1*Sin[1*(x/Sqrt[2])], 
    pbcs}, u, {t, 0, TMax}, {x, xMin, xMax},met][[1]];

Plot[hsol[TMax,x],{x,xMin,xMax},AxesLabel->{"Substrate","Film Profile"},Epilog->Inset[Framed[Style["Pendant Drop",12]],{7,2}]];
Animate[Plot[hsol[time TMax,x],{x,xMin,xMax},AxesLabel->{"Substrate","Film Profile"},PlotRange->{{xMin,xMax},{0,3.}}],{time,0,1,0.005},AnimationRate->.12,AnimationRepetitions->2]

Thermocapillary instabilities: The base case
(Subscript[h, T]) + (Subscript[S((h^3)(Subscript[h, xxx])), x]) -(Ga/3)(Subscript[((h^3)(Subscript[h, x])), x]) + (Subscript[(M/Pr) [((h^2)/((1 + Bi)))(Subscript[h, x])], x])= 0

Surface tension and gravity are stabilizing mechanisms while thermocapillary destabilizes the liquid film

Thermocapillary destabilization results in “finger structures” (primary and secondary).

Rupture is defined as that point in time when the problem is overly stiff to solve.

Clear[x,t,xMin, xMax,TMax,k];
xMin=-Pi/0.0677; xMax= Pi/0.0677;TMax=2000.;k=0.067;
Quiet[hbase= u /. NDSolve[{D[u[t, x], t] == -100*D[u[t, x]^3*D[u[t, x], x, x, x], x] + 
        (1/3)*D[u[t, x]^3*D[u[t, x], x], x] - 
        5*D[(u[t, x]/(1 + u[t, x]))^2*D[u[t, x], x], x], 
      u[0, x] == 1 - 0.1*Cos[k*x],pbcs}, u, {t, 0, TMax}, 
     {x, xMin, xMax},met][[1]];]

Animate[Plot[hbase[time TMax,x],{x,xMin,xMax},AxesLabel->{"Substrate","Film Profile"},PlotRange->{{xMin,xMax},{0,3.}},GridLines->Automatic],{time,0,1,0.02},AnimationRate->0.1,AnimationRepetitions->2]

Zero gravity: Surface tension and thermocapillarity
(Subscript[h, T]) + (Subscript[S((h^3)(Subscript[h, xxx])), x])+ (Subscript[(M/Pr) [((h^2)/((1 + Bi)))(Subscript[h, x])], x])= 0

Gravity is turned “OFF” by ignoring the relevant nonlinear term from the evolution equation

Clear[x, t, xMin, xMax, TMax, k]; 
xMin = -(Pi/0.0677); xMax = Pi/0.0677; TMax = 2000.; k = 0.067; 
  hbase0 = u /. NDSolve[{D[u[t, x], t] == 
         -100*D[u[t, x]^3*D[u[t, x], x, x, x], x] + 
          (0/3)*D[u[t, x]^3*D[u[t, x], x], x] - 
          5*D[(u[t, x]/(1 + u[t, x]))^2*D[u[t, x], x], x], 
        u[0, x] == 1 - 0.1*Cos[k*x], pbcs}, u, {t, 0, TMax}, 
       {x, xMin, xMax}, met][[1]];

Animate[Plot[hbase0[time TMax,x],{x,xMin,xMax},AxesLabel->{"Substrate","Film Profile"},PlotRange->{{xMin,xMax},{0,3.}},PlotStyle->{Thick,Red},GridLines->Automatic],{time,0,1,0.02},AnimationRate->.1,AnimationRepetitions->2]

Comparison: Base case vs Zero gravity

TMax=2000.;xMin=-Pi/0.067;xMax=Pi/0.067;
Animate[
Plot[{hbase[time TMax,x],hbase0[time TMax,x]},{x,xMin,xMax},PlotRange->{{xMin,xMax},{0,3}},PlotStyle->{{Thick,Black},{Thick,Red}},PlotLegends->{"Base","Zero Gravity"},ImageSize->Medium],
{time,0,1,0.05},AnimationRate->.07,AnimationRepetitions->2]

Without gravity, film destabilization proceeds quicker and film rupture occurs at an earlier time.

It is also observed that primary and secondary finger structures grow to larger amplitudes without gravity.

Base case with porous substrate
(Subscript[h, T]) + (Subscript[S((h^3)(Subscript[h, xxx])), x]) -(Ga/3)(Subscript[((h^3)(Subscript[h, x])), x]) + (Subscript[(M/Pr) [((h^2)/((1 + Bi)))(Subscript[h, x])], x])= S Pm((Subscript[h, xx]) - Ga h)

Clear[x,t,xMin, xMax,TMax,k,Pm,mult,Ga];
xMin=-Pi/0.0677; xMax= Pi/0.0677;TMax=2000.;k=0.067;
Pm=0.0025;Ga=1.; mult=1.;
hbasep= u /. NDSolve[{D[u[t, x], t] == -100*D[u[t, x]^3*D[u[t, x], x, x, x], x] + 
        (1./3.)*D[u[t, x]^3*D[u[t, x], x], x] - 
        5*D[(u[t, x]/(1 + u[t, x]))^2*D[u[t, x], x], x]+ mult*Pm*(D[u[t, x], x, x] - Ga*u[t, x]), 
      u[0, x] == 1 - 0.1*Cos[k*x],pbcs}, u, {t, 0, TMax}, 
     {x, xMin, xMax}, met][[1]];

Animate[Plot[hbasep[time TMax,x],{x,xMin,xMax},AxesLabel->{"Substrate","Film Profile"},PlotRange->{{xMin,xMax},{0,1.4}},PlotStyle->{Thick,Black},GridLines->Automatic,ImageSize->Large],{time,0,1,0.005},AnimationRate->0.04,AnimationRepetitions->2]

When an isothermal porous substrate is introduced, it drains the liquid film.

Effect of porous substrate on film evolving in zero gravity
(Subscript[h, T]) + (Subscript[S((h^3)(Subscript[h, xxx])), x])+ (Subscript[(M/Pr) [((h^2)/((1 + Bi)))(Subscript[h, x])], x])= S Pm((Subscript[h, xx]) - Ga h)

Clear[x,t,xMin, xMax,TMax,k,Pm,mult,Ga];
xMin=-Pi/0.0677; xMax= Pi/0.0677;TMax=8000.;k=0.067;
Pm=0.15;Ga=0.; mult=1.;
hbasep0= u /. NDSolve[{D[u[t, x], t] == -100*D[u[t, x]^3*D[u[t, x], x, x, x], x] +
        (Ga/3.)*D[u[t, x]^3*D[u[t, x], x], x] - 
        5*D[(u[t, x]/(1 + u[t, x]))^2*D[u[t, x], x], x]+ mult*Pm*(D[u[t, x], x, x] - Ga*u[t, x]), 
      u[0, x] == 1 - 0.1*Cos[k*x],pbcs}, u, {t, 0, TMax}, 
     {x, xMin, xMax},met][[1]];

Animate[Plot[hbasep0[time TMax,x],{x,xMin,xMax},AxesLabel->{"Substrate","Film Profile"},PlotRange->{{xMin,xMax},{0,3.}},PlotStyle->{Thick,Red},GridLines->Automatic],{time,0,1,0.02},AnimationRate->0.07,AnimationRepetitions->2]

The porous substrate leads to damping of thermocapillary fingers

The film also avoids rupture for a long duration ((Subscript[T, rup])≈4 x in case below) of time.  This time of rupture is directly related to porosity number

Multiple events: Initial evolution similar to base case and zero gravity case on non-porous substrate, intermediate evolution dominated by decceleration of film dynamics and damping of thermocapillary structures, long time film dynamics show arrestation of film dynamics as compared to base case or that with zero gravity.

Comparison of film dynamics of all cases (base, zero gravity, zero gravity with porous substrate)

xMin=-Pi/0.0677; xMax= Pi/0.0677;TMax=2000.;TMaxp=12000.;
Animate[
Plot[{hbase[fraction TMax,x],hbase0[fraction TMax,x],hbasep0[fraction TMax,x]},{x,xMin,xMax},PlotStyle->{{Thick,Black},{Thick,Red},{Thick,RGBColor[0.49,0.52,0.36]}},PlotLegends->{"Base case, non porous substrate","Zero Gravity, non porous substrate","Zero Gravity, porous substrate"},PlotRange->{{xMin,xMax},{0,3.2}},ImageSize->Medium]
,{fraction,0,1,0.02},AnimationRepetitions->3]

All three situations have the same initial conditions and periodic boundary conditions

The case with zero gravity evolution of thermocapillary instabilities on a non porous substrate is the first to rupture at (Subscript[T, rup])=2000

The case with zero gravity evolution of thermocapillary instabilities on a porous substrate has not yet ruptured at (Subscript[T, rup])=2000.

Initially, the case film on a porous substrate has similar dynamics to the base case (surface tension, gravity, thermocapillarity) and leads leads the latter.

However, the film on a porous substrate experiences a decceleration of film dynamics along with a damping of thermocapillary structures.

Exploring film dynamics with “Recurrence Plots”
What is a recurrence plot? Θ(x)(ϵ-(Abs[(Subscript[x, i])-(Subscript[x, j])]))

The differential equation: y'(x)==y(x) cos(y(x)+x), with boundary condition y(0)=1 is numerically solved.  The solution shows oscillation and evolution of periodic structures.

The power spectrum of the DFT of y is calculated and compared with the recurrence plot of the same.

Clear[ysol,dft,y,x];
ysol=NDSolveValue[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,200}];
dft=Abs[Fourier[Table[ysol[x],{x,0,200,1}]]]^2;
p1=Plot[ysol[x],{x,0,200},ImageSize->250{1,1},BaseStyle->{FontFamily->"SansSerif",FontSize->12},AxesLabel->{"x","y"},PlotRange->{All,{0,1.1}}];
p2=ListPlot[dft[[1;;100]],Filling->Axis,PlotRange->{All,{0,0.9}},Frame->True,FrameLabel->{"","","","PowerSpectrum"},ImageSize->Medium];
p3=MatrixPlot[UnitStep[0.00005-DistanceMatrix[ysol[Range[0,200,.2]]]],Frame->True,FrameLabel->{"x","x","","Recurrence Plot"},ImageSize->Medium];
p4=MatrixPlot[UnitStep[0.00005-DistanceMatrix[ysol[Range[100,200,.2]]]],Frame->True,FrameLabel->{"x","x","","Recurrence Plot \n(Long time evolution)"},ImageSize->Medium];

GraphicsGrid[{{p1,p2},{p3,p4}}]

Note on Recurrence Plots

The recurrence plot shows a tessellated structure that is an aspect of spatial (or temporal) periodic structures. Recurrence plots allow for qualitative and quantitative visualization of spatiotemporal evolution of dynamical systems.

A recurrence plot is a matrix plot of a square matrix that is calculated by the equation (θ(x) is the Heaviside function, ϵ is neighbourhood or cut-off distance within which recurrence is sought, (Subscript[x, i]) and (Subscript[x, j]) are states of the dynamical system, typically, at time=i > time=j.  In case of film dynamics, these are, for the current demonstration, spatial locations.  ||...|| is some sort of Norm.):

Film dynamics with zero gravity and porous substrate: RP
damping of thermocapillary structures

This simulation runs for 4(Subscript[xT, rup])(≈8000).

Animate[
Plot[hbasep0[np 8000,x],{x,xMin,xMax},AxesLabel->{"Substrate","Film Profile"},PlotRange->{{xMin,xMax},{0,3.}},PlotStyle->{Thick,Red},GridLines->Automatic,ImageSize->Medium,BaseStyle->{FontFamily->"SansSerif",Bold,FontSize->12}]

With[{stepSize=.2,end=xMax,fraction=np},
MatrixPlot[UnitStep[0.0005-DistanceMatrix[hbasep0[fraction 8000,Range[xMin,end,stepSize]]]],ColorFunction->"Monochrome",Frame->True,FrameLabel->{"","Base + 0g + porosity \n (Narendranath et al. 2016)" },ImageSize->Medium]
],{np,0,1,0.02},AnimationRate->.04,AnimationRepetitions->2]

Comparison of film dynamics with Recurrence Plots

Recurrence plots are now used to study film spatiotemporal dynamics

Provides insight on differing times at which periodic structures arise or dissipate.

Zero-g and Base case have a time lag between them.

Zero-g when compared to zero-g+porous substrate show that the floral pattern, indicative of multiple thermocapillary fingers forms in both at around 0.45 (Subscript[T, rup]). However, these thermocapillary fingers are held in “suspension” at (Subscript[T, rup]) with a porous substrate

xMin=-Pi/0.0677; xMax= Pi/0.0677;
Animate[g11=With[{stepSize=.25,end=xMax,fraction=np},
MatrixPlot[UnitStep[0.0005-DistanceMatrix[hbase[ fraction TMax,Range[xMin,end,stepSize]]]],ColorFunction->"Monochrome",Frame->True,FrameLabel->{"","Base \n(Krishnamoorthy et al. 1995)"},ImageSize->300{1,1}]
];
g12=With[{stepSize=0.25,end=xMax,fraction=np},
MatrixPlot[UnitStep[0.0005-DistanceMatrix[hbase0[fraction TMax,Range[xMin,end,stepSize]]]],ColorFunction->"Monochrome",Frame->True,FrameLabel->{"","Base + 0g \n(Narendranath et al. 2015)"},ImageSize->300{1,1}]];
g13=With[{stepSize=0.25,end=xMax,fraction=np},
MatrixPlot[UnitStep[0.0005-DistanceMatrix[hbasep0[fraction TMax,Range[xMin,end,stepSize]]]],ColorFunction->"Monochrome",Frame->True,FrameLabel->{"","Base + 0g + porosity \n (Narendranath et al. 2016)" },ImageSize->300{1,1}]
];
GraphicsGrid[{{g11,g12,g13}(*,{g21,g22,g23}*)}],{np,0,1,0.02},AnimationRate->0.05,AnimationRepetitions->2]

Linear Stability

Linear stability analysis is applied to find an expression for growth rate of instabilities and fastest growing wavelength (Similar to Davis et al., Burelbach et al.)

The growth rate for LW instabilities is first depicted along with the effect of magnitude of various mechanisms.

Clear[f,Ga,Sn,Bi,M,q,Pm,qmaxval];
f=-((Ga (q^2))/3)+((M (q^2))/(((1+Bi))^2))+Pm (-Ga-(q^2))-(q^4) Sn;
Bi=0.1;
qmaxval=2;
Manipulate[Plot[f/.{Ga->xg 1.,Sn->xs 100.,Pm->xp 1/2,M->xm 5.},{q,0,qmaxval},PlotRange->{{0,qmaxval},{0,8}},BaseStyle->{FontSize->18},PlotStyle->{Thick,Black},ImageSize->Large,GridLines->Automatic,Frame->True,FrameLabel->{"wavenumber \n(inverse wavelength)","Growth rate"}],{{xg,0.,"Gravity"},0,10},{{xs,.6,"Surf Ten"},0.1,1},{{xp,0.,"Porosity"},0,10},{{xm,5,"ThermoCap."},0,5}]

Conclusions

Realization that the multitude of nonlinear dynamics are coupled and control of- depends on most destructive effect and its cascade

Control strategies for porous substrate: multiple substrates of varying porosity?  Multi-dimensional effects of porous substrate?

Effect of evaporation? (Turn ON a term)

Analytical methods for solution: Soliton theory, perturbation method, differential transform method

Acknowledgement

Dr.Jeffrey S. Allen from Department of Mechanical engineering-engineering mechanics, Michigan Technological University

Dr.Robert W. Kolkka and Dr.Allan A. Struthers,Department of Mathematical sciences, Michigan Technological University

Dr. James C. Hermanson, Department of aeronautics and astronauticcs, University of Washington at Seattle

References

Burelbach, James P., Seymour G.Bankoff, and Stephen H.Davis."Nonlinear stability of evaporating/condensing liquid films." Journal of Fluid Mechanics 195 (1988) : 463 - 494.

Oron, Alexander. “Three-dimensional nonlinear dynamics of thin liquid films.” Physical review letters 85, no. 10 (2000): 2108.

Davis, S. H., and L. M. Hocking. “Spreading and imbibition of viscous liquid on a porous base.” Physics of Fluids (1994-present) 11, no. 1 (1999): 48-57.

Narendranath, Aneet D., James C. Hermanson, Robert W. Kolkka, Allan A. Struthers, and Jeffrey S. Allen. “The Effect of Gravity on the Stability of an Evaporating Liquid Film.” Microgravity Science and Technology 26, no. 3 (2014): 189-199.

Narendranath, Aneet, James Hermanson, Allan Struthers, Robert Kolkka, and Jeffrey Allen. “Thermocapillarity driven Instabilities in thin liquid layers subject to long-wave analysis.” In APS Division of Fluid Dynamics Meeting Abstracts, vol. 1, p. 9004. 2012.

Allan A. Struthers, Robert W. Kolkka, Jeffrey S. Allen., and Narendranath, Aneet D.,“Accurately Resolving Rupture in Fourth Order Partial Differential Equations Modeling Thin Films Using NDSolve”, Wolfram Technology Conference. 2010

Spatiotemporal Dynamics as depicted by a traditional Phase portrait

Clear[x,t,xMin, xMax,TMax,k];
xMin=-Pi/0.0677; xMax= Pi/0.0677;TMax=2000.;k=0.067;
s= NDSolve[{D[u[t, x], t] == -100*D[u[t, x]^3*D[u[t, x], x, x, x], x] + 
        (1/3)*D[u[t, x]^3*D[u[t, x], x], x] - 
        5*D[(u[t, x]/(1 + u[t, x]))^2*D[u[t, x], x], x], 
      u[0, x] == 1 - 0.1*Cos[k*x],pbcs}, u, {t, 0, TMax}, 
     {x, xMin, xMax},met];

ParametricPlot[Evaluate[{u[0.8TMax,ξ],D[u[0.8TMax,ξ],ξ]}/.s/.ξ->x],{x,xMin,xMax},AspectRatio->1/2]

Initialization of solution parameters and plot properties

Needs["HierarchicalClustering`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
SetOptions[Plot,ImageSize->Large,PlotStyle->{Black,Thick},BaseStyle->{FontFamily->"SansSerif",Bold,FontSize->24}];
SetOptions[MatrixPlot,ImageSize->Small,BaseStyle->{FontFamily->"SansSerif",Bold}];
SetOptions[ParametricPlot,ImageSize->Medium,BaseStyle->{FontFamily->"SansSerif",Bold,FontSize->19}];
SetOptions[ListPlot,ImageSize->Medium,BaseStyle->{FontFamily->"SansSerif",Bold},PlotStyle->{Thick,Black},PlotMarkers->None,Filling->Axis];
SetOptions[Plot3D,PlotPoints->20,AxesLabel->{"x","y","h"},ImageSize->Medium,BaseStyle->{FontSize->14}]

$Pre=(SetOptions[EvaluationNotebook[],ShowPredictiveInterface->True];#)&;
$Post=Function[expr,If[FreeQ[expr,SparseArray],expr,SetOptions[EvaluationNotebook[],ShowPredictiveInterface->False];expr],HoldAll];

(*Switch off some messages*)
Off[NDSolve::mxsst]; 
Off[NDSolve::ibcinc]; 
Off[NDSolve::eerr]; 
Off[NDSolve::ndsz]; 


pbcs={ Derivative[0, 1][u][t, xMin] == Derivative[0, 1][u][t, xMax], 
      Derivative[0, 2][u][t, xMin] == Derivative[0, 2][u][t, xMax], 
      Derivative[0, 3][u][t, xMin] == Derivative[0, 3][u][t, xMax], 
      u[t, xMin] == u[t, xMax]};
met={ MaxSteps -> 100000, Method -> {"MethodOfLines", 
       "Method" -> "LSODA", "TemporalVariable" -> t, "SpatialDiscretization" -> 
        {"TensorProductGrid", "MinPoints" -> 800, "MaxPoints" -> 1200, 
         "DifferenceOrder" -> 5}}};



Clear[Eq5,Complete,h,S,G,K1,M,mult,Pm]
Eq5[h_,{S_,G_,K1_,M_,mult_,Pm_}]:=(∂_t)h+Div[-(h^3) G Grad[h,{x,y}]+(h^3) S Grad[Laplacian[h,{x,y}],{x,y}]+M (h/(1+h))^2 Grad[h,{x,y}],{x,y}]+-(mult Pm) (Laplacian[h,{x,y}]-G h)==0;
SetCoordinates[Cartesian[x,y,z]];
Complete[S_,G_,K1_,M_,mult_,Pm_]:=Eq5[h[x,y,t],{S,G,K1,M,mult,Pm}];
TraditionalForm[Complete[S,G,K1,M,mult,Pm]];