# Low Reynolds Number Swimming Near Interfaces in Multi-Fluid Media
## 2 Model Equations and Numerical Method
### Immersed Boundary Method for Multi-Fluid Mixture  
Immersed Boundary (IB) Method is a powerful computational approach for handling the dynamic interactions between fluid and immersed elastic structures [17]. The simplicity and robustness of the IB method lead to its numerous applications in biological problems [18; 19; 20; 21; 22]. In [13], a two-phase IB scheme is developed for the simulation of interactions between elastic structures and a two-fluid mixture. Here we briefly describe the basic elements of the method. In the two-phase IB scheme, an Eulerian description is used for the fluid variables such as the velocity and pressure. A Lagrangian description is used for the immersed structures. Each immersed structure is represented by two immersed boundaries, denoted by \(\Gamma^{n}\) (network IB) and \(\Gamma^{s}\) (solvent IB), respectively. As shown in Figure 1, the spatial locations of IB points on \(\Gamma^{n}\) and \(\Gamma^{s}\) are represented by the vector functions \(\mathbf{X}_{n}(q,t)\) and \(\mathbf{X}_{s}(q,t)\), respectively, with \(\mathbf{X}_{n}(q,0)=\mathbf{X}_{s}(q,0)\). Here q is the parametrization variable and \(\mathbf{X}_{n}(q,0)=\mathbf{X}_{s}(q,0)\). In this work, both the solvent and the network satisfy non-slip condition on the immersed structure. That is, the IB points on \(\Gamma^{n}\) move with the local network velocity \(\mathbf{u}_{n}\). The IB points on \(\Gamma^{s}\) move with the local solvent velocity \(\mathbf{u}_{s}\). To prevent \(\Gamma^{n}\) and \(\Gamma^{s}\) from separating, the pair of corresponding IB points on them is connected by a stiff spring with zero rest length. Penalty forces equal in magnitude with opposite directions are generated at these two points when their spatial locations are different. The penalty force density at \(\mathbf{X}_{n}\) is computed by \(\mathbf{F}_{n}^{P}=k_{p}(\mathbf{X}_{s}-\mathbf{X}_{n})\), where \(k_{p}\) is the penalty spring constant. The penalty force density \(\mathbf{F}_{s}^{P}\) at \(\mathbf{X}_{s}\) is computed in a similar way. In addition to the penalty forces, other types of forces are exerted at each IB point to simulate locomotion problems in this paper. First, to specify the internal elastic property of the swimmer, each IB point is connected by elastic springs to its two neighboring points on the same immersed boundary. Second, each IB point is connected by a stiff spring to a corresponding \"tether\" point with imposed track of motion. The movements of the tether point describe the specific gait of swimming. In addition to the penalty force densities \(\mathbf{F}_{n}^{P}\) and \(\mathbf{F}_{s}^{P}\), we denote the densities of other force components (sum of elastic force and tether force) on the two boundaries by \(\mathbf{F}_{n}^{o}\) and \(\mathbf{F}_{s}^{o}\). Following the classical IB method, the coupling between the fluids and the immersed boundaries is through the integral relations:  
\[\mathbf{f}_{j}^{i}(\mathbf{x},t)=\int_{\Gamma^{j}}\mathbf{F}_{j}^{i}(q,t) \delta\big{(}\mathbf{x}-\mathbf{X}_{j}(q,t)\big{)}dq \tag{8}\]  
\[\frac{\partial\mathbf{X}_{j}(q,t)}{\partial t}=\int_{\Omega}\mathbf{u}_{j}( \mathbf{x},t)\delta\big{(}\mathbf{x}-\mathbf{X}_{j}(q,t)\big{)}dx \tag{9}\]Here \(i=o,p\) and \(j=n,s\). \(\Omega\) is the fluid domain. \(\mathbf{f}_{j}^{i}(\mathbf{x},t)\) represents the Eulerian force density at \(\mathbf{x}\) contributed by the Lagrangian force \(\mathbf{F}_{j}^{i}\). \(\delta(\mathbf{x})=\delta(\mathbf{x})\delta(\mathbf{y})\) is the two-dimensional Dirac delta function. Equation (8) describes how the Lagrangian forces at each IB point are transmitted to the corresponding fluid phase. Equation (9) reflects the fact that each IB object moves with the local fluid velocity.