A Higher Order Unfitted Space-Time Finite Element Method for Coupled Surface-Bulk problems

TL;DR

This work develops a higher-order unfitted space-time finite element method for coupled convection-diffusion problems on moving surface-bulk domains. It combines a bulk unfitted approach with a surface counterpart and employs a time-dependent isoparametric mapping $\Theta_h^{\text{st}}$ to achieve high-order geometric accuracy, enabling tensor-product space-time discretizations of arbitrary order. The method is stabilized by a space-time Ghost penalty and a surface normal gradient term, and is applicable to both linear (Henry) and nonlinear (Langmuir) couplings, with numerical experiments showing convergence of order $k+1$ in the final-time $L^2$ norm on bulk and surface. The approach provides a robust framework for accurately simulating evolving coupled systems (e.g., biological cells or multi-phase flows) while maintaining computational efficiency at high orders, and it opens avenues for rigorous analysis and broader applicability to other physical problems.

Abstract

We present a higher order space-time unfitted finite element method for convection-diffusion problems on coupled (surface and bulk) domains. In that way, we combine a method suggested by Heimann, Lehrenfeld, Preuß (SIAM J. Sci. Comput. 45(2), 2023, B139 - B165) for the bulk case with a method suggested by Sass, Reusken (Comput. Math. Appl. 146(15), 2023, 253-270) for the surface case. The geometry is allowed to change with time, and the higher order discrete approximation of this geometry is ensured by a time-dependent isoparametric mapping. The space-time discretisation approach allows for straightforward handling of arbitrary high orders. In that way, we also generalise results of Hansbo, Larson, Zahedi (Comput. Methods Appl. Mech. Engrg. 307, 2016, 96-116) to higher orders. The convergence of the proposed higher order discretisations is confirmed numerically.

Figures (5)

• Figure 1: Illustration of the decomposition of the fixed domain $\tilde{\Omega}$ into $\Omega_1(t)$ and $\Omega_2(t)$. For the surface $\Gamma(t) = \partial \Omega_2(t)$.
• Figure 2: Illustration of the discrete regions $\mathcal{E}(Q^{\text{lin}, n})$, $\mathcal{E}(\Gamma^{\text{lin}, n})$ in a one-dimensional example.
• Figure 3: Plots of the manufactured solution $u_B$ on the bulk, and $u_S$ on the surface. They respectively have their colour bar, indicating some coupling between higher/ lower concentration values. Left: $u_B(x,0)$, $u_S(x,0)$. Right: $u_B(x,T=0.5)$, $u_S(x,T=0.5)$.
• Figure 4: Numerical convergence results for the linear model system ($b_{BS} = 0$) and different polynomial orders $k=k_s=k_t=q_s=q_t$. Refinement by $i$ refines both space and time simultaneously. Left: The discrete error is measured in bulk and surface norms vs. refinements $i$. In all cases, the expected higher order convergence is observed. Right: The error on the bulk domain is plotted on the y-axis, whilst now the x-axis shows the runtime of the computation. The higher order methods lead to shorter runtimes when high accuracy is demanded.
• Figure 5: Numerical convergence results for the non-linear model system ($b_{BS} = 1$), where the nonlinearity is solved by a Newton solver. Left: Numerical errors in $L^2$ norms at final time on bulk and surface. Right: Numbers of Newton iterations required. For coarse meshes/ large timesteps, higher iteration numbers might be needed for higher order.