# An AMR-Based Liquid Film Simulation with Surfactant Transport Using PLIC-HF Method
## 3 Numerical Methods
### Finite Difference Method to Solve Fluid Calculations
In the components of the viscous stress tensor, if the directions of the velocity and the spacial difference are the same, such as in the term \(\frac{\partial}{\partial x}\left(\mu\frac{\partial u}{\partial x}\right)\), the discretization can be obtained by  
\[\left[\frac{\partial}{\partial x}\left(\mu\frac{\partial u}{\partial x}\right) \right]_{i-\frac{1}{2},j}=\frac{\left(\mu\frac{\partial u}{\partial x}\right) _{i,j}-\left(\mu\frac{\partial u}{\partial x}\right)_{i-1,j}}{\Delta x}. \tag{41}\]The stencil calculation contained in Equation (41) is shown in Figure 4a. In the term \(\frac{\partial}{\partial y}\left(\mu\left(\frac{\partial u}{\partial y}+\frac{ \partial v}{\partial x}\right)\right)_{i-\frac{1}{2},j}\), which is included in the time evolution equation of \(u\), the directions of the velocity and spacial difference are different, and discretization can be achieved by  
\[\left[\frac{\partial}{\partial y}\left(\mu\left(\frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}\right)\right)\right]_{i-\frac{1}{2},j}=\frac{ \left[\mu\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)\right]_{i-\frac{1}{2},j+\frac{1}{2}}-\left[\mu\left(\frac{\partial u} {\partial y}+\frac{\partial v}{\partial x}\right)\right]_{i-\frac{1}{2},j- \frac{1}{2}}}{\Delta y}. \tag{42}\]  
Each term in Equation (42) can be obtained with the stencil calculation, as shown in Figure 4b. The velocity divergence term included on the right side of the pressure evolution equation, Equation (2), is evaluated by the second-order central difference, as follows:  
\[\left(\
abla\cdot\mathbf{u}\right)_{i,j}=\frac{u_{i+\frac{1}{2},j}-u_{i-\frac{1}{2 },j}}{\Delta x}+\frac{v_{i,j+\frac{1}{2}}-v_{i,j-\frac{1}{2}}}{\Delta y} \tag{43}\]  
In Equation (27), the surface tension is defined at the cell faces, so the surface tension coefficient \(\sigma(f)\) and curvature \(\kappa\) need to be interpolated from the cell center values:  
\[\sigma(f_{i+\frac{1}{2}})=\frac{\sigma(f_{i})+\sigma(f_{i+1})}{2} \tag{44}\]  
\[\kappa_{i+\frac{1}{2}}=\frac{\kappa_{i}+\kappa_{i+1}}{2} \tag{45}\]  
Figure 4: Stencil computations of the viscous term where **(a)** the directions of the derivative and velocity are the same and **(b)** the directions of the derivative and the velocity are different. **(a)** Stencil computation where the directions of the derivative and the velocity are the same. **(b)** Stencil computation where the directions of the derivative and the velocity are different. The pressure Poisson equation is solved in a few initial time steps to suppress acoustic waves. It is spatially discretized as  
\[\frac{\partial p_{i,j}}{\partial t}=-\rho_{i,j}c_{s}^{2}\bigg{(}\frac{u_{i+\frac{1 }{2},j}-u_{i-\frac{1}{2},j}}{\Delta x}+\frac{v_{i,j+\frac{1}{2}}-v_{i,j-\frac{1 }{2}}}{\Delta y}\bigg{)} \tag{46}\]  
and solved with the parallel-computing-suitable red/black successive over relaxation (SOR) method [43]. The advection and re-initialization equation of the level set function is solved with the third-order HJ-WENO scheme [20].