4.3.1 Discrete skew-adjoint operators
In this subsection, we show that the spreading operator S acting on the elastic tension and the surface divergence operator ∇s acting on the velocity are also skew-adjoint in the discrete sense. That is, we shall prove the numerical identity for Eq. (4.8). To proceed, we first define the corresponding discrete inner product on the fluid grid Ωh and the interfacial grid Γh in the following
hu, viΩh =X
where the second summation is nothing but the mid-point rule for the second integral of (4.7). We also define the discrete spreading operator Sh acting on the discrete elastic tension σh as
Sh(σh) =
One should notice that this discrete skew-adjoint property is crucial to our IB formu-lation for solving Eqs. (4.9) to (4.11). Due to the fact that discrete surface divergence of the velocity is zero in Eq. (4.11), we can re-scale this constraint to make the resul-tant matrix obtained by Eq. (4.11) is the transpose of the resulresul-tant matrix obtained by the discrete spreading operator of the tension. One can also verify this symmetric property by expressing the terms explicitly. The detail is given in the Appendix B.
We now are ready to write down the linear system of equations for Eqs. (4.9) to (4.11), and develop a numerical algorithm to solve the system.
4.3.2 Existence of the solution
By using the staggered grid for the discretization of fluid variables, it is well-known that the matrix obtained by the discrete divergence operator of the fluid velocity can be written as the transpose of the discrete gradient operator of the pressure. As discussed in previous subsection, the resultant matrix obtained by the discrete surface divergence of velocity can be written as the transpose of the matrix obtained by the discrete spreading operator of the tension. Thus, the linear system for Eqs. (4.9)-(4.11) is symmetric and can be written as
where the sub-matrix L, G and S are represented the discrete Laplacian µ∆h, discrete gradient −∇h, and the discrete spreading operator Sh. The sub-matrix size of L, G and S are ((m−1)n+m(n−1))×((m−1)n+m(n−1)), ((m−1)n+m(n−1))×(mn) and (m − 1)n + m(n − 1) × M, respectively. The right-hand side vector [bc1, bc2, bc3]T of Eq. (4.17) consists only of the velocity boundary conditions since the pressure does not need the boundary condition in staggered grid formulation.
Let us discuss the existence of the linear system of Eq. (4.17). From now on, we denote the matrix in (4.17) by A. As is known, without the effect of the inextensible interface, the linear system becomes pure Stokes flow as
"
Let us denote the matrix in Eq. (4.18) by eA. It is also well-known that the nullity of eA equals to one since the pressure is unique up to a constant, and the existence of a solution can be verified by using the discrete incompressible constraint (4.10).
To be precisely, since the rank of deficiency of eA is only one, based on the algebraic structure of the sub-matrix G, the kernel of eA is
ker( eA) = span{[ 0 · · · 0| {z }
And for any vector z ∈ ker( eAT) = ker( eA), we have
= 0, (by the discrete incompressibility (4.10)) (4.20) which shows the compatibility condition for the existence of a solution.
If the effect of the inextensible interface is added, the matrix eA is augmented by S and ST to become A as in (4.17). Since the matrix S comes from the discrete spreading operator of the tension, the entries of S depend on the number of moving Lagrangian markers, their positions and tangents. It is unlikely to show rigorously that the nullity of A is exactly equal to one; however, we have checked the above statement to be true in our numerical experiments. So the apparent kernel will be
ker(A) = span{[ 0 · · · 0| {z }
and the existence of a solution for the linear system (4.17) follows the equality of (4.20) immediately.
4.3.3 Fractional step method
In this subsection, we follow the idea of fractional step method developed by Taira and Colonius [64] to solve the resultant linear system of equations (4.17). In [64], the authors applied the IB method to simulate the incompressible flow over solid bodies with prescribed body surface motion. Unlike the previous approaches in [13, 57, 30], they introduce the boundary force as another Lagrange multiplier to enforce the no-slip constraint for the velocity at the immersed boundary. From this point of view, the present approach shares the similar spirit as in [64] by introducing the elastic tension as a new Lagrange multiplier to enforce the surface divergence free constraint (4.11) along the interface. Since the pressure can be regarded as a Lagrange multiplier for the fluid divergence free constraint (4.10), one can group those two Lagrange multipliers as a new column vector φ = [p, σ]T and combine the sub-matrices G and S as Q = [G, S], the linear system (4.17) now becomes
As in [64], we perform a block LU decomposition to the above equation to obtain
Note that, the above decomposition is possible (L−1 exists) since the matrix L arising from the discrete Laplacian operator is symmetric and negative definite. This matrix decomposition is also referred to Uzawa method. One can further split the above matrix equation into following steps by introducing an intermediate velocity vector u∗ as
Lu∗ = bc1, (4.24)
(−QTL−1Q)φ = eb − QTu∗, (4.25)
u = u∗− L−1Qφ. (4.26)
Recall that the original matrix denoted by A in Eq. (4.17) is singular due to the pressure value is unique up to a constant, thus the singularity cannot be removed from applying the block LU decomposition. In fact, the matrix (−QTL−1Q) in Eq. (4.25) is symmetric and singular with rank mn + M − 1 since G is singular. We now provide an existence of solution for Eq. (4.25) as follows.
For any vector y ∈ ker((−QTL−1Q)T) = ker(−QTL−1Q), then the vector y satis-fies yT(−QTL−1Q)y = (Qy)T(−L−1)(Qy) = 0. This implies that Qy = 0 since −L−1 is a positive definite matrix; thus, y ∈ ker(Q). By writing L · 0 + Qy = 0, we can immediately obtain the vector [ 0 y ]T ∈ ker(A) = ker(AT). Notice that, such y also satisfies yT(eb− QTu∗) = 0 since yTeb = 0 due to the discrete incompressible constraint (4.10). Therefore, the right-hand side vector of Eq. (4.25) belongs to the range of the matrix (−QTL−1Q) in which a solution exists. One can also immediately see from the structure of sub-matrix G that the solution in Eq. (4.25) is unique up to a constant.
Now we are ready to describe the detailed numerical implementation for solving Eqs. (4.24)-(4.26). It is a common practice to avoid the direct computation of the inverse of the matrix L since it is too expensive. In [64], a second-order approximation for L−1based on Taylor expansion is implemented for solving similar equations as our Eqs. (4.25)-(4.26), and conjugate gradient method is applied to solve those equations iteratively. However, this leads to another time step constraint related to the viscosity and the eigenvalues of the discrete Laplacian. In this chapter, since we are working on the Stokes flow rather than the Navier-Stokes, we are unable to approximate L−1 using Taylor’s expansion. Although we do not approximate the L−1 directly, we still can solve Eqs. (4.25)-(4.26) efficiently thanks to the fast Poisson solver developed in public software package FISHPACK [1]. (The present matrix L is nothing but the discrete Laplacian operator.) The detailed steps for solving Eqs. (4.24)-(4.26) are as follows.
Step 1. Solve Eq. (4.24) by two fast Poisson solvers to obtain intermediate velocity field u∗.
Step 2. Since the matrix (−QTL−1Q) is symmetric and positive semi-definite, the conju-gate gradient method can be applied. In each iteration, a matrix-vector product (−QTL−1Q) ϕ is needed; fortunately, this can be done by letting z = L−1Qϕ, and solving Lz = Qϕ. Once it is done, we multiply z by −QT to obtain the product needed. Again, the time-consuming cost in each iteration is one fast Poisson solver.
Step 3. Find the velocity field u from Eq. (4.26). Since φ is solved via Step 2, by solving Lw = Qφ, we then obtain u = u∗− w. Again, this involves applying one fast Poisson solver.
Therefore, the overall cost in Step 1-3 for our present numerical algorithm can be counted in terms of the number of fast Poisson solver applied. In next section, we shall show the numbers of fast Poisson solver used in the Stokes flow for different grid resolutions.