Journal of Samara State Technical University, Ser. Physical and Mathematical SciencesJournal of Samara State Technical University, Ser. Physical and Mathematical Sciences1991-86152310-7081Samara State Technical University4199610.14498/vsgtu1763Original ArticleAn undamped oscillation model with twodifferent contact angles for a spherical dropletimpacting on solid surfaceChenShiPhD, Associate professorCongBozhongZhangDongqiLiuXiaohuaShenShengqiangDalian University of Technology3107202024239040004082020Copyright © 2020, Samara State Technical University2020In order to further elucidate the dynamic theory of droplet oscillating on solid surface, a new handling method of contact angle of the droplet during the process of the oscillation was founded, which is based on the spherical model. The influence of gravity on the contact angle andspreading radius was discussed. Thus, an equation between the spreading radius of the dropletand time flow was founded. The results of theoretical calculation were compared with smoothednumerical results.droplet oscillationcontact angletheoretical analysisspectrum analysisколебание каплиугол контактатеоретический анализспектральный анализ\centerline{\bf\large 1. Introduction} \smallskip As a common phenomenon in nature, droplets impacting on solid surfaces exists in many fields, such as engineering, chemical industry, agriculture, aerospace and energy. For examples, the liquid fuel injection in the internal combustion engine [1], the ink jet of the printer [2], the process of seawater evaporation in desalination system [3], the spraying of the refrigerant in the evaporator in the refrigeration system [4], etc. The study of relationships among the dynamic parameters, initial parameters of droplets and initial can be beneficial to predict the processes of droplets impacting on solid flats. For a droplet impacting on solid surfaces, according to the initial parameters, there may be oscillations, bounces, splashes, etc. Compared with other circumstances, the oscillations have more value and possibility to be studied. Since T. Young [5] presented that the contact angle between liquid and solid surface was constant, and P.S. Laplace [6] inferred the relationship between additional pressure and radius of curvature, many studies have been done by experiments, numerical simulations and theoretical analysis. In the aspect of experimental studies, M. Marengo et al. [7] studied the effects of the impact parameters on the droplet impingement, and found that the surface wettability had a strong influence on the spreading of droplet in the later stages of the process. I.S. Bayer [8] studied the dynamic characters of contact angle between smooth surface and droplets with different wettability. M. Remer et al. [9] studied the variation of droplet dynamical contact angle after impacting on three kinds of coating surfaces. In the aspect of numerical simulation, S.F. Lunkad et al. [10] simulated the drop impact and spreading process on horizontal and inclined surfaces using the volume of fluid (VOF) method, and investigated the effects of surface inclination, surface wetting characteristics, liquid properties and impact velocity on the droplet oscillation by using static contact angle (SCA) and dynamic contact angle (DCA) models. Y. Yao et al. [11] analyzed droplets oscillation with VOF method and introduced a model of dynamical contact angle to improve the accuracy of simulation. S. Sikalo et al. [12] carried out the numerical simulations of a single drop impacting onto a dry, partially wettable substratum. In the aspect of theoretical analysis, S. Vafaei et al. [13, 14] explained the dependence of contact angle on the size of liquid droplets on smooth solid substrates, and demonstrated that for sessile droplets on smooth surfaces, the contact angle can be uniquely determined for given droplet mass (or volume) and liquid/solid/gas properties. I.V. Roisman et al. [15] studied the normal impact of a liquid drop on a dry solid surface theoretically, and introduced a strictly theoretical model to predict the evolution of the drop diameter. The purpose of this paper is to establish an amended theoretical model of undamped droplet oscillation.Compared with original model, in which droplet contact angle keeps constant value, the droplet contact angle in new model varies with spreading radius and droplet volume. The waveform of droplet spreading radius obtained from new model, original model and numerical simulation were compared. % \bigskip \newpage \centerline{\bf\large 2. Theoretical Model} \smallskip \centerline{ \slshape 2.1. The Oscillating Equation of Spherical Segment Droplet} \centerline{ \slshape with Changing Contact Angle} \smallskip For a droplet oscillating on a horizontal solid surface, when the size of droplet is short enough that the influence of gravity on the droplet shape can be ignored. The shape of the droplet can be considered as a sphere, shown in Fig.~\ref{shichen:fig1}. \begin{figure}[h!] \centering \includegraphics[scale=1]{shichen_fig1} \caption{The geometrical model of droplet \label{shichen:fig1}} \end{figure} At a certain moment of oscillation, spreading radius of droplet in horizontal direction is $r$, the radius of the sphere is $R$, the height is $z$ and the contact angle between the droplet and the solid surface is $\theta$. According to the geometrical relationship, the equation between $r$, $R$, and $\theta$ can be written as: \begin{gather} \label{shichen:eq1} \frac{3V}{\pi r^3} =(1-\cos \theta)^2(2+\cos\theta), \label{shichen:eq3} R=\frac{r}{\sin \theta}, \label{shichen:eq3} z= r \frac{1-\cos \theta}{\sin \theta}, \end{gather} where $V$ is the droplet volume. When a droplet impacts on a solid surface, the shape of droplet is strongly irregular and hard to describe in beginning several oscillation periods after impact. But the oscillation will be gradually stabilized with close periods and amplitudes. Thus, it can be assumed that the droplet keeps the shape of spherical cap in the hereafter oscillation. In previous study [16], the oscillation of droplet could be seen as the result of interaction of the surface tension, the internal pressure and the inertial force. A differential element is selected in the segment droplet of the angle, shown in Fig.~\ref{shichen:fig2}. \begin{figure}[h!] \centering \includegraphics[scale=1]{shichen_fig2} \caption{The force analysis of droplet differential. The surface tension and the internal pressure equally and vertically act on each vertical section of the droplet. The surface tension points to the inside of differential element and the internal pressure points to the outside of differential element\label{shichen:fig2}} \end{figure} And the oscillating equation of spherical segment droplet can be written as~[16]: \begin{multline} \label{shichen:eq4} \Bigl[ - \frac{\pi \rho R^4\theta}{4r} + \frac{\pi\rho (R-z)}{12}(2r^2+3R^2) \Bigr] \frac{d^2r}{dt^2} - - 2\pi\sigma (l-r\cos \varphi) + \frac{4\pi\sigma A}{R} + 2\pi\rho g \Bigl( RA -\frac{r^3}{3} \Bigr)=0, \end{multline} where $\varphi$~--- the intrinsic contact angle of the droplet; $\theta$~--- the dynamic contact angle of the droplet;\footnote{The angle $\varphi$ is determined by the involved surface energies. When the gravity is absent and the droplet is static on the horizontal solid surface, $\varphi=\theta$. The angle $\theta$ can be influenced by inertial force and gravity of the droplet.}\!\! $\sigma$~--- the surface tension coefficient; \begin{equation} \label{shichen:eq5} l =RA \end{equation} --- arc length of the great circle of droplet; \begin{equation} \label{shichen:eq6} A=\frac12 \bigl[Rl - r(R-z) \bigr] \end{equation} --- the sectional area of the droplet in the vertical plane. Eq. \eqref{shichen:eq4} can be also written as: \begin{gather} \label{shichen:eq7} P \frac{d^2r}{dt^2} +Q + W=0 , \label{shichen:eq8} P =- \frac{\pi \rho R^4\theta}{4r} + \frac{\pi\rho (R-z)}{12}(2r^2+3R^2), \label{shichen:eq9} Q = - 2\pi\sigma (l-r\cos \varphi)+\frac{4\pi\sigma A}{R} , \label{shichen:eq10} W= 2\pi\rho g \Bigl( RA -\frac{r^3}{3} \Bigr), \end{gather} where $P$~--- the factor of inertia, $Q$~--- the factor of surface tension and internal pressure, $W$~--- the factor of gravity. When substituted equations of geometrical relationship (Eq.~\eqref{shichen:eq1}--\eqref{shichen:eq3}, and Eq.~\eqref{shichen:eq5}, \eqref{shichen:eq6}), $P$, $Q$, and $W$ can be written as: \begin{gather} \label{shichen:eq11} P = \pi \rho r^3 \Bigl( \frac{\cos\theta}{6\sin \theta} +\frac{\cos\theta}{4\sin^3\theta} \frac{\theta}{4\sin^4\theta} \Bigr), \label{shichen:eq13} Q= 2\pi\sigma r (\cos\varphi -\cos\theta), \label{shichen:eq13} W= \pi \rho g r^3 \Bigl(\frac{\theta}{\sin^3\theta} -\frac{\cos \theta}{\sin^2\theta} -\frac 23 \Bigr). \end{gather} In this model, the intrinsic contact angle is a constant valve, which determines direction of the contact force from the solid surface. And the dynamic contact angle determines the shape of droplet. The geometrical relationship between them is shown in Fig.~\ref{shichen:fig3}. \begin{figure}[h!] \centering \includegraphics[scale=1]{shichen_fig3} \caption{The geometrical relationship between the intrinsic contact angle and the dynamic contact angle \label{shichen:fig3}} \end{figure} The deviation between the intrinsic contact angle and the dynamic contact angle increases with the increase of acceleration of the droplet. And the expression of the contact force from the solid surface in Fig.~\ref{shichen:fig2}, ``$\sigma\cos\theta$'' can be replaced by ``$\sigma\cos\varphi$''. \smallskip \centerline{ \slshape 2.2. The Influence of Gravity} \smallskip When the spreading radius $r$ reaches the balance value (or the droplet keeps static): \begin{equation} \label{shichen:eq14} P \frac{d^2r}{dt^2} +Q + W=0. \end{equation} Substituted Eq. \eqref{shichen:eq11}--\eqref{shichen:eq13} and multiplied: \begin{equation} \label{shichen:eq15} \cos \varphi =\cos\theta - \frac{\rho g r^2}{2 \sigma} \Bigl( \frac{\theta}{\sin^3\theta} - \frac{\cos\theta}{\sin^2\theta}-\frac{2}{3} \Bigr). \end{equation} Obviously, in the absence of gravity, $\theta=\varphi$. In normal gravity condition that the droplet can still keep the shape of spherical segment (called micro gravity in followings), the dynamic contact angle $\theta$ can be calculated by iteration from Eq.~\eqref{shichen:eq1} and Eq.~\eqref{shichen:eq5}. The influence of gravity is shown in Fig.~\ref{shichen:fig4}. \begin{figure}[p!] \centering \includegraphics[scale=.9]{shichen_fig4} \caption{(top) the influence of droplet volume on contact angle $\theta$ (here, the gravity of droplet depends on the droplet volume); (bottom) the influence of gravity on balance radius \label{shichen:fig4}} \bigskip \includegraphics[scale=.9]{shichen_fig5} \caption{The frequency spectrum of numerical data \label{shichen:fig5}} \end{figure} With the increase of droplet volume, the deviation between the micro gravity solutions and the absence gravity solutions on dynamic contact angle and spreading radius becomes larger, but that of spreading radius is quite small. \newpage \centerline{\bf\large 3. Numerical Simulation} \smallskip \centerline{ \slshape 3.1. Original Case} \smallskip For the numerical solution,the processes of drop oscillation are simulated by using the VOF model of Fluent. The body force weighted model is used to calculate of pressure and gravity balance. The PISO algorithm is used to couple the droplet velocity and pressure in the momentum equation. The slip boundary condition is used and the shear stress is 0~Pa. The viscosity of water is $1.003\cdot 10^{-5}$~Pa${}\cdot{}$s which is multiplied by~0.01 in order to simulate the undamped oscillation of droplet. The time step is $4\cdot 10^{-6}$~s, and the residual error is $1\cdot 10^{-5}$. Others keep the default algorithm. \smallskip \centerline{ \slshape 3.2. Spectral Analysis} \smallskip The data from original case contains much ingredient of noise. This is because the oscillation is influenced by many minimal factors. After eliminating several data points with extreme deviation, the data was resampled by FFT operation. The spectral analysis of the original data is shown in Fig.~\ref{shichen:fig5}. In Fig.~\ref{shichen:fig5}, the main frequency of the oscillation is concentrate upon [250~Hz, 550~Hz]. The corresponding frequency value of the peak value in figure is 426~Hz, which means the periods of the oscillation is about~2.35~ms with the same initial parameters shown in Tab.~\ref{shichen:tab1}. \begin{table}[h!] \small \centering \caption{The initial parameters of theoretical and numerical calculation\label{shichen:tab1}} \vspace{-3mm} \begin{tabular}{c|c|c|c|c} \hline \footnotesize Initial Speed \footnotesize Droplet Volume \footnotesize Gravity \footnotesize Density \footnotesize Surface Tension \hline \rule{0mm}{12pt}% 0.8 m/s $1.1310 \cdot 10^{-10}$ m$^3$ 9.8 m/s 1003 kg/m$^3$ 0.073 N/m \hline \end{tabular} \end{table} \smallskip \centerline{ \slshape 3.3. Data Smoothing} \smallskip To eliminate the noise, a Butterworth filter was applied. The response type is bandpass with the pass band [400~Hz, 500~Hz], the fluctuate is less than 1~dB. On the both side of pass band, the signal is decreased to 10~dB. The frequency spectrum after smoothing is shown in Fig.~\ref{shichen:fig6}. \begin{figure}[p!] \centering \includegraphics[scale=.9]{shichen_fig6} \caption{The frequency spectrum after smoothing \label{shichen:fig6}} \bigskip \centering \includegraphics[scale=.9]{shichen_fig7} \caption{The relationship between spreading radius and time flow ($\zeta=(r-r_b)/r_b$, $r_b$ is the theoretical balance radius) \label{shichen:fig7}} \bigskip \centering \includegraphics[scale=.9]{shichen_fig8} \caption{The relationship between spreading radius and its accelerated velocity \label{shichen:fig8}} \end{figure} After smoothing, the data can be regarded as a superposition of several sine curves, with close frequency, amplitude and phase position. \bigskip \centerline{\bf\large 4. Comparison} \smallskip Eq.~\eqref{shichen:eq7} is solved by using the fourth-order Runge--Kutta method in MATLAB R2016a. The relationship between spreading radius and time flow of theoretical and numerical results (after smoothing) is shown in Fig.~\ref{shichen:fig7}, with the same initial parameter shown in Tab.~\ref{shichen:tab1}. The period of theoretical results is 2.37~ms, which is almost equal to the average periods of numerical results. The amplitude of theoretical results is about 10\,\% larger than the average value of the numerical results. Compared with constant contact angle model [16], the average periods were more close to the numerical data (1.8~ms). The relationship between spreading radius and its accelerated velocity is shown in Fig.~\ref{shichen:fig8}. When $\zeta$ is approaching or higher than the balance value (0), the deviation between theoretical and numerical results is relatively small, but when $\zeta$ is far below the balance value, the deviation is considerable. This is because an obvious distortion would happen on the shape of droplet and the above-mentioned relationships (Eq.~\eqref{shichen:eq11}--\eqref{shichen:eq13}) would lose efficacy. Besides, because the relationship between spreading radius and its acceleration in present model is not linear, the oscillation period is related to the initial velocity. \bigskip \centerline{\bf\large 5. Discussion} \smallskip The process of undamped oscillation for a droplet impacting on solid surface was theoretically described by an equation based on spherical segment model. Followings are several conclusions. The contact of the droplet is changing in the oscillation and geometrically calculated. The influence of gravity on contact angle and spreading radius was analyzed quantitatively, the deviation was lower than the constant contact angle model. However, the result was based on the geometrical model of spherical segment of droplet, which might be inaccurate under the influence of gravity. More study was needed in this content. After spectrum analyzing and smoothing by a filter, the numerical solutions can be regarded as a superposition of several sine curves with close frequency, amplitude and phase positions. The relationship between spreading radius and its accelerated velocity was founded by polynomial fit, and the curve was symmetrical about the balance point and strongly linear. By comparing the results of theoretical and numerical calculation, when the spreading radius of the droplet approaches the balance value, the deviation between two solutions are small, but when the spreading radius is far from the balance value, the deviation can be considerable. And the deviation of far below from the balance value is larger than that of far above from the balance value. This is possibly because when the spreading radius is far from the balance value, the shape of the droplet would also deviate from spherical segment. More work about the amendment of the shape of the droplet in oscillation is needed in later study.[Nakayama Y., Kidokoro T., Sakurai K., Fuel injection control system of an internal combustion engine, US Patent no. US9169758B2, 2015][Slater S.D., Clippingdale A.J., Newcombe G.C.F., Printing process and liquid ink jet ink, US Patent no. US9156256B2, 2015][Qi C.H., Feng H.J., Lv H.Q., Miao C., "Numerical and experimental research on the heat transfer of seawater desalination with liquid film outside elliptical tube", Int. J. Heat Mass Transfer, 93 (2016), 207-216][Hartfield J.P., Sanborn D.F., Falling film evaporator with refrigerant distribution system, Canada Patent no. CA2219676A1, 1995][Young T., "An essay on the cohesion of fluids", Phil. Trans. Roy. Soc. London, 95 (1805), 65-87][Laplace P.S., "Sur l'action capillaire. Supplement à la theorie de l'action capillaire", Traite de mecanique celeste, v. 4, Supplement 1, Livre X, Gauthier–Villars et fils, Paris, 1805, 771–777][Šikalo Š., Marengo M , Tropea C., Ganic E.N., "Analysis of impact of droplets on horizontal surfaces", Experimental Thermal and Fluid Science, 25:7 (2002), 503-510][Bayer I. S., Megaridis C. M., "Contact angle dynamics in droplets impacting on flat surfaces with different wetting characteristics", J. Fluid Mechanics, 558 (2006), 415-449][Remer M., Psarski M., Gumowski K., Rokicki J., Sobieraj G., Kaliush M., Pawlak D., Celichowski G., "Dynamic water contact angle during initial phases of droplet impingement", Colloids and Surfaces A: Physicochemical and Engineering Aspects, 508 (2016), 57-69][Lunkad S. F., Buwa V. V., Nigam K.D.P., "Numerical simulations of drop impact and spreading on horizontal and inclined surfaces", Chem. Eng. Sci., 62:24 (2007), 7214-7224][Yao Y., Meng S., Li C., Chen X., Yang R., "Droplet oscillation after impact on a solid surface", International Mechanical Engineering Congress and Exposition, 7, Fluids Engineering (2016), IMECE2016-66025][Šikalo Š. , Wilhelm H.-D., Roisman I. V. , Jakirlic S., Tropea C., "Dynamic contact angle of spreading droplets: Experiments and simulations", Phys. Fluids, 17:6 (2005), 062103][Vafaei S., Podowski M. Z., "Theoretical analysis on the effect of liquid droplet geometry on contact angle", Nuclear Eng. Design, 235:10–12 (2005), 1293-1301][Vafaei S., Podowski M. Z., "Analysis of the relationship between liquid droplet size and contact angle", Adv. Colloid Interface Sci., 113:2–3 (2005), 133-146][Roisman I. V., Rioboo R., Tropea C., "Normal impact of a liquid drop on a dry surface: model for spreading and receding", Proc. Royal. Soc. A, 458 (2002), 1411–1430][Chen S., Zhang D., Shen S., Liu X., Chen Y., "Spherical drop impact on solid surfaces: Un-damped oscillation theoretical model", AIP Conf. Proc., 1984:1 (2018), 020032]