On immersed boundary kernel functions: a constrained quadratic minimization perspective

In the immersed boundary (IB) approach to fluid-structure interaction modeling, the coupling between the fluid and structure variables is mediated using a regularized version of Dirac delta function. In the IB literature, the regularized delta functions, also referred to IB kernel functions, are either derived analytically from a set of postulates or computed numerically using the moving least squares (MLS) approach. Whereas the analytical derivations typically assume a regular Cartesian grid, the MLS method is a meshless technique that can be used to generate kernel functions on complex domains and unstructured meshes. In this note we take a viewpoint that IB kernel generation, either analytically or via MLS, is a constrained quadratic minimization problem. The extremization of a constrained quadratic function is a broader concept than kernel generation, and there are well-established numerical optimization techniques to solve this problem. For example, we show that the constrained quadratic minimization technique can be used to generate one-sided (anisotropic) IB kernels and/or to bound their values.

READ FULL TEXT VIEW PDF

Authors

page 9

04/15/2021

A one-sided direct forcing immersed boundary method using moving least squares

This paper presents a one-sided immersed boundary (IB) method using kern...
05/30/2021

On the Lagrangian-Eulerian Coupling in the Immersed Finite Element/Difference Method

The immersed boundary (IB) method is a non-body conforming approach to f...
03/22/2022

Isoparametric singularity extraction technique for 3D potential problems in BEM

To solve boundary integral equations for potential problems using colloc...
03/29/2021

An image-incorporated immersed boundary method for diffusion equations

A novel sharp interface ghost-cell based immersed boundary method has be...
11/30/2021

An Analysis of the Numerical Stability of the Immersed Boundary Method

We present a numerical stability analysis of the immersed boundary(IB) m...
08/25/2016

Minimizing Quadratic Functions in Constant Time

A sampling-based optimization method for quadratic functions is proposed...
01/22/2021

The art of coarse Stokes: Richardson extrapolation improves the accuracy and efficiency of the method of regularized stokeslets

The method of regularised stokeslets is widely used in microscale biolog...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The immersed boundary (IB) method was originally proposed by Peskin to model blood flow around heart valves Peskin1977 . In the IB formulation, the fluid occupying the entire computational domain is described by Eulerian variables, whereas the immersed structure that is allowed to freely cut the background Eulerian grid is described by Lagrangian variables. The presence of the immersed structure leads to an additional body force term in the momentum equation. The Lagrangian structure moves with the local fluid velocity and satisfies the no-slip boundary condition on the fluid-structure interface. The Eulerian-Lagrangian coupling in the IB formulation is mediated using regularized integral kernels. Specifically, IB forces defined on the marker points are spread to the Eulerian grid from the Lagrangian structure, and conversely, the fluid velocity is interpolated from the Eulerian mesh to the marker points. The use of an isotropic and regularized IB kernel diffuses

the fluid-structure interface to nearby grid cells and this smearing of the interface is proportional to the width of the kernel. This also implies that the (usual) isotropic IB kernels couple the structural degrees of freedom to fluid variables on both sides of the interface.

In the literature, IB kernels have been derived either analytically or computed numerically via the moving least squares (MLS) approach. Peskin Peskin02

used a set of postulates to derive analytical expressions for several IB kernels such as the four-, five-, and six-point function. These postulates aim to achieve several desirable features in an IB kernel, such as polynomial reproducing conditions (also known as the discrete moment conditions), grid translational invariance property, and compact kernel support. More recently, Bao et al. 

Bao2016

used a weaker second moment condition to derive improved versions of the standard five-, and six-point IB kernels. A weaker second moment condition allowed the authors to remove the negative tails of the standard five- and six-point kernels, and to enhance their gird translational invariance property. A main limitation of the analytical approach is the assumption of a regular Cartesian grid and the tensor-product form of the kernel function in dimensions higher than one. To circumvent these limitations, Vanella and Balaras introduced the MLS approach to generate IB kernel functions 

Vanella09 . Although the authors in Vanella09 demonstrated the IB/MLS approach on a regular Cartesian grid, the technique is easily extendable to unstructured grid-based IB methods as shown by Krishnan et al. Krishnan17 and Saadat et al. Saadat18 . In this work we take a viewpoint that IB kernel generation is a constrained quadratic minimization problem and both analytical and MLS formulations can be described using the language of optimization.

In the prior IB/MLS works Vanella09 ; Li15 ; Tullio16 ; Krishnan17 ; Le2017 ; Saadat18 ; Haji19 , a standard formulation of the moving least squares method was used to transform a cubic spline weight function into a regularized IB kernel function that satisfies zeroth- and first-order discrete moment conditions. The standard MLS formulation minimizes the weighted -norm of the difference between the interpolant and the given data values. A lesser known formulation of the moving least squares problem (in the IB literature) is the Backus-Gilbert theory, which was developed in the context of geophysics Backus1968resolving . In the Backus-Gilbert formulation, polynomial reproducing conditions are used in the MLS problem definition. However, both formulations ultimately lead to the same IB kernel function as shown in this work. The connection between the two MLS formulations was first pointed out by Bos and Šalkauskas in Bos1989moving and re-iterated here in the context of constrained quadratic minimization. We also show that the constrained quadratic minimization approach can be used to generate discrete moment conditions satisfying anisotropic IB kernels that effectively couple structural degrees of freedom to fluid variables on only one side of the fluid-structure interface. This can help to reduce spurious feedback forcing and “leakage” that are typically observed in the diffuse-interface IB models. Furthermore, the constrained quadratic minimization method can be used to bound the weights of the IB kernels.

In the following sections, we first review the two formulations of the moving least squares problem and describe them using the language of optimization. Next, we review the set of postulates used by Peskin to derive the four-point IB kernel and frame it into the language of constrained quadratic minimization. Finally, we describe methods to generate one-sided and/or bounded IB kernels and demonstrate this process using a two-dimensional example.

2 Moving least squares method

2.1 The standard formulation

The standard formulation of the moving least squares method considers the following approximation

problem. Given a (column) vector of data values

, defined on distinct data sites by a (smooth) function , find the quasi-interpolant to

(1)

at an evaluation point using the polynomial approximation space . Typically, and the coefficient vector in Eq. (1) is found by minimizing the weighted -norm of the error function

(2)

with respect to . In the equation above, is a given weight function which decreases in magnitude away from the evaluation point and is taken to be non-negative. The error function is minimized using the essential condition of extremum, , which leads to the normal equation of the form

(3)

Here, denotes the weighted inner product of two quantities. If denotes the polynomial matrix with entries , and is the diagonal matrix containing weights (with as the main diagonal of ), then in matrix notation Eq. (3) reads as

(4)

Here, is the positive-definite Gram matrix and is the right-hand side vector . Note that the coefficient vector , Gram matrix , and the right hand side vector , all depend on the evaluation point . Consequently, the normal Eq. (4) is solved anew each time the evaluation point is changed or “moved.”

In the context of the immersed boundary method, the standard MLS formulation is used to interpolate velocity from Eulerian grid nodes to a Lagrangian marker point and to spread IB force from a Lagrangian marker point to the surrounding Eulerian nodes Vanella09 . This is achieved by considering a particular Lagrangian marker point as the current evaluation point and the surrounding Eulerian grid nodes (to the Lagrangian point) as the distinct data sites that provide velocity or force data. The IB/MLS kernel function , also known as the generating function, is obtained by combining Eqs. (1) and (3)

(5)

Therefore, from the equation above we get

(6)

Comparing the two forms of the quasi-interpolant in Eqs. (1) and (6), it is clear that the latter form employs a kernel function that depends only on the evaluation point . In contrast, the former form uses a coefficient vector , which depends on both the evaluation point and the data values , i.e., . For IB methods, it is therefore more convenient to work with directly, as it allows the same kernel function to be used in both velocity interpolation and force spreading operations. It can be shown that using the same IB kernel function in the two operators conserve energy during Lagrangian-Eulerian coupling Peskin02 .

Although, it is not explicitly enforced in the error minimization of , the standard MLS method will reproduce the polynomial basis functions , i.e.,

(7)

We shall prove this property of standard MLS formulation in Sec. 2.2.

Next, we show that the standard MLS formulation can be viewed as a constrained quadratic minimization problem of the form

(8)

Note that Problem A is a minimization problem because is a (symmetric) positive-definite matrix. The above constrained optimization problem can be solved by introducing the Lagrangian function , with the help of a Lagrange multiplier

(9)

Using the essential conditions of extrema, and , we get the following equalities

(10)
(11)

respectively. Solving the above two equations, we obtain . Hence, the solutions to the standard MLS formulation and Problem A are the same.

2.2 The Backus-Gilbert formulation

The Backus-Gilbert formulation of the moving least squares method seeks the quasi-interpolant of the type

(12)

that satisfies the polynomial reproducing conditions

(13)

The polynomial reproduction constraints correspond to the discrete moment conditions for the function . Eq. (13) in matrix form reads as , in which is the same polynomial matrix as defined in the previous section 2.1. The desired generating function in the Backus-Gilbert formulation is obtained by solving the following constrained quadratic minimization problem:

(14)

The solution to Problem B is obtained by minimizing the Lagrangian function , given by

(15)

in which is the Lagrange multiplier. The Lagrangian is minimized by using the essential conditions of extrema, and , which yields the following equalities

(16)
(17)

respectively. Solving the above two equations, we obtain

(18)
(19)

Comparing Eqs. (19) and (6), it is clear that both formulations of the moving least squares method yield the same quasi-interpolant to function . Since the polynomial reproduction constraints are directly included in the quadratic function minimization used to solve Problem B, the equivalence of the two formulations implies that the standard moving least squares method also reproduces the polynomial basis functions.

In the IB literature, constant and linear polynomials are typically employed in the moving least squares formulation. Defining the polynomials relative to the evaluation point , the set of basis functions for IB problems consists of

(20)

Note that if the weight function already satisfies the discrete moment conditions on the support of the IB/MLS kernel, then the constrained quadratic minimization would result in a trivial solution . Consequently, the product in Eq. (19) produces a vector whose entries are all one.

3 Peskin’s IB kernel

If the IB simulations are performed on a regular Cartesian grid (which is often the case), then Peskin’s delta functions can be used instead of the IB/MLS kernels described earlier. Peskin’s delta functions satisfy the polynomial reproducing conditions, can be expressed compactly in a piecewise analytical manner, and furthermore, they do not explicitly depend upon a particular evaluation/marker point. The latter property is due to the fact that Peskin’s delta functions are constructed in terms of relative distance between the data site and the evaluation point, i.e, and .

Taking the example of widely popular four-point delta function from the IB literature, we list the four postulates used by Peskin for its construction below:

Note that the even-odd condition in the second postulate implies the zeroth moment condition (or the constant reproducing condition), and therefore, the first equation is not explicitly used in deriving . The first moment condition in the third postulate is the linear polynomial reproducing condition, whereas the sum-of-squares condition in the last postulate is a weaker condition on grid translational invariance.

Similar to the IB/MLS kernel formulation, Peskin’s four-point delta function can be formulated as a constrained quadratic minimization problem. The minimization problem reads as

(21)

in which is a (column) vector containing the values one and zero at even- and odd-indexed data sites, respectively. is defined analogously. Note that Problem C has a nonlinear equality constraint, which would require a nonlinear optimizer, such as the interior point algorithm Wachter2006 or the sequential quadratic programming technique Nocedal2006 .

4 One-sided and bounded IB kernels

Next, we show that the constrained quadratic minimization approach can help to generate IB kernels that effectively couple structural degrees of freedom to fluid variables on only one side of the fluid-structure interface. Such interpolation kernels have been termed “one-sided IB kernels” by Bale et al. Bale2021 , who have recently proposed them. The one-sided IB kernels help to reduce spurious feedback forcing and internal flows that are typically observed in IB models that use isotropic kernel functions which couple the structure to fluid degrees of freedom on both sides of the interface Bale2021 . However, such kernels have large positive and negative weights, which can lead to numerical instability over time. To mitigate the numerical instability the authors in Bale2021 mollified the one-sided IB kernel weights at the expense of the first moment condition. In contrast, here we show that the numerical values of the one-sided kernel function can be bounded, e.g., made positive or constrained within a smaller range using the constrained quadratic minimization approach, while still satisfying the first moment condition.

(a) Support of an isotropic IB kernel
(b) Support of a one-sided IB kernel
Figure 1: Representative Lagrangian markers () on a circular interface embedded in a Cartesian grid.LABEL: Full support of a typical isotropic IB kernel and LABEL: restricted support of one-sided IB kernel. The Eulerian support of the maker points are shown by circles ().

We demonstrate the one-sided IB kernel construction by considering a circular interface that is embedded in a regular Cartesian grid. The interface is discretized using Lagrangian marker points, four of which are highlighted by red markers in Fig. 1. Typical isotropic IB kernels consider both sides of the interface to define the support region; this is shown by black circles for the four marker points in Fig. 1(a). For the one-sided kernel, the Eulerian support region is restricted to only one side of the interface. For example, in Fig. 1(b) the kernel support is restricted to the exterior of the cylinder, which is denoted region. Analogously, the kernel support can be defined interior to the cylinder ( region); here we discuss the approach in the context of the former case. The one-sided, bounded, and discrete moment (polynomial reproducing) conditions satisfying IB kernels can be constructed numerically by solving a constrained quadratic minimization problem of the type

(22)

In the Problem D above, is the restricted version of the weight function that is defined to be

(23)

The generating function in Eq. (22) is bounded on the lower end by values and on the upper end by values. In lack of any particular requirement, all weights are bounded alike, i.e., and , in which and are two scalars and is a column vector with all entries as one. Note that we have included an inequality constraint, along with the equality constraint in Problem D. Without the inequality constraint, Problem D is same as Problem A or Problem B.

Next, we solve Problem D for the circular interface example by taking constant and linear polynomial basis functions in (as in Eq. (20)). A linear field of the type

(24)

is defined on the cell centers of a uniform Cartesian mesh that discretizes the computational domain of extents . The uniform cell size is taken to be . As the circular interface partitions the domain into inner and outer regions, . The cylinder is centered around the point and has a radius of . The four marker points on the periphery of the cylinder make an angle of , , and with the -axis. A tensor-product form of the six-point spline kernel is used as a weight function in Problem D. The one-dimensional version of function reads as

(25)

In Eq. (25) above, and . Note that the weight function satisfies constant and linear polynomial reproducing conditions when the support of the Lagrangian marker is included on both sides of the interface.

(a) Case 1
(b) Case 2
(c) Case 3
(d) Case 4
Figure 2: Color plot of for four representative Lagrangian markers () placed on a circular interface of radius . is obtained using the constrained quadratic minimization formulation of Problem D. The interface is embedded in a computational domain of extents . Interpolation weights when LABEL: Eulerian support on both sides of the interface and only equality constraints are used; LABEL: Eulerian support restricted to region and only equality constraints are used; LABEL: Eulerian support restricted to region and both equality and inequality constraints are used; and LABEL: Eulerian support restricted to region and both equality and inequality constraints are used.

The linear field is interpolated onto the four marker points using

  • Case 1: Eulerian support on both sides of the interface and only equality constraints in Problem D;

  • Case 2: Eulerian support restricted to region and only equality constraints in Problem D;

  • Case 3: Eulerian support restricted to region and both equality and inequality constraints in Problem D; and

  • Case 4: Eulerian support restricted to region and both equality and inequality constraints in Problem D.

The three cases listed above are solved using the quadprog routine in MATLAB using the default solver options. The quadprog routine requires defining the cost function in terms of a matrix ( in this case) and lower and upper bound vectors ( and , respectively in this case) for the variable to be minimized ( in this case).

Table 1 lists the relative error in the interpolated values of and Fig. 2 shows the color plot of the interpolation weights around the four marker points. Fig. 2(a), which plots the weights for Case 1 shows that Problem D generates the usual isotropic IB kernel because the support of the marker points is defined on both sides of the interface. In fact for Case 1, Problem D did not modify the underlying weight function because the weight function already satisfies the polynomial reproducing conditions, as mentioned earlier. This can be verified by computing the - and -norms of the difference , which come out to be and , respectively.

Interpolation weights for Case 2 are shown in Fig. 2(b). In this case we observe large positive weights near the marker points, which are followed by large negative values; ranges from -0.3627 to 0.9178 in Case 2. This kind of weight distribution (having large positive and negative values) was shown to have a destabilizing effect when used within a direct-forcing IB method Bale2021 . To make use of the one-sided IB kernel, the authors in Bale2021 mollified the weights through a renormalization procedure. However, the mollified weights satisfy only the zeroth moment condition, i.e, they are only first-order accurate.

Next, we show that Problem D can be used to systematically mollify the non-smooth weights obtained in Case 2 by prescribing the upper and lower bounds, without having to forfeit the first moment condition. The resulting smooth weights for Case 3 are shown in Fig. 2(c), wherein we take and , which is a subset of the range observed for in Case 2. To quantitatively access the accuracy of the IB kernels, we tabulate the interpolation error at the four marker locations in Table 1. From the table, it can be inferred that close to machine precision error norms are obtained for Case 1 and 2, which indicates that the linear field is interpolated exactly onto the Lagrangian marker points. This is not a surprise, because the linear polynomial reproducing constraint is imposed in the problem formulation itself, but it demonstrates the procedure to obtain discrete moment conditions satisfying one-sided kernels using quadratic programming, which is different compared to the moving-least squares technique of Problem A or B. For Case 3, is in the range of . We remark that close to machine precision error norms can be achieved in this case by further decreasing (increasing) the lower (upper) bound of . However, this will not lead to any practical difference in the results of an IB simulation and loosing some digits of interpolation accuracy should be preferred over accurate but non-smooth weights, such as those shown in Fig. 2(b).

Lastly, we consider Case 4, in which our objective is to completely eliminate the negative weights of Case 2. To achieve this, we take and . Table 1 shows the error norms for this case. Here, we observe that although the error norms are small (in the range of ), they are not close to machine precision, indicating that the linear field is not interpolated exactly. This is because the quadprog routine is unable to find weights that can satisfy the equality constraints exactly, while simultaneously respecting the imposed bounds; changing the solver settings did not affect the final convergence of quadprog (data not shown). On the other hand, if we observe the min/max limits of the color bar in Fig. 2(d), for Case 4 stays within the imposed limits. Note that although the one-sided weights of Case 4 are positive, they are not smooth and are highly localized near the interface. This type of weight distribution can also lead to numerical instability for the direct-forcing IB method Uhlmann2005 . Apparently, the weight-shifting technique of Bale et al. to completely eliminate the negative weights at the expense of first moment condition is more suitable for the direct-forcing IB approach than the inequality constraint employed in Case 4. Nevertheless, if we allow small magnitude of negative weights in the constrained quadratic formulation, we can achieve smoothness and as well as the desired accuracy in the one-sided IB kernels (Case 3).

Case Relative error at marker location
40 140 230 310
1
gray!10 2 gray!10 gray!10 gray!10 gray!10
3
gray!10 4 gray!10 gray!10 gray!10 gray!10
Table 1: Relative error for the three cases of Problem D. For Case 3, and and for Case 4, and

5 Summary and Conclusions

In this note we took a viewpoint that IB kernel generation, either analytically or via MLS, is a constrained quadratic minimization problem. The extremization of a constrained quadratic function is a broader concept than kernel generation and there are well-established numerical optimization techniques to solve this problem. We also reviewed the two formulations of the moving least squares problem and described them using the language of optimization. Further, the set of postulates used by Peskin to derive the four-point IB kernel was also described using this language. Finally, we described methods to generate one-sided and/or bounded IB kernels and demonstrated the process using a two-dimensional example. Four constrained quadratic minimization problems were described, namely, the Problem A, B, C, and D. All but Problem C involved linear constraints and can be solved using efficient optimization techniques such as the quadratic programming technique. For Problem C, a nonlinear optimizer, such as the interior point algorithm or the sequential quadratic programming technique is required.

Acknowledgements

A.P.S.B acknowledges helpful discussion with Rahul Bale.

Bibliography

References

  • (1) C. S. Peskin, Numerical analysis of blood flow in the heart, Journal of computational physics 25 (3) (1977) 220–252.
  • (2) C. S. Peskin, The immersed boundary method, Acta Numer 11 (2002) 479–517.
  • (3) Y. Bao, J. Kaye, C. S. Peskin, A gaussian-like immersed-boundary kernel with three continuous derivatives and improved translational invariance, Journal of Computational Physics 316 (2016) 139–144.
  • (4) M. Vanella, E. Balaras, Short note: A moving-least-squares reconstruction for embedded-boundary formulations, Journal of Computational Physics 228 (18) (2009) 6617–6628.
  • (5) S. Krishnan, E. S. Shaqfeh, G. Iaccarino, Fully resolved viscoelastic particulate simulations using unstructured grids, Journal of Computational Physics 338 (2017) 313–338.
  • (6) A. Saadat, C. J. Guido, G. Iaccarino, E. S. Shaqfeh, Immersed-finite-element method for deformable particle suspensions in viscous and viscoelastic media, Physical Review E 98 (6) (2018) 063316.
  • (7) D. Li, A. Wei, K. Luo, J. Fan, An improved moving-least-squares reconstruction for immersed boundary method, International Journal for Numerical Methods in Engineering 104 (8) (2015) 789–804.
  • (8) M. D. de Tullio, G. Pascazio, A moving-least-squares immersed boundary method for simulating the fluid–structure interaction of elastic bodies with arbitrary thickness, Journal of Computational Physics 325 (2016) 201–225.
  • (9) D.-V. Le, B.-C. Khoo, A moving-least-square immersed boundary method for rigid and deformable boundaries in viscous flow, Communications in Computational Physics 22 (4) (2017) 913–934.
  • (10) M. Haji Mohammadi, F. Sotiropoulos, J. Brinkerhoff, Moving least squares reconstruction for sharp interface immersed boundary methods, International Journal for Numerical Methods in Fluids 90 (2) (2019) 57–80.
  • (11) G. Backus, F. Gilbert, The resolving power of gross earth data, Geophysical Journal International 16 (2) (1968) 169–205.
  • (12) L. Bos, K. Salkauskas, Moving least-squares are backus-gilbert optimal, Journal of Approximation Theory 59 (3) (1989) 267–275.
  • (13) A. Wächter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical programming 106 (1) (2006) 25–57.
  • (14) J. Nocedal, S. J. Wright, Sequential quadratic programming, Numerical optimization (2006) 529–562.
  • (15) R. Bale, A. P. S. Bhalla, B. E. Griffith, M. Tsubokura, A one-sided direct forcing immersed boundary method using moving least squares, Journal of Computational Physics 440 (2021) 110359.
  • (16) M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, Journal of Computational Physics 209 (2) (2005) 448–476.