A finite-element model for computing fluid flow inside a sessile evaporating droplet on a solid surface
AA finite-element model for computing fluid flow inside a sessileevaporating droplet on a solid surface
Manish Kumar, Rajneesh Bhardwaj ∗ Department of Mechanical Engineering, Indian Institute of Technology Bombay, Mumbai, 400076, India
Abstract
A finite element model was developed to compute the fluid flow inside a sessile evaporating dropleton hydrophilic substrate in ambient conditions. The evaporation is assumed as quasi-steady andthe flow is considered as axisymmetric with a pinned contact line. The Navier-Stokes equationsin cylindrical coordinates were solved inside the droplet. Galerkin weight residual approach andvelocity pressure formulation was used to discretise the governing equations. Six node triangularmesh and quadratic shape functions were used to obtain higher accuracy solutions. Radial velocityprofiles in axial directions calculated by the FEM solver were compared with a existing analyticalmodel and were found in excellent agreement. The contours of velocity magnitude and streamlinesshow the characteristic flow i.e. radially outward inside the evaporating droplet.
1. Introduction
In recent times, investigation of sessile droplet evaporation has been an interesting area of re-search owing to its many industrial applications such as inkjet printing, spray cooling, makingbioassay etc. To understand heat flow inside an evaporating droplet of pure liquid or predict parti-cle deposition pattern in evaporating colloidal droplet, flow inside the droplet needed to be studiedand calculation of flow field becomes the necessity. There are two modes of droplet evaporation;(1) constant contact radius (CCR) mode, in which droplet contact line remains pinned during theevaporation, (2) constant contact angle (CCA) mode, in which contact angle of the droplet remainsconstant and contact line recedes.Ghasemi and Ward [1] experimentally investigated the e ff ect of thermal capillary flow andenergy transported by it inside droplet during evaporation and reported it a dominating mode ofenergy transport near the three-phase line of the droplet, but at the apex of the droplet, conduc-tion was the dominating mode of energy transfer. Using COMSOL Multiphysics V4.3a software,Barmi and Meinhart [2] simulated the convective flow inside the pinned evaporating droplet whileassuming axisymmetric flow field, temperature, and vapor concentration distribution. They con-cluded that flow inside the droplet is responsible for particle deposition patterns but it does notinfluence heat transfer much and conduction remains the major mode of heat transfer. ∗ Corresponding author
Email address: [email protected] (Rajneesh Bhardwaj )
Preprint submitted to arXiv September 8, 2020 a r X i v : . [ phy s i c s . f l u - dyn ] S e p n the present work, we consider di ff usion-limited and quasi-steady droplet evaporation in CCRmode and presents modeling to calculate the flow inside the evaporating droplet. As it reportedin previous studies [3, 4], droplet evaporates unevenly along the liquid-gas interface, minimum atthe top and maximum at the edge. For a micro-liter size droplet surface tension force dominates, ittries to maintain its spherical cap shape. So when there is more loss of liquid at the edge compareto top due to evaporation, droplet shape gets changed and in order to retain its spherical shape,fluid from the top rushes to edge side to compensate the loss. This process generates the fluid flowinside the droplet. The flow inside the droplet is also responsible for di ff erent types of depositionpatterns on complete drying of the droplet [5, 6].
2. Modeling
In this section, the Finite Element Modeling (FEM) of continuity and Navier-Stokes equationwill be explained to calculate flow field inside the evaporating droplet. As it explained in previoussection, droplet evaporates unevenly along the liquid-gas interface, minimum at the top and max-imum at the edge. For a micro-liter droplet surface tension force dominates, it tries to maintainits spherical cap shape. So when there is more loss of liquid at the edge compare to top due toevaporation, droplet shape gets changed and in order to retain its spherical shape, fluid from thetop rushes to edge side to compensate the loss. This process generates the fluid flow inside thedroplet.
When droplet is very small and fluid flow is very slow, in that case, surface tension forcesdominates over shear and normal stress in the droplet and retains a spherical cap shape. For thesmall droplet of size 1 mm and height of 0.4mm and having the flow velocity of the order of 1 µ m / s, the Bond number (accounts for gravitational force versus surface tension force) is 0.04 andcapillary number (accounts for viscous stress versus surface tension forces) is the order of 10 − [3].Therefore, the approximation of spherical cap is justified. Reynolds number is also small ( Re = . r , z ) cylindrical coordinates as follows:1 r ∂ ( ru ) ∂ r + ∂ v ∂ z = r ∂ ( r σ rr ) ∂ r + ∂σ rz ∂ z − σ θθ r = r ∂ ( r σ zr ) ∂ r + ∂σ zz ∂ z = σ rr = µ ∂ u ∂ r − P (4)2 igure 1: Schematic of droplet domain showing boundary conditions at respective boundary. σ zz = µ ∂ v ∂ z − P (5) σ θθ = µ ur − P (6) σ rz = σ zr = µ (cid:32) ∂ v ∂ r + ∂ u ∂ z (cid:33) (7)The boundary conditions (Figure 1) for the above equation can be expressed as follows: At thebottom of the droplet (at z = u = v =
0) was applied. At r = In our FEM formulation, we have utilized Galerkin weighted residual approach, in which,weight function ( w i ) are the same as shape function or approximation function. Using weightfunction following weak form of Navier-Stokes equation was obtained to lower the requirementof choosing higher order approximation function, which is as follows: (cid:90) Ω e w i (cid:34) r ∂ ( r σ rr ) ∂ r + ∂σ rz ∂ z − σ θθ r (cid:35) d Ω = (cid:90) Ω e w i (cid:34) r ∂ ( r σ zr ) ∂ r + ∂σ zz ∂ z (cid:35) d Ω = d Ω = π rdrdz , di ff erential volume for axisymmetric coordinates and Ω e representsa finite element. After solving and rearranging di ff erential term in above integrals, we get thefollowing equations: (cid:90) Ω e (cid:34) ∂ w i ( σ rr ) ∂ r + ∂ w i σ rz ∂ z + w i σ θθ r (cid:35) rdrdz = (cid:90) Γ e w i ( t r ) rds (10) (cid:90) Ω e (cid:34) ∂ w i ( σ zr ) ∂ r + ∂ w i σ zz ∂ z (cid:35) rdrdz = (cid:90) Γ e w i ( t z ) rds (11)where, t r = σ rr n r + σ rz n z and t z = σ zr n r + σ zz n z are the boundary stress component in radialand axial directions, respectively. Γ e represents the boundary of element and di ff erential surfacearea ( d Γ ) of axisymmetric geometry is 2 π rds . where, ds is di ff erential arc length of the boundary. In penalty function formulation pressure ( P ) is eliminated from Navier-Stokes equations. Theelimination of pressure leads to constrained problem and where constrained equation is the conti-nuity equation (eq. 1). Using penalty function method, constrained problem is reformulated as anunconstrained problem and pressure P is replaced by following function [7]: P = − γ (cid:32) ∂ u ∂ r + ∂ v ∂ z + ur (cid:33) (12)Where, γ is a penalty factor which should be a large arbitrary value [7], i.e. µ × . Note that,After calculating velocities pressure ( P ) can be recovered using above equation (eq. 12). Now, onputting the values of σ rr , σ zz , σ θθ and P from eq. 4, 5, 6, 7 and 12 respectively in eq. 10 and 11,following equation can be obtained: (13) (cid:90) Ω e (cid:34) µ ∂ w i ∂ r ∂ u ∂ r + µ ∂ w i ∂ z (cid:32) ∂ u ∂ z + ∂ v ∂ r (cid:33) + µ w i r ur + γ ∂ w i ∂ r (cid:32) ∂ u ∂ r + ∂ v ∂ z + ur (cid:33) + γ w i r (cid:32) ∂ u ∂ r + ∂ v ∂ z + ur (cid:33)(cid:35) rdrdz = (cid:90) Γ e w i ( t r ) rds (14) (cid:90) Ω e (cid:34) µ ∂ w i ∂ r ∂ u ∂ z + µ ∂ w i ∂ r ∂ v ∂ r + µ ∂ w i ∂ z ∂ v ∂ z + γ ∂ w i ∂ z (cid:32) ∂ u ∂ r + ∂ v ∂ z + ur (cid:33)(cid:35) rdrdz = (cid:90) Γ e w i ( t z ) rds To discretize the above continuous equations within each finite element, the velocities ( u , v )are approximated by the following trial solution: u = n (cid:88) j = ψ j ( r , z ) u j v = n (cid:88) j = ψ j ( r , z ) v j (15)4here u j and v j are the radial and axial velocity at the nodal points of the element. n and ψ j are the number of nodes in the element and shape functions respectively. As explained before, inGalerkin weighted residual approach, weighting functions are same as shape functions, therefore w i = ψ i . On putting trial solution and weighting function, we can write eq. 13 and 14 in compactmatrix notation as follows: (16)(2 µ [ S rr ] + µ [ S zz ] + µ [ M ] + γ ([ S rr ] + [ S ro ] + [ S or ] + [ M ])) { u j } + ( µ [ S zr ] + γ [ S rz ] + γ [ S oz ]) { v j } = { F r } (17)( µ [ S rz ] + γ [ S zr ] + γ [ S zo ]) { u j } + ( µ [ S rr ] + µ [ S zz ] + γ [ S zz ]) { v j } = { F z } Further, we can write above two equations in more simplified way in matrix form as follows:(18) (cid:34) [ K rr ] [ K rz ][ K zr ] [ K zz ] (cid:35) (cid:40) { u j }{ v j } (cid:41) = (cid:40) { F r }{ F z } (cid:41) where, [ K rr ] = µ [ S rr ] + µ [ S zz ] + µ [ M ] + γ ([ S rr ] + [ S ro ] + [ S or ] + [ M ])[ K rz ] = µ [ S zr ] + γ [ S rz ] + γ [ S oz ][ K zr ] = µ [ S rz ] + γ [ S zr ] + γ [ S zo ][ K zz ] = µ [ S rr ] + µ [ S zz ] + γ [ S zz ] (19)where the element coe ffi cient matrices are defined as follows: (20) M i j = (cid:90) Ω e ψ i ψ j r drdzS rri j = (cid:90) Ω e ∂ψ i ∂ r ∂ψ j ∂ r rdrdz S rzi j = (cid:90) Ω e ∂ψ i ∂ r ∂ψ j ∂ z rdrdz (21) S zri j = (cid:90) Ω e ∂ψ i ∂ z ∂ψ j ∂ r rdrdz S zzi j = (cid:90) Ω e ∂ψ i ∂ z ∂ψ j ∂ z rdrdz (22) S roi j = (cid:90) Ω e ∂ψ i ∂ r ψ j drdz S ori j = (cid:90) Ω e ψ i ∂ψ j ∂ r drdz (23) S zoi j = (cid:90) Ω e ∂ψ i ∂ z ψ j drdz S ozi j = (cid:90) Ω e ψ i ∂ψ j ∂ z drdz (24) F ri = (cid:90) Γ e psi i ( t r ) rds F zi = (cid:90) Γ e psi i ( t z ) rds (25)5 .2.2. Velocity-Pressure formulation Velocity-pressure formulation is a natural and direct formulation [7]. In this formulation, theweak form of continuity equation is obtained using weighting function, one order less than usedfor Navier-Stokes equations. Both Navier-Stokes equation physically represents force, hence sameweighting function can be used. However, continuity equation represents the volume change.Volume change occur under the action of hydrostatic pressure, hence weight function ( w k ) forcontinuity, equation should like the pressure ( P ) or the shape function of pressure. The weak formof the continuity equation can be written as follows: (26) − (cid:90) Ω e w k (cid:34) ∂ u ∂ r + ∂ v ∂ z + ur (cid:35) rdrdz = (cid:90) Ω e (cid:34) µ ∂ w i ∂ r ∂ u ∂ r + µ ∂ w i ∂ z (cid:32) ∂ u ∂ z + ∂ v ∂ r (cid:33) + µ w i r ur − ∂ w i ∂ r P − w i r P (cid:35) rdrdz = (cid:90) Γ e w i ( t r ) rds (28) (cid:90) Ω e (cid:34) µ ∂ w i ∂ r ∂ u ∂ z + µ ∂ w i ∂ r ∂ v ∂ r + µ ∂ w i ∂ z ∂ v ∂ z − ∂ w i ∂ z P (cid:35) rdrdz = (cid:90) Γ e w i ( t z ) rds (29) − (cid:90) Ω e (cid:34) w k ∂ u ∂ r + w k ∂ v ∂ z + w k ∂ u ∂ r (cid:35) rdrdz = u , v ) areapproximated by the trial solution defined in eq. 15 and pressure ( P ) is approximated by followingtrial solution, which one order less than used for velocities: P = m (cid:88) l = φ l ( r , z ) P l (30)where P l is the pressure at the corner nodal points of the element. m and φ l are the numberof corner nodes in the element and shape functions respectively. As explained before, in Galerkinweighted residual approach, w i = ψ i and w k = φ l . On putting trial solution and weighting function,we can write eq. 27, 28 and 29 in compact matrix notation as follows: (31)(2 µ [ S rr ] + µ [ S zz ] + µ [ M ]) { u j } + ( µ [ S zr ]) { v j } − ([ G ro ] + [ G oo ]) { P l } = { F r } (32)( µ [ S rz ]) { u j } + ( µ [ S rr ] + µ [ S zz ]) { v j } − ([ G zo ]) { P l } = { F z } (33) − (cid:16) [ G ro ] T + [ G oo ] T (cid:17) { u j } − (cid:16) [ G zo ] T (cid:17) { v j } = [ K rr ] [ K rz ] [ K ro ][ K zr ] [ K zz ] [ K zo ][ K or ] [ K oz ] [0] { u }{ v }{ P } = { F r }{ F z }{ } where, [ K rr ] = µ [ S rr ] + µ [ S zz ] + µ [ M ][ K rz ] = µ [ S zr ][ K ro ] = − ([ G ro ] + [ G oo ])[ K zr ] = µ [ S rz ][ K zz ] = µ [ S rr ] + µ [ S zz ][ K zo ] = − [ G zo ][ K or ] = − (cid:16) [ G ro ] T + [ G oo ] T (cid:17) [ K oz ] = − [ G zo ] T (35)where the [ S ], [ M ] and { F } element coe ffi cient matrices are defined as before from eq. 20 toeq. 25. The remaining element coe ffi cient matrices are defined as follows: (36) G rokl = (cid:90) Ω e ∂ψ k ∂ r φ l rdrdz (37) G zokl = (cid:90) Ω e ∂ψ k ∂ z φ l rdrdz (38) G ookl = (cid:90) Ω e ψ k φ l drdz2.2.3. Evaluation of element coe ffi cient matrices In this model, a six-node triangular element (Figure 2) was chosen over three-node elementas it improves the solution accuracy and reduce the overall solution time by reducing the numberelement required. Quadratic shape functions were used to evaluate velocity field and liner shapefunctions were used to evaluate pressure field. Six-node triangular element contains a node on themid-side of each of the three side of the triangle, which defines the boundary of the element.In order to evaluate the integrals of element coe ffi cient matrices, numerical integration wasemployed. To do numerical integration, first we need to map real element in ( r , z ) coordinatesto a master or parent element in generic ( ξ , η ) coordinates over which integration will performedusing Gaussian quadrature points. As in our real element all three nodes located at side of triangleare located at mid-point of each triangle side, so it is su ffi cient to use linear shape function to7 igure 2: Six-node triangular element. do coordinate transformation. The following expression facilitates the ( r , z ) to ( ξ , η )coordinatetransformation: r = (cid:88) k = r k φ k ( ξ, η ) z = (cid:88) k = z k φ k ( ξ, η ) (39)where r k and z k is k th radial and axial coordinates of corner nodes of the triangle element. Shapefunctions φ k can be defined as follows: φ = ξφ = ηφ = − η − ξ (40)To evaluate integrals, the following expression is required to map infinitesimal area of realelement to corresponding area of parent element: drdz = det [ J ] d ξ d η (41)where, det [ J ] is the determinant of the Jacobian ([ J ]) for the element in consideration and isdefined as follows: (42)[ J ] = (cid:34) ∂ r ∂ξ ∂ z ∂ξ∂ r ∂η ∂ z ∂η (cid:35) det [ J ] is equal to 2 times area of triangle in consideration. Next, we needto convert partial derivative of ( r , z ) into derivatives of ( ξ , η ) which can be done using followingrelationship: (43) (cid:40) ∂ψ i ∂ r ∂ψ i ∂ z (cid:41) = [ J ] − ∂ψ i ∂ξ∂ψ i ∂η where, shape function ψ i is quadratic shape function for a six-node triangle element as men-tioned before. Quadratic shape can be defined as follows: ψ = ξ − ξψ = η − ηψ = − ξ − η + ξη + ξ + η ψ = ξηψ = η − ξη − η ) ψ = ξ − ξη − ξ ) (44)Using eq. 41, 43 and 44, the integrals specified in equations from eq. 20 to 24 can be convertedto ( ξ , η ) coordinates and numerical integration can be performed using Gaussian quadrature asfollows: (cid:90) Ω e f ( ξ, η ) d ξ d η ∼ = gp (cid:88) l = w l f ( ξ l , η l ) (45)where f ( ξ, η ) is any function in ( ξ , η ) coordinates and ( ξ l , η l ) are the Gaussian points whichare well document in any FEM textbook as given in [7]. gp is the number of Gaussian points usedto evaluate the integral. Three and one-point Gaussian quadrature was used to evaluate integrals.To calculate coe ffi cient of force vector on the right hand side of eq. 18 or to apply stressboundary condition on the boundary of the domain, we need to transform the boundary or openside of the triangle to one dimensional line coordinate ( s ) whose origin at first point of boundary.As we are using six-node or quadratic triangular element, we have three nodes on each side of thetriangle, so we need to use quadratic shape function ( ψ bi ). These three shape function ( ψ bi ) can bedefined as follows: ψ b = (cid:18) − sh (cid:19) (cid:32) − sh (cid:33) ψ b = sh (cid:18) − sh (cid:19) ψ b = − sh (cid:32) − sh (cid:33) (46)where, h is the length of the side or boundary of the element. The elements which are inside thedomain, for them, the value of coe ffi cients in force vector are need not to be calculated as the forceacting on the side or boundary of the neighboring elements cancels out each other. Therefore, we9eed to calculate coe ffi cient of force vector only at the boundary of elements that are located atboundary of the domain and some external force or stress is applied. At the bottom of the droplet (at z = u = , v =
0) was applied,which is a Dirichlet boundary condition. At r =
0, axisymmetric boundary condition was applied,which is u = τ t ) and kinematicboundary condition in normal direction ( u n ). The liquid-gas interface boundary condition can beapplied using two di ff erent approach and both were tested. These approaches are described asfollows: Direct approach
As the shape of droplet is spherical cap, the liquid-gas interface is not parallel to r or z axis.To apply the both boundary conditions precisely, coordinate rotation must be done for the nodes atthe liquid-gas interface of the droplet [7], which can be done with the help of rotation matrix [ Q ].Eq. 18 can be written in more compact form as follows:[ K ] { u } = { F } (47)To perform the coordinate rotation on the nodes at the liquid-gas interface of the droplet fol-lowing mathematical operation can be done:[ Q ] T [ K ] [ Q ] { u } = { F } (48)where, rotation matrix [Q] can be defined as follows: (49)[ Q ] = . . . · · · · · · · · ·· · · n r , i n z , i · · ·· · · − n z , i n r , i · · ·· · · · · · · · · . . . where, n r , i and n z , i are the normal vector in radial and axial direction for i th node on the liquid-gas interface. n r and n z can be calculated as follows: n r = r sin θ Rn z = (cid:112) − ( n r ) (50)where r is the radial distance of the considered node and θ is the contact angle of the droplet.In the normal direction of the liquid-gas interface, kinematic boundary condition can be definedas follows [3]: u n = j ρ − n z ∂ h ∂ t (51)10here u n is the normal velocity at the liquid-gas interface, j /ρ is outflow velocity of dropletevaporation and ∂ h /∂ t is the surface motion of spherical cap to retain its spherical shape. ∂ h /∂ t can be defined as follows [3]: ∂ h ∂ t = ˙ h d h d − R h d (cid:115)(cid:32) h d + R h d (cid:33) − r + R − h d h d (52)where ˙ h d is the velocity of the top of the evaporating droplet, which can be calculated as follows[2] ˙ h d = m ρπ ( R + h d ) (53)In the tangential direction, we apply zero shear stress boundary condition (i.e. τ t = t z = Penalty approach
In this approach, kinetic boundary condition is treated as constrained equation, which can bewritten as follows: u r n r + u z n z = u n (54)where u r and u z are the velocity component in r and z directing at the liquid-gas interface.Using penalty function method governing equations were solved and the constraint equation (eq.54) is made to satisfy, which yields solution in following form [7]:([ K ] + [ K P ]) { u } = { F } + { F P } (55)where, (56)[ K P ] = · · · · · · · · · · · ·· · · γ p n r γ p n r n z · · ·· · · γ p n r n z γ p n z · · ·· · · · · · · · · · · · (57) { F P } = · · · γ p u n n r γ p u n n z · · · Where, γ is a penalty factor which should be a large arbitrary value [7], i.e. max ([ K ]) × .Thus, a modification of coe ffi cient in [ K ] and { F } matrix associated with boundary nodes will yieldthe desired solution with the constraint kinetic boundary condition.Zero shear stress boundary condition can be applied by dividing the tangential shear stress into r and z components, i.e. t r = τ t n r = t z = τ t n z = (m) × -3 z ( m ) × -4 Figure 3: Six-node mesh in droplet domain. Red dots represents nodes. .2.5. Mesh generation As the whole code in written in MATLAB, we were more inclined to use MATLAB functionsto save ourselves from any compatibility and data import issues. To generate six-node triangularelement mesh, first we generated a three-node triangular mesh inside the droplet domain usingMATLAB function initmesh. To the best of author knowledge, there is no MATLAB functionavailable to generate six-node triangular mesh. Therefore, to generate six-node mesh, we devel-oped an in-house code, which places a node at the exact midpoint of the each side of triangularelement and updates the point matrix, element connectivity matrix, boundary matrix. The finalmesh generated using in-house code and MATLAB function ‘initmesh is shown in Figure 3.
3. Results and Validation of code
We have used two formulations (Penalty function formulation and Velocity-pressure formu-lation) and two approaches (Direct and Penalty approach) of boundary condition application tosolve flow filed inside the droplet, which makes four combinations of methods to solve flow filed.Here, we are presenting validation using direct approach of boundary condition application forboth formulations to avoid repeatability of the similar results of penalty approach application ofboundary condition. ff erent formulations To validate the code radial velocity at di ff erent radial location along the height of the droplethas been compared with results of model of Hu and Larson [3]. Wetting angle and wetting radiusof droplet considered are 40 ◦ and 1 mm. Figure 4 shows the radial velocity ( u ) profiles calculatedalong axial direction ( z ) at di ff erent radial locations, r = Fluid velocity inside the droplet increases as flow moves from droplet interior to the three-phase contact line. This phenomenon is evident from three velocity profiles shown in Figure 4.At radial location r = . × − mm / s and as fluid move radiallyoutward at location r = . × − mm / s and1 . × − mm / s, respectively. However, in axial direction, from z = (mm/s) z ( mm ) Hu & Larson, r = 0.1Hu & Larson, r = 0.6Hu & Larson, r = 0.9VPF, r = 0.1VPF, r = 0.6VPF, r = 0.9PF, r = 0.1PF, r = 0.6PF, r = 0.9
Figure 4: Validation of Velocity-pressure and Penalty function formulation: comparison of radial velocity profiles (atdi ff erent radial locations) obtained from present model with the model of Hu and Larson [3]. (mm) z ( mm ) -03 -03 -03 -03 -03 -03 -04 -04 -04 +00 u (mm/s) Figure 5: Streamlines and u -radial velocity contour inside the evaporating droplet. r (m) × -3 z ( m ) × -4 ( m / s ) × -6 -0.500.511.522.533.54 Figure 6: Axial velocity ( v ) field inside the evaporating droplet. (m) × -3 z ( m ) × -4 ( N / m ) × -3 -4.5-4-3.5-3-2.5-2-1.5-1-0.50 Figure 7: Pressure field ( P ) inside the evaporating droplet. the interior of the droplet rushes towards three-phase contact line to maintain droplet spherical capshape. Streamline shown in Figure 5 shows the same behaviour and characteristics and corroboratethe above explanation regarding fluid flow inside the droplet.Figure 6 shows the axial velocity ( v ) contour inside the evaporating droplet. It can be con-cluded from the figure that magnitude of v increases away from the droplet-air interface to droplet-substrate contact. Figure 7 shows the pressure ( P ) contour inside the droplet. it can be observedthat pressure inside the droplet doesn’t vary much, it remains almost same except at the contactline.
4. Conclusions
Present work explains the finite element modeling of fluid flow inside evaporating sessiledroplet resting on hydrophilic substrate. In this modeling, constant contact radius (CCR) mode ofevaporation was considered. Finite element code was implement in MATLAB. Galerkin weightedresidual approach was used to formulate weak form of the numerical equations. Velocity-pressureformulation was applied to discretise the Navier-Stokes equation in cylindrical coordinate system.Six-node triangular mesh was used in the simulations for higher accuracy of the solution. Resultsobtained from the simulations were compared with the model of Hu and Larson [3] and goodagreement was found in present work and the model of Hu and Larson.16 eferences [1] H. Ghasemi, C. Ward, Energy transport by thermocapillary convection during sessile-water-droplet evaporation,Physical review letters 105 (13) (2010) 136102.[2] M. R. Barmi, C. D. Meinhart, Convective flows in evaporating sessile droplets, The Journal of Physical ChemistryB 118 (9) (2014) 2414–2421.[3] H. Hu, R. G. Larson, Analysis of the microfluid flow in an evaporating sessile droplet, Langmuir 21 (9) (2005)3963–3971.[4] Y. O. Popov, Evaporative deposition patterns: spatial dimensions of the deposit, Physical Review E 71 (3) (2005)036313.[5] R. D. Deegan, Pattern formation in drying drops, Physical review E 61 (1) (2000) 475.[6] Y. Li, C. Lv, Z. Li, D. Qu´er´e, Q. Zheng, From co ff ee rings to co ff ee eyes, Soft Matter 11 (23) (2015) 4669–4673.[7] J. Reddy, An Introduction to the Finite Element Method, 3rd Edition, Tata Mc Graw Hill, New Delhi, 2005.[8] H. Hu, R. G. Larson, Evaporation of a sessile droplet on a substrate, The Journal of Physical Chemistry B 106 (6)(2002) 1334–1344.ee eyes, Soft Matter 11 (23) (2015) 4669–4673.[7] J. Reddy, An Introduction to the Finite Element Method, 3rd Edition, Tata Mc Graw Hill, New Delhi, 2005.[8] H. Hu, R. G. Larson, Evaporation of a sessile droplet on a substrate, The Journal of Physical Chemistry B 106 (6)(2002) 1334–1344.