Investigation of behavior of the dynamic contact angle on the basis of the Oberbeck-Boussinesq approximation of the Navier-Stokes equations

Flows of a viscous incompressible liquid with a thermocapillary boundary are investigated numerically on the basis of the mathematical model that consists of the OberbeckBoussinesq approximation of the Navier-Stokes equations, kinematic and dynamic conditions at the free boundary and of the slip boundary conditions at solid walls. We assume that the constant temperature is kept on the solid walls. On the thermocapillary gas-liquid interface the condition of the third order for temperature is imposed. The numerical algorithm based on a finite-difference scheme of the second order approximation on space and time has been constructed. The numerical experiments are performed for water under conditions of normal and low gravity for different friction coefficients and different values of the interphase heat transfer coefficient.


Introduction
The problems of the fluid flows with a dynamic contact angle are very important at present. Such problems occur in many applications by study of the convective multi-phase flows, flows through porous medium, of evaporative liquid drops, wetting hydrodynamics etc. (see, for example, [1]). The contact angle is defined as an angle between the two surfaces, where the free boundary touches the solid boundary. The non-stationary problems remain to be rather difficult for investigations in connection with an open problem of the dynamic contact angle, due to an appearing problem of motion of the contact line of three phases, which does not have unique closure (see [2][3][4][5][6][7]). The main question of such types of the problems relates to a non-consistency of the condition on the free boundary and the classical no-slip boundary condition on the solid wall in a neighborhood of the contact point. From the mathematical point of view this non-consistency has been considered in [8].
If the boundary walls of a liquid containing domain are moving with a constant speed, the contact angle is changing with the velocity [4,5,9]. Changing of this angle will also depend on a character of thermal regimes on the free and solid boundaries. Interaction of the gravitational and thermocapillary convection and their influence on behavior of the contact angle and shape of the free boundary should be investigated.

Statement of the problem
Two-dimensional motion of a viscous incompressible liquid with a thermocapillary boundary is investigated numerically on the basis of the mathematical model that consists of the Oberbeck-Boussinesq approximation of the Navier-Stokes equations, kinematic and dynamic conditions at free boundary [5,9]. The slip conditions (conditions of proportionality of the tangential stress to the difference of the tangential velocities of liquid and wall) are prescribed on the solid boundaries of the channel supporting by constant temperature. On the thermocapillary gas-liquid interface the condition of the third order for temperature is imposed. The problem is studied in the quasi-stationary formulation when the coordinate system is moving with the liquid. The flow is assumed to be symmetric that has to be taken into account in determination of the boundary conditions. Let ȍ be the flow region with the boundary ȍ which consists of three solid ī0, īs and one open part īf (free thermocapillary boundary). Herewith the Cartesian coordinate system is chosen so that the gravity acceleration vector g is directed along the longitudinal axis (here Ox, see Fig. 1a). Let the linear size of the flow domain in the y-direction l be the characteristic length. We introduce the characteristic values for the problem of the fluid flow as follows: u* is the characteristic velocity, T* is the characteristic characteristic temperature drop, p* = ȡu*Ȟ/l is the characteristic pressure, ȡ is the density, Ȟ is the coefficient of kinematic viscosity.
After reformulation the boundary value problem in the terms of new unknown functions ȥ (stream function) and Ȧ (vorticity) we will have the following governing equations: is the velocity vector, T is the temperature, Re = ȣ*l/Ȟ is the Reynolds number, Gr = ȕT*gl 3 /Ȟ 2 is the Grasgof number, Pr = Ȟ/Ȥ is the Rrandtl number, Ȥ and ȕ are the coefficients of thermal diffusivity and thermal expansion, respectivelly, g=|g|.
T = Ts on Ƚs at y = 0 and at y = 1, T = T0 on Ƚ0 (x = x0), Here x=f(y) is the parametrization of the free boundary, , Ca = ȡȞȣ*/ı0 is the capillary number, ı is the non-dimensional surface tension (ı = 1 -MaCa(T-T0)), S is the contact point velocity, ț is the thermal conductivity, ı0 is the surface tension value (dimensional) at some reference temperature, ıT is the temperature coefficient of surface tension, Ȗ, Ȗ0 are the friction coefficients of the liquid and solid walls, į is the interphase heat transfer coefficient. Ts, T0 and Tex are the given constant values of the temperature at solid walls and on the outside. The indexes n and τ denote the derivatives with respect to the normal and tangential directions. The conditions (2)-(6) written in the terms of ȥ and Ȧ are: (2) is the nonpermeability condition including the kinematic condition at free boundary, (3) are the dynamic conditions at free boundary, (4) are the slip conditions on solid walls, which follow from the assumption of momentum conservation and denote the proportionality of shear stresses to the tangential velocity on the walls). The thermal boundary regime is determined by (5) and (6). It is assumed that the constant temperature (5) is kept on the solid walls, on the free boundary the heat exchange with the atmosphere is prescribed by the Newton's law (6).

Numerical results
The numerical algorithm for the problem (1)-(6) based on a finite-difference scheme of second-order approximation in time and space has been constructed (see [5,7]). This algorithm includes computation of the free boundary position f and the initial position f0 of Γf defined by a rest state (at S=0, v=0, T=const) at a given value of the static contact angle Ɏ0. This value Ɏ0 is determined by the type of a liquid, material of the solid wall and external gas medium. Starting with an approximation of f0 we fix a sufficiently small ΔS>0 and organize "step-by-step" computations to reach a necessary value of S. In each step we organize the iteration processes to compute the unknown functions ȥ, Ȧ, T.  [1]). The other value (Nu=17) corresponds to the interphase heat transfer coefficient (δ) that exceeds the values in [1], but allows us to compare the effects of intensity of the thermal regime at a free boundary.
The vortex structure of the flow (see Fig. 1a) and the behaviour of the free boundary (Fig. 1b) are investigated. We observe two-vortex symmetric flow structure. For higher values of the contact point velocity S the flow topology is characterized by displacement of the vortex centres to the free boundary (as shown in Fig. 1a). The flow picture (Fig. 1a) and the examples of the free boundary shape (Fig. 1b) are presented in the case of thermal boundary regime characterized by the nondimensional values Ts = 8, T0 = 1 (or of type "8x1") and Tex = 8, Nu=0.17 (see (5), (6)) under conditions of low (g=9.81(cm/s 2 ); Gr=230) and normal (g=981(cm/s 2 ); Gr=23000) gravity.
The influence of changing of the values of friction coefficients γ0, γ coefficients (γ0 = γ = 100 and γ0 = γ = 1 (g/(cm 2 s))), gravity acceleration g (g=9.81 and g=981 (cm/s 2 )) and Nusselt number Nu (Nu=0.17 and Nu=17) on the dynamics of the contact angle has been studied numerically (Fig. 2). The presented results demonstrate the differences in the contact angle behavior with respect to the different values of the contact point velocity, friction coefficients γ0 = γ = 100 and γ0 = γ = 1 (g/(cm 2 s)), gravity acceleration g=981 and g=9.81 (cm/s 2 ) and to an intensity of the thermal boundary regimes (Fig. 2a, 2b). We observe the increasing character of the contact angle dependence on the contact point velocity. In the case of investigating boundary regimes with rather high values of the external temperature Tex = 8 and high temperature on the lateral walls īs (Fig. 2a) we have found slight discrepancies in the contact angle dynamics under conditions of normal (Fig. 2a, dashed and dashthree-dots lines) and low (Fig. 2a, solid and dash-dot lines) gravity. The similar behavior is observed both for Nu=0.17 and Nu=17 in Fig. 2a and in the case of more moderate thermal regime presented in Fig. 2b for Nu=17 (here Tex = 1 and the thermal regime with Ts = 4, T0 = 1 or of type "4x1" on solid walls has been modeled). At the same time we can speak about more intensive increase of the dynamic contact angle with respect to the contact point velocity in the case of larger values of the friction coefficient (here for γ0 = γ = 100 in comparison with γ0 = γ = 1 (g/(cm 2 s)), see Fig. 2b).