A Cartesian grid-based boundary integral method for moving interface problems

This paper presents a stable and efficient Cartesian grid-based boundary integral method that reformulates elliptic and parabolic PDEs into boundary integral equations solved via matrix-free GMRES and finite difference-based integral evaluation, utilizing $\theta-L$ variables to simplify mesh preservation and enable robust time-stepping for complex moving interface problems like Hele-Shaw flow and Stefan solidification.

The Gist

Imagine you are watching a drop of oil spread through water, or a piece of ice growing in a glass of super-cooled water. In both cases, there is a moving "frontier" or boundary separating two different worlds (oil vs. water, ice vs. liquid). This is what scientists call a moving interface problem.

The challenge is that this boundary is constantly changing shape, and the physics happening on one side of the boundary affects the other side in a complex, back-and-forth dance. Calculating this on a computer is notoriously difficult because the boundary is jagged, moving, and often creates "stiff" math problems that crash standard calculators.

This paper introduces a new, clever way to solve these problems using a Cartesian grid-based boundary integral method. Here is a simple breakdown of how it works, using some everyday analogies.

1. The Problem: The "Shape-Shifting" Boundary

Think of the moving boundary (the edge of the oil drop or the ice crystal) as a shapeshifting rubber band.

• The Old Way: Traditional methods try to build a custom net (mesh) that fits perfectly around this rubber band every single time it moves. As the band twists and turns, the net gets tangled, and you have to constantly rebuild the whole net. It's like trying to take a photo of a dancing jellyfish by constantly reshaping your camera lens to fit its wiggles. It's slow and messy.
• The New Way: This paper suggests keeping a fixed, rigid grid (like a chessboard or graph paper) underneath the action. The rubber band just floats over the grid. We don't need to reshape the grid; we just need to figure out what the rubber band is doing at the points where it crosses the grid lines.

2. The Core Trick: "Kernel-Free" Boundary Integrals

Usually, to solve these problems, mathematicians use a "magic formula" (called a Green's function) that tells them how a single point of influence spreads out. But these formulas are often messy, involving "singularities" (mathematical infinities) that are hard to calculate, like trying to measure the exact temperature of a single, infinitely hot point.

The authors developed a "Kernel-Free" method.

• The Analogy: Imagine you want to know the temperature of a room, but you don't have the formula for how heat spreads. Instead of calculating the formula, you just solve the heat equation on your fixed grid.
• How it works: Instead of doing difficult, messy integrals (summing up infinite tiny points), the method turns the problem into a standard puzzle that a computer can solve very fast using tools like the Fast Fourier Transform (FFT). Think of FFT as a super-fast translator that turns a complex, wiggly sound wave into a simple list of notes. The method uses this translator to skip the messy math and get straight to the answer.

3. The "Stiffness" Problem: The Rubber Band's Snap

When surface tension (the "skin" of the fluid) is involved, the math becomes "stiff."

• The Analogy: Imagine the rubber band is made of a super-tight spring. If you try to move it with a simple, step-by-step approach (explicit time-stepping), the spring snaps back so violently that your simulation explodes. You'd have to take microscopic steps (like moving a snail) to keep it stable, which takes forever.
• The Solution: The authors use a technique called Small-Scale Decomposition (SSD) combined with a $\theta-L$ formulation.
  • $\theta-L$ Formulation: Instead of tracking the $x$ and $y$ coordinates of the rubber band, they track its length ($L$) and its angle ($\theta$). It's like describing a dancer by their speed and the angle of their spin, rather than the exact position of their left foot and right foot. This keeps the "dancer" from getting tangled up.
  • SSD: This technique separates the "easy" part of the movement from the "hard, stiff" part. They treat the stiff part (the springy tension) with a special "semi-implicit" step that stabilizes the spring, allowing them to take giant steps in time without the simulation blowing up.

4. What Did They Test?

They tested this new "fixed grid + smart math" system on two classic scenarios:

1. Hele-Shaw Flow: Watching an air bubble push through oil in a thin gap. They simulated the bubble growing, twisting, and forming "fingers" (viscous fingering). The method handled long, complex simulations without the grid getting confused.
2. The Stefan Problem: Watching ice grow in water. They simulated "dendritic growth" (snowflake patterns). They showed that their method could handle ice growing with or without water flowing around it, and even with buoyancy (hot water rising, cold water sinking).

The Bottom Line

This paper presents a universal toolkit for simulating moving boundaries.

• It's Fast: It uses fixed grids and fast math tricks (FFT) instead of rebuilding nets.
• It's Stable: It uses a special "angle and length" language to stop the math from exploding.
• It's Accurate: It avoids the messy, error-prone integrals that usually plague these simulations.

In short, they figured out how to watch a complex, shape-shifting dance on a fixed stage without ever having to move the stage or break the camera. This allows scientists to simulate things like crystal growth, oil recovery, and fluid mixing much faster and more accurately than before.

Technical Summary

1. Problem Statement

The paper addresses moving interface problems, a class of free boundary problems where the interface separating two phases evolves over time based on the solution of underlying Partial Differential Equations (PDEs). The authors focus on two representative examples:

• Hele-Shaw Flow: An elliptic problem describing the flow of viscous fluid in a thin gap, governed by the Poisson equation for pressure. It involves surface tension effects modeled by the Laplace-Young condition.
• Stefan Problem: A parabolic problem modeling solidification or melting (diffusion-driven interface motion). It involves coupled heat diffusion and fluid flow (Navier-Stokes) with convection, governed by the heat equation and modified Stokes equations.

Key Challenges:

• Complex Geometry: The interface evolves, making domain-fitted mesh methods (like FEM) computationally expensive due to frequent remeshing.
• Singularities: Traditional Boundary Integral Equation (BIE) methods require explicit Green's functions and complex quadrature to handle singular and nearly singular integrals.
• Stiffness: High-order curvature terms in the interface evolution (e.g., surface tension) introduce severe stiffness, imposing strict stability constraints on explicit time-stepping schemes.

2. Methodology

The authors propose a unified numerical framework combining Cartesian grid-based solvers with Boundary Integral Equations (BIE) and Small-Scale Decomposition (SSD).

A. Reformulation as Boundary Integral Equations

Instead of solving bulk PDEs on complex domains, the authors reformulate the problems into BIEs defined only on the interface $\Gamma$:

• Hele-Shaw: The Poisson equation is split into a singular source term and a harmonic function, leading to a Fredholm integral equation of the second kind for the dipole density.
• Stefan: Time-dependent PDEs are discretized in time (semi-implicit scheme) to reduce them to elliptic problems (Modified Helmholtz and Modified Stokes). These are then split into a regular part (solved via finite differences) and an interface part (solved via BIE).

B. Kernel-Free Boundary Integral (KFBI) Method

To solve the resulting BIEs without evaluating singular integrals:

• Equivalent Interface Problems: The boundary integral operators (single-layer and double-layer potentials) are treated as solutions to equivalent interface problems with jump conditions across $\Gamma$.
• Cartesian Grid Solvers: These interface problems are solved on a fixed Cartesian grid using Finite Difference Methods (FDM) with correction terms near the interface to handle jumps in derivatives.
  • For the Modified Helmholtz equation, a corrected 5-point central difference scheme is used.
  • For the Modified Stokes equation, a corrected Marker-and-Cell (MAC) scheme on a staggered grid is employed.
• Matrix-Free GMRES: The linear systems arising from the discretization of the BIEs are solved using the Generalized Minimal Residual (GMRES) method. Matrix-vector products are computed by solving the equivalent interface PDEs using fast solvers (FFT for Helmholtz, Geometric Multigrid for Stokes), avoiding the explicit formation of dense matrices.
• Kernel-Free: The method does not require the analytical form of the Green's function for the integral evaluation, only for imposing artificial boundary conditions on unbounded domains.

C. Interface Evolution and Small-Scale Decomposition (SSD)

To evolve the interface stably and efficiently:

• $\theta-L$ Formulation: The interface is parameterized by the tangent angle $\theta$ and arc-length $L$ rather than Cartesian coordinates $(x, y)$. This prevents mesh clustering and maintains uniform point distribution.
• Small-Scale Decomposition: The normal velocity $U$ is decomposed into a stiff linear part (dominated by curvature) and a non-stiff nonlinear part.
  • The stiff part (involving Hilbert transforms and high-order derivatives) is treated implicitly in the Fourier space.
  • The non-stiff part is treated explicitly (Adams-Bashforth).
• Semi-Implicit Time Stepping: This approach removes the severe time-step restrictions imposed by surface tension stiffness while avoiding the complexity of solving fully nonlinear algebraic systems.

3. Key Contributions

1. Unified Framework: The first Cartesian grid-based framework that successfully combines Kernel-Free Boundary Integral (KFBI) methods with Small-Scale Decomposition (SSD) for both elliptic (Hele-Shaw) and parabolic (Stefan) moving interface problems.
2. Avoidance of Singular Quadrature: By using the KFBI approach, the method eliminates the need for evaluating singular and nearly singular integrals, which are notoriously difficult in traditional BIE methods.
3. Efficiency and Scalability: The use of fast PDE solvers (FFT and Multigrid) on fixed Cartesian grids allows for $O(N \log N)$ or $O(N)$ complexity, significantly reducing computational cost compared to remeshing methods.
4. Robust Stability: The semi-implicit $\theta-L$ formulation with SSD enables large time steps, making long-time simulations of complex phenomena (like dendritic growth and viscous fingering) feasible and stable.

4. Numerical Results

The method was validated through extensive numerical experiments:

• Convergence: Both Hele-Shaw and Stefan problems demonstrated second-order accuracy in space and time.
• Hele-Shaw Simulations:
  • Successfully simulated bubble relaxation (area-preserving, curve-shortening).
  • Captured unstable viscous fingering with varying surface tension, showing the competition between driving forces and surface tension.
  • Performed long-time simulations (up to $T=203$) of large bubbles with intricate finger-like patterns, maintaining mesh quality.
• Stefan Simulations:
  • Grid Refinement: Showed superior convergence and reduced grid-induced anisotropy compared to Level Set and Front Tracking methods.
  • Stability: The semi-implicit scheme remained stable with time steps orders of magnitude larger than explicit schemes (e.g., $\tau=0.01$ vs. $\tau=2.5\times 10^{-5}$).
  • Dendritic Growth: Accurately reproduced dendritic patterns with 4-fold and 6-fold anisotropy. Tip velocities matched theoretical solvability predictions.
  • Convection Effects: Simulated dendritic growth with external flow and buoyancy-driven flow, correctly capturing asymmetric growth patterns and temperature distributions caused by fluid transport.

5. Significance

This work represents a significant advancement in computational fluid dynamics and materials science simulation. By bridging the gap between the efficiency of Cartesian grid methods and the accuracy of boundary integral formulations, it provides a robust tool for simulating complex free boundary problems. The ability to handle variable coefficients, anisotropic surface tension, and coupled fluid-thermal systems without remeshing or singular quadrature makes this method highly versatile for applications ranging from crystal growth to biological tumor modeling. The authors note that while currently limited to 2D, the framework lays the groundwork for future 3D extensions.