In this paper, the MRT multi-relaxation factor model and the large density ratio pseudopotential model are mainly used to simulate the process and motion characteristics of water droplets hitting stationary walls. The collision evolution of the basic MRT multi-relaxation model can be expressed as follows:
$${f}_{i}\left(x+{e}_{\alpha }{\delta }_{t},t+{\delta }_{t}\right)={f}_{i}\left(x,t\right)-{\Lambda }\left({f}_{i}\left(x,t\right)-{f}_{i}^{eq}\left(x,t\right)\right)+{\delta }_{t}(I-0.5\stackrel{-}{{\Lambda }})S$$
1
where \({f}_{i}\) is the distribution function, \({f}_{i}^{eq}\) is the equilibrium distribution function, \(\stackrel{-}{{\Lambda }}\) is the collision diagonal matrix, \(I\) is the identity matrix, and \(S\) is the external force term. Through orthogonal matrix transformation, the distribution function is converted from velocity space to matrix space, including equilibrium distribution function, external force term, etc: \(\text{m}=\text{M}\bullet \text{f},{m}^{eq}=M\bullet {f}^{eq},\stackrel{-}{{\Lambda }}=M{\Lambda }{M}^{-1}.\)For the D2Q9 model \({m}^{eq}={\rho (1,-2+3{\left|{\Lambda }\right|}^{2},1-3{\left|{\Lambda }\right|}^{2},{v}_{x},-{v}_{x},{v}_{y},-{v}_{y},{v}_{x}^{2}-{v}_{y}^{2},{v}_{x}{v}_{y})}^{T}\), where \(v\) is the equilibrium velocity. The collision evolution equation is transformed by an orthogonal matrix as follows
$${m}^{*}=m-\stackrel{-}{{\Lambda }}\left(m-{m}^{eq}\right)+{\delta }_{t}(I-\frac{\stackrel{-}{{\Lambda }}}{2})\stackrel{-}{S}$$
2
where \({m}^{*}\) is the evolved distribution function and \(\stackrel{-}{S}\) is the force term in matrix space. The collision diagonal matrix in the D2Q9 model is: \(\stackrel{-}{{\Lambda }}=diag({\tau }_{\rho }^{-1},{\tau }_{e}^{-1},{\tau }_{\varsigma }^{-1},{\tau }_{j}^{-1},{\tau }_{q}^{-1},{\tau }_{j}^{-1},{\tau }_{q}^{-1},{\tau }_{\upsilon }^{-1},{\tau }_{\upsilon }^{-1})\),\({\tau }_{\rho }^{-1},{\tau }_{e}^{-1},{\tau }_{\varsigma }^{-1},{\tau }_{j}^{-1},{\tau }_{q}^{-1},{\tau }_{\upsilon }^{-1}\) in the equation are the multiple relaxation factors of the matrix space, and reasonable settings are required to achieve the stability of the model. The pseudopotential force format is used as follows:
$$\stackrel{-}{S}={\left[\text{0,6}\left({u}_{x}{F}_{x}+{u}_{y}{F}_{y}\right)+\frac{12\varpi {\left|{F}_{m}\right|}^{2}}{{\psi }^{2}{\delta }_{t}({\tau }_{e}-0.5)},-6\left({u}_{x}{F}_{x}+{u}_{y}{F}_{y}\right)-\frac{12\varpi {\left|{F}_{m}\right|}^{2}}{{\psi }^{2}{\delta }_{t}\left({\tau }_{\varsigma }-0.5\right)},{F}_{x},{F}_{y},{F}_{x},{F}_{y},2\left({u}_{x}{F}_{x}-{u}_{y}{F}_{y}\right),\left({u}_{x}{F}_{y}-{u}_{y}{F}_{x}\right)\right]}^{T}$$
3
In the pseudopotential force format, the velocity is \(u=v+\varpi F/\left(\right(\tau -0.5\left){\psi }^{2}\right)\), the interparticle force is \({F}_{m}=-G\psi (x,t)\left(\sum _{\alpha }w\left({\left|{e}_{\alpha }\right|}^{2}\right)\psi (x+{e}_{\alpha },t){e}_{\alpha }\right)\), the component of the total force is \(F=\left({F}_{x}{,F}_{y}\right)\), where the \(\varpi\) is 0.114 according to the mechanical stability. The surface tension in the pseudopotential model can be adjusted by constructing an external force term in the collision evolution equation as follows:
$${m}^{*}=m-\stackrel{-}{{\Lambda }}\left(m-{m}^{eq}\right)+{\delta }_{t}\left(I-\frac{\stackrel{-}{{\Lambda }}}{2}\right)\stackrel{-}{S}+{\delta }_{t}C$$
4
In the formula:\(\text{C}={\left[\text{0,1.5}{\tau }_{e}^{-1}\left({Q}_{xx}+{Q}_{yy}\right),-1.5{\tau }_{\varsigma }^{-1}\left({Q}_{xx}+{Q}_{yy}\right),\text{0,0},\text{0,0},-{\tau }_{\upsilon }^{-1}\left({Q}_{xx}-{Q}_{yy}\right),-{\tau }_{\upsilon }^{-1}{Q}_{xy}\right]}^{T}\), where\({Q}_{xx},{Q}_{yy},{Q}_{xy}\), can be calculated by the following formula:
$$\text{Q}=\text{k}\frac{G}{2}\psi (\text{x},\text{t})\left[\sum _{i=1}^{8}w\left({\left|{e}_{i}\right|}^{2}\right)(\psi \left(x+{e}_{i},\text{t}\right)-\psi \left(x,\text{t}\right)){e}_{i}{e}_{i}\right]$$
5
The value of the surface tension coefficient k is 0.35. Due to mechanical instability and thermodynamic inconsistency: \({\int }_{{\rho }_{g}}^{{\rho }_{l}}\left({P}_{0}-\rho {C}_{s}^{2}-\frac{G{c}^{2}}{2}{\psi }^{2}\right)\frac{{\psi }^{\prime }}{{\psi }^{1+\epsilon }}d\rho =0\),\({\int }_{{\rho }_{g}}^{{\rho }_{l}}\left({P}_{0}-{P}_{EOS}\right)\frac{{\psi }^{\prime }}{{\psi }^{1+\epsilon }}d\rho =0\). Therefore, the pseudopotential function adopts the radical model: \({\int }_{{\rho }_{g}}^{{\rho }_{l}}\left({P}_{0}-{P}_{EOS}\right)\frac{{\psi }^{\prime }}{{\psi }^{1+\epsilon }}d\rho =0\), and the gas equation of state adopts the C-S equation of state: \({P}_{EOS}=\rho RT\frac{1+b\rho /4+{\left(b\rho /4\right)}^{2}-{\left(b\rho /4\right)}^{3}}{{(1-b\rho /4)}^{3}}-a{\rho }^{2}\), critical temperature and critical pressure: \({T}_{c}=0.3773a/bR\), where \(a\)=0.25, \(b\)=1, \(R\)=1, temperature ratio is 0.5. The wall wettability model uses:
$${\rho }_{w}\left(x\right)=\left\{\begin{array}{c}\phi {\rho }_{avg} \phi \ge 1,for decreasing \theta \\ {\rho }_{avg}-△\rho △\rho \ge 0,for increasing \theta \end{array}\right.$$
6
In the formula, \({\rho }_{avg},\phi ,△\rho\) are the average density, hydrophilic coefficient and hydrophobic coefficient of the fluid at the wall, where \({\rho }_{avg}=\frac{{\sum }_{\alpha }{\omega }_{\alpha }\rho (x+{e}_{\alpha }{\delta }_{t}){S}_{w}(x+{e}_{\alpha }{\delta }_{t})}{{\sum }_{\alpha }{\omega }_{\alpha }{S}_{w}(x+{e}_{\alpha }{\delta }_{t})}\), \({S}_{w}(x+{e}_{\alpha }{\delta }_{t})\) is the judgment function, 1 for fluids and 0 for solids. The following equation is used to calculate the macrodensity, macroscopic velocity and equilibrium velocity: \(\rho =\sum {f}_{\alpha }\text{、}\rho v=\sum {e}_{\alpha }{f}_{\alpha }+\frac{{\delta }_{t}F}{2}\).