# An AMR-Based Liquid Film Simulation with Surfactant Transport Using PLIC-HF Method
## 3 Numerical Methods
### Finite Difference Method to Solve Fluid Calculations  
To solve Equations (1) and (2), undefined velocities at cell centers and cell corners are linearly interpolated from the nearest points:  
\[u_{i,j}=\left(u_{i-\frac{1}{2},j}+u_{i+\frac{1}{2},j}\right)/2 \tag{29}\]  
\[u_{i-\frac{1}{2},j+\frac{1}{2}}=\left(u_{i-\frac{1}{2},j}+u_{i-\frac{1}{2},j+1 }\right)/2 \tag{30}\]  
\[v_{i-\frac{1}{2},j+\frac{1}{2}}=\left(v_{i-1,j+\frac{1}{2}}+v_{i,j+\frac{1}{2 }}\right)/2 \tag{31}\]  
Undefined scalar values at cell faces and cell corners are interpolated in the same way. For example, the density at the cell face \(i,j-\frac{1}{2}\) is determined by linear interpolation, as follows:  
\[\rho_{i,j-\frac{1}{2}}=\frac{1}{2}\left(\rho_{i,j}+\rho_{i,j-1}\right) \tag{32}\]  
In the above discretization, which considers the spatial change of density, numerical instability is likely to occur because the differential operation is performed near the interface,  
Figure 3: Configuration of the staggered grid system. Vector components are defined at cell faces, and scalars are defined at cell centers. where physical quantities such as density undergo drastic changes. In order to improve the numerical stability, the gradient restriction function is applied to the gradient obtained by the third-order WENO (Weighted Essentially Non-Oscillatory) scheme [41]:  
\[f_{x}=f_{x}^{1 st}+\Psi(r)\left(f_{x}^{WENO}-f_{x}^{1 st}\right) \tag{33}\]  
where \(f\) is a variable and \(r\) is a smooth indicator that is expressed as  
\[r_{i-\frac{1}{2}}= \begin{cases}\frac{f_{i+\frac{1}{2}}-f_{i-\frac{1}{2}}}{f_{ i-\frac{1}{2}}-f_{i-\frac{3}{2}}}&u_{i-\frac{1}{2}}>&0\\ \frac{f_{i+\frac{1}{2}}-f_{i-\frac{1}{2}}}{f_{i-\frac{1}{2}}-f_{i-\frac{3}{2}} }&u_{i-\frac{1}{2}}\leq&0\end{cases}  \tag{34}\]  
The minmod limiter function [42] was used for the restriction function \(\Psi\). The derivative value \(\left(\frac{\partial f}{\partial x}\right)_{i}\) of a certain grid point \(i\) can be obtained as follows, assuming that the lattice width is \(\Delta x\). In the third-order precision HJ-WENO scheme, the sum of the central difference and the numerical viscosity can be expressed as follows:  
\[\left(\frac{\partial f}{\partial x}\right)_{i}=\frac{f_{i+1}-f_{i-1}}{2\Delta x }+\frac{f_{x}^{\pm}}{1+2r^{2}} \tag{35}\]  
The numerical viscosity was constructed by taking a 4-point stencil depending on the upwind direction with respect to the advection velocity \(u\), as follows:  
\[\left\{ \begin{array}{l}\Delta t=\Delta x\\ f_{pp}=f_{i+2},f_{p}=f_{i+1},\quad(u\leq 0)\\ f_{m}=f_{i-1}\end{array} \right. \tag{36}\]  
\[\left\{ \begin{array}{l}\Delta t=-\Delta x\\ f_{pp}=f_{i-2},f_{p}=f_{i-1},\quad(u>0)\\ f_{m}=f_{i+1}\end{array} \right. \tag{37}\]  
\[f_{x}^{\pm}=\frac{-f_{pp}+3f_{p}-3f_{i}+f_{m}}{2\Delta h}, \tag{38}\]  
\[r=\frac{f_{i}-2f_{p}+f_{pp}+\epsilon}{f_{p}-2f_{i}+f_{m}+\epsilon}. \tag{39}\]  
In Equation (39), \(\epsilon\) is a very small constant value that prevents the denominator from being zero. \(r\) is an index of the smoothness of the function, and when \(r=1\), the scheme shown above becomes the same as the third-order upwind difference method. The pressure gradient \(\
abla p\) included in the time evolution equation of velocity is calculated by the second-order central difference method at the cell face \(i+\frac{1}{2}\), as follows:  
\[\left(\frac{\partial p}{\partial x}\right)_{i+\frac{1}{2}}=\frac{p_{i+1}-p_{ i}}{\Delta x} \tag{40}\]  
The viscous stress term \(\
abla\cdot\mathbf{\tau}=\
abla\cdot\left\{\mu\left(\
abla\mathbf{u}+\
abla\mathbf{u}^{T} \right)\right\}\) is also calculated by the second-order central difference method. \(\mu\) is the spacial difference when the viscous stress term is discretized.