## 3 Immersed Boundary Method
### Diffuse Interface IBM
#### 3.1.2 Penalty Forcing IBM  
The classic IBM computes the boundary force \(\mathbf{F}(s,t)\) by Equation (9) for elastic structures, where the function \(S(\mathbf{X}(s,t),t)\) is determined by a certain problem. The feedback scheme proposed by [PERSON] et al. [11] established a unified model for boundary force, which inspired a new IBM, the so-called penalty forcing IBM. Dimensional constants \(\alpha\) and \(\beta\) are introduced in the feedback loop:  
\[F(s,t)=\alpha\int_{0}^{t}\mathbf{u}(s,t^{\prime})dt^{\prime}+\beta\mathbf{u}(s,t), \tag{11}\]  
where \(\mathbf{u}\) is the fluid velocity at the boundary surface. To satisfy the no-slip boundary condition, the constants \(\alpha\) and \(\beta\) are set negative and large enough that then \(\mathbf{u}\) will approximately converge to zero for stationary objects. Inspired by this, the penalty forcing IBM has been widely studied and applied from then on [12; 13; 15]. However, it is difficult to select appropriate values with physical meaning for the two constants. For different flow problems, the optimal values may be different. Therefore, several penalty forcing IBMs are developed using particular cases of Equation (11). By choosing \(\alpha=0\) and \(\beta=\mu/K\)[72; 73], the external force becomes \(\mathbf{F}=\mu\mathbf{u}/K\). This treatment is suitable for flows with a porous medium with permeability \(K\). So far, the porous medium is rarely introduced in biological motion simulations, but it may be useful in more detailed simulations of the feather or fin in the future. By setting \(\beta=0\)[13], the boundary force can be rewritten as:  
\[\mathbf{F}(s,t)=-\kappa(\mathbf{X}(s,t)-\mathbf{X}^{\epsilon}(s,t)), \tag{12}\]  
where the position of boundary points \(\mathbf{X}(s,t)\) is connected to the corresponding equilibrium points \(\mathbf{X}^{\epsilon}(s,t)\) with a spring stiffness constant \(\kappa\). For the rigid boundary, \(\mathbf{X}^{\epsilon}(s,t)\) is the desired position, and the spring model is defined with large stiffness \(\kappa\gg 1\). The movement of the boundary points can be constrained to some degree. A famous application of the penalty forcing IBM for rigid bodies is the simulations of particle-laden flows. [PERSON] and [PERSON] [14] define a displacement \(\xi(s,t)=\mathbf{X}(s,t)-\mathbf{X}^{\prime}(s,t)\) for each Lagrangian point, where \(\mathbf{X}^{\prime}\) is the reference boundary point position that can be computed using the position of mass center \(\mathbf{X}_{c}\) and the instantaneous rigid rotational matrix \(R\):  
\[\mathbf{X}^{\prime}(s,t)=\mathbf{X}_{c}(t)+R(t)(\mathbf{X}^{\prime}(s,0)-\mathbf{X}_{c}(0)). \tag{13}\]  
The restoration force \(F(s,t)\) is needed to restore the boundary point back to its prescribed position from \(\mathbf{X}\) to \(\mathbf{X}^{\prime}\). The expressions of \(F(s,t)\) is given by  
\[F(s,t)=-\kappa\xi(s,t), \tag{14}\]  
Figure 3: Four typical kernel functions.  
which is a linear spring relationship and similar to Equation (12). Another piecewise expression is written as  
\[\mathbf{F}(s,t)=\begin{cases}0,&\|\xi(s,t)\|=0\\ -\frac{c_{f}}{e_{d}}\frac{\xi(s,t)}{\|\xi(s,t)\|},&\|\xi(s,t)\|>0^{\prime}\end{cases} \tag{15}\]  
where the coefficients \(c_{f}\) and \(e_{d}\) are scale factors, corresponding to force and stiffness, respectively.  
Compared with the original IBM, this spring-like model for the solid boundary can behave well in simulating elastic structures [13; 23] and preliminarily apply to rigid body simulations. However, the key issue for the rigid boundary is choosing the value of stiffness \(\kappa\). Although a large value can maintain the rigidity of an immersed body, it leads to severe time-step limitations [13]. Conversely, a small value mitigates the rigidity of the boundary by the spurious elastic effects.