CC BY
26
RESPONSE OF A STRATIFIED VISCOUS HALF-SPACE TO A PERTURBATION OF THE FREE SURFACE
*S.A. Chivilikhin, 2A.S. Amosov, F. Melikhov
1 St. Petersburg National Research University of Information Technologies, Mechanics and Optics, St. Petersburg, Russia 2 Corning Incorporated, Corning, USA chivilikhin@gmail.com, amosovas@corning.com, ivan.melikhov@gmail.com
PACS PACS 47.15.G-
The flow of a highly viscous liquid in a half-space due to the deformation of the free surface is investigated. The viscosity of the layer adjoining to the free surface is different from the viscosity of the remaining half-space. In the framework of small perturbation theory, the relationship between the deformation of the free surface and the deformation of the layer/half-space interface is obtained. It was demonstrated that the volume and geometrical center of the perturbation on the interface and on the free surface are the same. The dependence of the perturbation's amplitude and width on layer thickness was investigated. The results of numerical and analytical calculations are close, even for moderate free surface perturbations. Keywords: quasisteady-state Stokes approximation, free surface, linear approximation, evolution of perturbation.
1. Introduction
In this paper, we consider an infinite viscous half-space covered by a fluid layer with a different viscosity. Using a Fourier transform-based approach, we calculated the deformation of the layer/half-space interface due to the perturbation of the free surface.
A Newtonian viscous incompressible fluid was considered. Viscosities of the layer and the half-space may not necessarily be the same, but they are assumed to be constant within each domain. We considered the localized perturbation of the fluid near the free surface, assuming it decays to zero infinitely far from the origin. Two formulations were developed: i) 2D axisymmetric to describe flow perturbation caused by a small particle (fig. 1-a) and ii) 2D planar to represent perturbation due to a long groove (fig. 1-b). In both cases, the shape of the initial perturbation of the free-surface is described by a smooth analytical function.
In section 2, we analytically derive the relation between the shape and amplitude of the free surface perturbation and the corresponding deformation of the layer/half-space interface.
In section 3, the analytical solution is compared to the numerical solution obtained with COMSOL Multiphysics finite element package.
2. Analytical model
The analytical model was developed under the following assumptions:
• Planar or axisymmetric viscous flow
• The fluid dynamics can be described by the Stokes equations ( [6], [1], [7], [3], [4],
[9], [2])
• Perturbation amplitude is small compared to the layer thickness ( [5], [8], [10], [11])
• Fluid viscosity is constant within each domain (layer and half-space)
2.1. The problem statement
In the analytical model, a fluid layer of the thickness H covering an infinite halfspace is considered. The layer has viscosity ¡i and the half-space has viscosity ¡i\. Consider the Cartesian coordinate system as shown on the figure 2 below. The axes x\ and x2 form a non-deformed planar boundary and the axis z is oriented through thickness having origin on the free surface of the layer.
/
(tmm
=1= / ^
Fig. 2. Coordinate system used in the analysis
The Stokes equations for the surface layer take the form:
#LVZ + dl,Vz = I dz P
d^zVa + df^Va = - daP (1)
1
dz Vz + df Vf = 0,
where Va, Vz are the velocity components, P is the pressure, ¡i is the viscosity of the fluid in the layer. The summation over the repeated indexes is assumed. From (1), it can be seen that pressure is a harmonic function:
d2zP + dffP = 0. (2)
Zero tangential stress and a given normal velocity are assumed on the external free boundary. Corresponding boundary conditions are:
Vz \z=0 = Vz0, l(da Vz + dz Va^ |z=0 = 0. (3)
Similarly to (1), the Stokes equations for the half-space domain become:
CVz + dffVz = - dz P ¡1
82zzVa + d^fVa = - daP (4)
¡¡1
dz Vz + df Vf = 0 32ZZP + dffP = 0. Zero perturbation infinitely far from the free surface is required:
= 0. (5)
The velocity field must be continuous through the interface:
Vz lz=H-0 = Vz |z=H+0 , Va\z=H-0 = Va\z=H+0 . (6)
Also, both the normal and tangential stresses must be continuous:
(-P + 2idzVz)\z=H-0 =(-P + 2lidzVz)\z=H+0 ,
¡(daVz + dzVa)\z=H-0 = Il(daVz + dzVa)\z=H+0 .
The equations of the evolution of the free boundary and the interface are the following:
dth0 = Vz \z=0 , dth = Vz\z=h , (8)
where h0(xl,x2,t) and h(xl,x2,t) are the ¿-displacements of the free boundary and interface respectively.
2.2. Calculation of the perturbation decay
Applying the Fourier transform, defined as (9) over the xl and x coordinates to the equations (1)-(8):
fk = i f (x)e-ikx d2x (9)
Response of a stratified viscous half-space one can get the following equations:
d2V7,k- k2V7
-dz Pk ¡
d^zVak - k2Vak = ^Pk
zz ¡
dz Vzk + ikf Vfk = 0
d2Pk - k2Pk = 0
zk
0 < z < H,
(10)
Vzk\z=0 = Vzo k, l(ikaVzk + dz Vak)\z=0 = 0,
dlzVzk - k 2Vzk = — dz Pk
dlVak - k2Va
¡l
ikn
ak
Pk
¡l
dz Vzk + ikf Vfk = 0
d2Pk - k2Pk = 0
z> H,
(11)
(12)
Vzk\z=H-0 = ^^k\z=H+0 , Vak\z=H-0 = Vak\z=H+0 ,
(-Pk + 2idz Vzk)\z=H-0 = (-Pk + 2ll dz Vzk)\z=H+0
l(ikaVzk + dzVak)\z=H-0 = Il(ikaVzk + dzVak)\z=H+0 .
dth0k = Vzk\z=0 , dthk = Vzk\z=H .
It is convenient to separate the longitudinal V1 and transversal V^ components velocity Vak in Fourier space:
(13)
(14)
(15)
(16)
of the
k
V = ka VII
= k k + Vak,
VI = kfV
Vak = kVfk,
V
ak
x kakf I V af--I Vfk.
The equations (10)-(16) can be rewritten in terms of VI and V^:
d2zzVzk - k2Vzk = 1 dzPk
dlVk1 - k2V}}
¡
ik
dlVatk - k2V^
Pk
ak
¡
0
0 < z < H,
(17)
dz Vzk + ikVfk = 0 d2zz Pk - k2Pk = 0
Vzk\
z=0
Vz
zo k,
l(ikVz
zk
z=0
(-Pk + 2ldz Vzk)\z =H-0 = fzHk
l(ikVz
zk + dzVk
f
Hk
z=H-0 Idz Vak I z =H-0 f aHk, dth0k = Vzk\z=0 , dthk = Vzk\z=H
(18)
(19)
d2zzVzk - h2Vzk = — dzPk
dizVÏ - k2V£ = -Pk
dlVl- kV
Jzzy ak
2t/-l ak
Hl
ik
Hi 0
0 < z < H,
(21)
dz Vzk + ikVk = 0
d2zzPk - k2Pk
0
(-Pk + 2hi dz Vzk)lz=
z=H+0
fz
zHk
Hi(ikVzk + dz VÏ)
z=H+0 Hldz Vak \z=H+0 = faHk
fVzk\
f
Hk
(22)
Vk V k Pk
0.
(23)
From the continuation equation (the third one in (17) and (21)) one can obtain:
VÏ
k
dz Vzk.
(24)
Taking into account the equation for Vzk, we get the following boundary conditions for the tangential stress:
(2Hk 2 Vzk + dz Pk)\z=0 = 0,
(2Hk 2Vzk + dz Pk)\=H_ 0 = (2Hik 2Vzk + dz Pk)\ =
(25)
z=H+0
The requirement of normal stress continuity on the interface becomes:
(-Pk + 2idzVzk)\z=n-o = (-pk + 2iidzVzk)\z=n+o ■ (26)
The equations of the evolution of the free boundary and the interface (20) involve only the Vzk term. To obtain it in the layer and half-space domains, we write out the equations and boundary conditions using (17), (21), (23), (25) and (26):
1
d2Vk - k2VZk = -dzPk d2zPk - k2Pk = 0
H
dly7k - k2Vzk = — dz Pk
Hi
d2 Pk — k2Pk
Vzk|
z=0 = Vz0k, (2Hik2Vzk + dz Pk)\z=0
0 < z < H,
z> H,
0,
Vzklz=H-0 = Vzklz=H+0 , dzVzklz=H-0 = dz^^k|z=H+0 ,
(-Pk + 2Hdz Vz
z Vzk)lz=H-0
(-Pk + 2Hldz Vzk)lz=H+0
(2Hk2 Vzk + dz Pk)\ = H_ 0 = (2Hi k2Vzk + dz Pk)\ =
z=H-0
Vzk Pk
z=H+0
0.
(27)
(28)
(29)
(30)
(31)
(32)
(33)
First, consider the half-space domain z > H. To simplify the system, we assume Vzk to be represented in the form:
zH
Vzk = V,k + —-Pk. (34)
2li
In that case:
d2zzVzk - k2Vzk = 0 8lzPk - k2Pk = 0 z > H, (35)
Kk = VzHke-k(z-H ), Pk = P+ke-k(z-H),
V= {V-Hk + ^¿k) e-k'z-H
> ^ k — ^ h k
z - H„ + \ -^-„) (36)
l Hk
In the layer domain 0 < z < H, we assume the velocity Vzk is represented as:
zH
Vzk = Vzk + Pk. (37)
2i
Therefore,
d2zzVzk - k2Vzk = 0 d2zzPk - k2Pk = 0 0 <z<H. (38)
The solution of these equations is:
p p- sinh(kz) sinh(k(H — z))
Pk = PHk sinh(kH) + Pok sinh(kH) :
z — H \ Sinh(kz) ^ z \ Sinh(k(H — z))
Vzk = I VzHk +--~-PHk . i /, m + I Vz0k + 7T"PHk
(39)
iZZ
2i " HV sinh(kH) ' \'z"k ' 2^ Hk) sinh(kH) .
Using (36), (39) and the corresponding boundary conditions, the equations connecting VzHk, Vz0k, P0k, P-k, P+k may be summarized as:
H
(l cosh(kH) + ¡l sinh(kH)) VzHk - IVz0k = — P0k 2(i - ¡l)k sinh(kH)VzHk + cosh(kH)P-k + sinh(kH)P+k = P0k
f H\ 1 1 (40)
ke-k H Vz Hk - klVz0k + — J P0k + — sinh(kH )P- = ^ sinh(kH)P+k
2ik sinh(kH)Vz0k + P-k - cosh(kH)P0k = 0. This system of equations allows one to find the relation between the VzHk and Vz0k:
VzHk = YVz0k (41)
or
hk = Yh0k. (42)
The attenuation coefficient y is the ratio of the amplitudes of the displacement of the interface to the displacement of the external free boundary:
kH(sinh(kH) + m cosh(kH)) + cosh(kH) + m sinh(kH) 1 = kH (m - m) + (cosh(kH) + m sinh(kH)) (cosh(kH) + m sinh(kH)), ( )
here H is the layer thickness, k is the length of the Fourier coordinate vector, m = ¡/¡, ¡ is the layer fluid viscosity, ¡ l is the half-space fluid viscosity. Applying the inverse Fourier transform to (41) and integrating the result over time, the interface deformation can be obtained:
h(x, y) = /fb(\/k2 + ky) h0k(kx, ky)ei(kxX+kyy)j dkxdky (44)
With a given deformation shape h0, one can evaluate the normalized interface deformation amplitude using the explicit expressions:
= 1 r
= (2n)2 AoJ J-„
Y
(y/kl + k2) 11 ho(x, y)e-i(kxX+kyy) dxdy
dk xdky
in the axisymmetric case, and:
1 r
S
2nA
0 J-r
f+œ
7 (\k\) ho(x)e x dx
dk
in the planar case.
2.3. Subsequent analytical results
Using the approach described above, a number of conclusions correlating the parameters of the free-surface perturbation and interface deformation can be made.
(1) Asserting k = 0 in (41) and applying 7\k=0 = 1, one can get:
or
Vzh d2x = / Vz0 d2x
dth0 d2x = / dthd2x.
(45)
(46)
Integrating (46) over the time and asserting zero initial deformations of the free surface and interface, one can obtain the equality of the interface deformation volume and the free surface deformation volume:
J h0 d2x = J h d2x = V. (2) Differentiate (42) by ka and apply k = 0. Taking into account:
(47)
d7
dk
0,
dhk
k=0
one can obtain:
or
dka
—i xah(x) d2x,
dh
0k
k=0
dka
—i xah0(x) d2x
k=0
where:
xah(x) dx = J xah0(x) dx
{xa) {xa)o,
(xa) = 1 xah(x) d2x, (xa)0 = 1 xa^(x) d2x
(48)
are the coordinates of the geometrical centers of the interface and free surface respectively.
(3) Differentiate (42) by ka and hp and apply k = 0. Taking into account
d27
dkadkf
= H 2 Sa f,
d 2hk
k=0
dkadkf
— xaxf h(x) d2x,
k=0
d2h
0k
dkadkf
one can obtain the following:
— xaxfhfl(x) d2x
k=0
xaxf h(x) d2x = xaxf h0(x) d2x + VH2S,
f
or
(xaxf h) = (xaxf K)0 + H2Saf,
(49)
where
^ h) = 11 ^ rn SX, ^ ft)o = / x„x h„(x) A.
Asserting a = /3 = 1 leads to relation of the mean-square width of the deformations on the interface and free surface:
X) = {x\)o + H2 (xl) = x )o + h2 (50)
Consider the squared cross sectional dimensions R2 = xf + x2. Using (50) one can get the following expression for the dimensional change of the perturbation:
(R2) = (R2)o + 2H2. (51)
(4) Integrate (42) over the two-dimensional k-space:
J hk d2k = J y(k)hok d2k.
Taking into account / hk d2k = (2n)2 h|x=0, the next relation can be obtained:
h|x=o =1^1 Y(k)hok d2k. (52)
Consider a free space deformation of the small cross sectional dimensions comparing to the layer thickness H. In this case, the Fourier image of such deformation is localized in a domain of bigger scale than H-1. On the other hand, Y(k) attenuates rapidly with k > H-1. Therefore, (52) can be represented as:
h|x=o = / Y(k) [hok]fc=o d2k = Y(k) d2k. (53)
Since k and H enters in Y(k) as product kH, (53) can be rewritten in the form:
V
H2
where
h|x=o = — (54)
0(m) = i y Uh) d% m = ± (55)
(2 n)2J VH^ ^
(5) Using the expressions ylk=o = 1, §t fc=o = 0, a§k^ ■ o = -H2^, obtained above
we can expand y in Taylor series:
Y =1 - 1 k2 H2 + ... Equation (42) can be expressed in the form:
hk = (l - 1 k2H2 + ..^hok. (56)
Applying the inverse Fourier transform to (56) and truncating the series after the second term one can obtain:
h(x) = (l + ±H2A^ ho(x), (57)
where A = df + d% is the Laplace operator by the plane coordinates. Formula (57) is valid when the perturbation width is much larger than the layer thickness H.
As an example of (57), consider the deformation of the free surface described by h0(x) = Aexp(-alx2 - a2x2), alH2 c 1, a2H2 c 1. In this case, the interface deformation becomes:
h(x) = (1 + H2[al(2alx\ - 1) + a2(2a2x2 - 1)]) Aexp(-alx\ - a2x2).
It can be shown that the interface deformation amplitude is smaller than the one for the free surface:
h\x=0 = [1 - H2(a\ + «2)] h0U .
On the other hand, with \xl\ > and \x2\ > , the interface deformation is larger than the one for the free surface.
3. Numerical model
To obtain a numerical solution, COMSOL Multiphysics was used. Two-dimensional planar and axisymmetric models were implemented. These models have only minor differences in their setups. Both use Laminar Flow and Moving Mesh interfaces.
A domain of layer and half-space is described by two rectangles with different Fluid Properties. Symmetry boundary conditions on the side faces are set (Axial Symmetry and Symmetry in axisymmetric case). The slip boundary condition is applied along other boundaries.
The deformation of the bottom boundary is prescribed by an analytical expression which is linear with respect to time (Prescribed Mesh Displacements is used). The side faces are free for vertical deformation and locked for horizontal. Displacement of the interface is described by Prescribed Mesh Velocity: velocity of the mesh is equal to the local fluid velocity.
Since we model incompressible flow in a small finite domain, we have to keep the volume to satisfy mass conservation. In the case of the two-dimensional problem, it can be said that we have to keep the area of the domain. There are three ways to solve this problem. The first one is to choose deformation function which satisfies the conditions described above. For example, it can be a sinusoidal function. The second option is to add negative deformation far enough from the area of interest. Since the phenomenon under study is local, its effect should be negligible. The third option is to allow the top border deformation, i.e. assume the free boundary condition on the top border. Again, since the perturbation is local, the effect of the free boundary shouldn't be significant.
The parameters of the model are layer thickness, layer/half-space viscosity ratio and the boundary deformation shape. The primary result of the model is the shape of the interface after the deformation.
Parametric Sweep was used to compute the model for different layer thicknesses and viscosities. The half-space thickness was kept constant at 15 m. The layer thickness was varied from 0.5 m to 2 m in 0.1 m increments. Viscosity ratios of 0.1, 1 and 10 were
considered. The free surface perturbation was described by functions h0(x) = A0 cos ^20x) (only for planar case) and h0(x) = A0e xo .
Linearity with respect to perturbation amplitude The analytical model asserts that for small perturbations, the normalized interface deformation doesn't depend on the surface defect amplitude. Numerical computation shows that this is correct, even for quite large amplitudes compared to the layer thickness (fig. 4).
Normalized interface amplitude A0 = 0.05 m (solid), A0 = 0.25 m (dashed)
0.4
\\
V
0.05
0
0.50 0.70 0.90 1.10 1.30 1.50 1.70 1.90
Layer thickness, m
Fig. 4. Normalized interface deformation due to different amplitudes of the free surface perturbation (axisymmetric model)
3.1. Result comparision
A number of numerical experiments were made for different parameters of the free surface perturbation, layer thickness and viscosity ratio. As a result, the dependences of the interface perturbation amplitude on the layer thickness was obtained. First, consider
planar harmonic free surface perturbation described by h0(x) = A cos ^. Since the
Fourier transform of this function is:
hok = An [S(k - 2n/Ac) + 5(k + 2n/Ao)] ,
it's easy to find analytical expression for the interface displacement h(x):
h(x) = ho(x)Y ^■
Comparison with the numerical results gave very good agreement (fig.5). On all figures below, dashed lines correspond to the analytical approach and dots correspond to the numerical solution.
Fig. 5. Interface perturbation due to harmonic deformation of the free surface (planar case)
In the case of Gaussian deformation h0(x) = A0e xo , the analytical approach remains accurate even for large deformation amplitudes A0 (figures 6-7).
4. Conclusion
Numerical and analytical models were developed to analyze the propagation of the free surface perturbation in the layered viscous half-space. Analytical formulas for perturbation amplitude and width on the interface between the half-space and the layer were obtained. The results may be summarized as follows:
A x ivy [ ■ n t K" t M i.
Aj=0,Z5 m, \, = 0.5 m
fir
I.
s 1 2
-D <3J
| s
z o.i
l*i -
„
V *
s e -
-t *
-••
ato a» 0« Layer thickness m 1 »0
2
Fig. 6. Interface perturbation due to Gaussian deformation of the free surface (axisymmetric case)
Fig. 7. Interface perturbation due to Gaussian deformation of the free surface (planar case)
• Perturbation of the interface between the layer and the remaining half-space decreases rapidly with the layer thickness increase
• Propagation depends on the ratio of the viscosities: higher viscosity of the halfspace fluid helps to reduce perturbation on the interface
• Propagation depth increases with the wavelength of the initial perturabtion on the free surface
• For the same wavelength, a smaller amplitude of the the initial perturbation causes a smaller deformation of the interface. Normalized perturbation remains almost the same.
References
[1] Happel S.J. and Brenner H. Low Reynolds Number Hydrodynamics. Upper Saddle River: Prentice-Hall, 553 p. (1965).
[2] Chivilikhin S. and Amosov A. Planar Stokes flows with free boundary. In book "Hydrodynamics -Advanced Topics". Book edited by: Dr. Harry Edmar Schulz, Andre Luiz. P. 77-92. Intech (2011).
[3] Hopper R. W. Stokes flow of a cylinder and half-space driven by capillarity. J. Fluid. Mech., 243, P. 171-181 (1992).
[4] De Voeght F. and Joos P. Damping of a disturbance of a liquid surface. Journal of Colloid and Interface Science, 98(1), P. 20-32 (1984).
[5] Mullins W. W. Flattening of a nearly planar solid surface due to capillarity. J. Appl. Phys., 30(77), P. 77-83 (1959).
[6] Lamb H. Hydrodynamics. Cambridge Univerity press, Cambridge, 768 p. (1993).
[7] Levich V. G. Physicochemical hydrodynamics.Longman, Harlow, 700 p. (1962).
[8] Zatsepin A.G., Kostyanoi A. G., and Shapiro G.I. Slow spreading of a viscous liquid over a horizontal surface. Dokl. Akad. Nauk SSSR. 265(1), P. 193-198 (1982).
[9] Jeong J. T., Moffatt H. K. Free-surface cusps associated with flow at low Reynolds number. J. Fluid Mech., 241, P. 1-22 (1992).
[10] Chivilikhin S.A. Relaxation of small perturbations of highviscous liquids planar surface. Nanosystems: Physics, chemistry, mathematics, 3(4), P.54-66 (2012).
[11] Chivilikhin S.A. Relaxation of a small local perturbation of the surface of a viscous fluid in the Stokes approximation. Mekhanika Zhitdkosti i Gaza, 3, P. 133-137 (1985).
