A spatially adaptive high-order meshless method for fluid-structure interactions

TL;DR

This paper introduces a high-order meshless method with adaptive refinement for fluid-structure interactions, effectively handling singularities and complex geometries to improve simulation accuracy.

Contribution

It develops a spatially adaptive GMLS scheme with an a posteriori refinement strategy for high-fidelity fluid-structure interaction simulations.

Findings

Achieves optimal convergence near singularities.

Effectively resolves complex colloid geometries.

Demonstrates high accuracy in fluid-structure simulations.

Abstract

We present a scheme implementing an a posteriori refinement strategy in the context of a high-order meshless method for problems involving point singularities and fluid-solid interfaces. The generalized moving least squares (GMLS) discretization used in this work has been previously demonstrated to provide high-order compatible discretization of the Stokes and Darcy problems, offering a high-fidelity simulation tool for problems with moving boundaries. The meshless nature of the discretization is particularly attractive for adaptive h-refinement, especially when resolving the near-field aspects of variables and point singularities governing lubrication effects in fluid-structure interactions. We demonstrate that the resulting spatially adaptive GMLS method is able to achieve optimal convergence in the presence of singularities for both the div-grad and Stokes problems. Further, we present a series of simulations for flows of colloid suspensions, in which the refinement strategy efficiently achieved highly accurate solutions, particularly for colloids with complex geometries.

Figures

Equations

$$\{\begin{array}{ll}-\nabla \cdot (\kappa \nabla \varphi )=f& \text{for }\mathbf{x}\in {\Omega }_f\\ \varphi ={\varphi }_{\Gamma }& \text{for }\mathbf{x}\in \Gamma \end{array},$$

$$\{\begin{array}{ll}\frac{\nabla p}{\rho }-\nu {\nabla }^2\mathbf{u}=\mathbf{f}& \text{for }\mathbf{x}\in {\Omega }_f\\ \nabla \cdot \mathbf{u}=0& \text{for }\mathbf{x}\in {\Omega }_f\\ \mathbf{u}=\mathbf{w}& \text{for }\mathbf{x}\in {\Gamma }_w\\ \mathbf{u}={\dot{\mathbf{X}}}_n+{\dot{\mathbf{\Theta }}}_n\times (\mathbf{x}-{\mathbf{X}}_n)& \text{for }\mathbf{x}\in {\Gamma }_n,n=1,\dots ,N_c,\end{array}$$

$$\{\begin{array}{ll}\frac{\nabla p}{\rho }+\nu \nabla \times \nabla \times \mathbf{u}=\mathbf{f}& \text{for }\mathbf{x}\in {\Omega }_f\\ \frac{{\nabla }^{2}p}{\rho }=\nabla \cdot \mathbf{f}& \text{for }\mathbf{x}\in \Omega \\ \mathbf{u}=\mathbf{w}& \text{for }\mathbf{x}\in {\Gamma }_w\\ \mathbf{u}={\dot{\mathbf{X}}}_n+{\dot{\mathbf{\Theta }}}_n\times (\mathbf{x}-{\mathbf{X}}_n)& \text{for }\mathbf{x}\in {\Gamma }_n,n=1,\dots ,N_c\\ \mathbf{n}\cdot \frac{\nabla p}{\rho }+\nu \mathbf{n}\cdot \nabla \times \nabla \times \mathbf{u}=\mathbf{n}\cdot \mathbf{f}& \text{for }\mathbf{x}\in \Gamma .\end{array}$$

$$\{\begin{array}{l}{M}_n{\ddot{\mathbf{X}}}_n={\int }_{{\Gamma }_n}\mathbf{\sigma }\cdot d\mathbf{A}+{\mathbf{F}}_{e,n}\\ I_n{\ddot{\mathbf{\Theta }}}_n={\int }_{{\Gamma }_n}\left(\mathbf{x}-{\mathbf{X}}_n\right)\times (\mathbf{\sigma }\cdot d\mathbf{A})+{\mathbf{T}}_{e,n},\end{array}$$

$$\begin{array}{rl}& {\int }_{{\Gamma }_n}\mathbf{\sigma }\cdot d\mathbf{A}+{\mathbf{F}}_{e,n}=0\\ & {\int }_{{\Gamma }_n}\left(\mathbf{x}-{\mathbf{X}}_n\right)\times (\mathbf{\sigma }\cdot d\mathbf{A})+{\mathbf{T}}_{e,n}=0.\end{array}$$

$$J({\mathbf{x}}_i)=\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}{\left[{\psi }_j-{\mathbf{P}}_i^{⊺}({\mathbf{x}}_j){\mathbf{c}}_i\right]}^2{W}_{ij}.$$

$${\psi }_h(\mathbf{x})={\mathbf{P}}_i^{⊺}(\mathbf{x}){\mathbf{c}}_i^{\ast }\text{with}{\mathbf{c}}_i^{\ast }={\left(\sum _{k\in {\mathcal{N}}_{{ϵ}_i}}{\mathbf{P}}_i({\mathbf{x}}_k)W_{ik}{\mathbf{P}}_i^{⊺}({\mathbf{x}}_k)\right)}^{-1}\left(\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}{\mathbf{P}}_i({\mathbf{x}}_j)W_{ij}{\psi }_j\right).$$

$$D^{\alpha }\psi (\mathbf{x})\approx D_h^{\alpha }{\psi }_h(\mathbf{x}):=(D^{\alpha }\mathbf{P}(\mathbf{x}){)}^{⊺}{\mathbf{c}}^{\ast }.$$

$$W_{ij}=\frac{W_{{ϵ}_i}(r_{ij})+W_{{ϵ}_j}(r_{ij})}{2}.$$

$${\mathbf{P}}_i^{div}(\mathbf{x})=\left[\left(\genfrac{}{}{0ex}{}{1}{0}\right),\left(\genfrac{}{}{0ex}{}{0}{1}\right),\left(\genfrac{}{}{0ex}{}{0}{-2\frac{x-x_i}{{ϵ}_i}}\right),\left(\genfrac{}{}{0ex}{}{\frac{x-x_i}{{ϵ}_i}}{-\frac{y-y_i}{{ϵ}_i}}\right),\left(\genfrac{}{}{0ex}{}{2\frac{y-y_i}{{ϵ}_i}}{0}\right),\left(\genfrac{}{}{0ex}{}{0}{-3(\frac{x-x_i}{{ϵ}_i}{)}^2}\right),\left(\genfrac{}{}{0ex}{}{(\frac{x-x_i}{{ϵ}_i}{)}^2}{-2\frac{x-x_i}{{ϵ}_i}\frac{y-y_i}{{ϵ}_i}}\right),\left(\genfrac{}{}{0ex}{}{2\frac{x-x_i}{{ϵ}_i}\frac{y-y_i}{{ϵ}_i}}{-(\frac{y-y_i}{{ϵ}_i}{)}^2}\right),\left(\genfrac{}{}{0ex}{}{3(\frac{y-y_i}{{ϵ}_i}{)}^2}{0}\right)\right].$$

$$\begin{array}{rl}& {\mathbf{u}}_h({\mathbf{x}}_i)=({\mathbf{P}}_i^{div}{)}^{⊺}({\mathbf{x}}_i){\mathbf{c}}_i^{\ast },\text{with}{\mathbf{c}}_i^{\ast }={\mathbf{M}}_i^{-1}\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}({\mathbf{P}}_i^{div}{)}^{⊺}({\mathbf{x}}_j)\cdot {\mathbf{u}}_j{W}_{ij},{\mathbf{M}}_i=\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}{\mathbf{P}}_i^{div}({\mathbf{x}}_j)W_{ij}({\mathbf{P}}_i^{div}{)}^{⊺}({\mathbf{x}}_j),\\ & {\nabla }_h{\mathbf{u}}_h({\mathbf{x}}_i)=(\nabla {\mathbf{P}}_i^{div}{)}^{⊺}({\mathbf{x}}_i){\mathbf{c}}_i^{\ast },(\nabla \times \nabla \times {)}_h{\mathbf{u}}_h({\mathbf{x}}_i)=(\nabla \times \nabla \times {\mathbf{P}}_i^{div}{)}^{⊺}({\mathbf{x}}_i){\mathbf{c}}_i^{\ast }.\end{array}$$

$${\mathbf{P}}_i(\mathbf{x})={\left[\frac{1}{\alpha !\beta !}{\left(\frac{x-x_i}{{ϵ}_i}\right)}^{\alpha }{\left(\frac{y-y_i}{{ϵ}_i}\right)}^{\beta }\right]}_{|\alpha +\beta |\le m}.$$

$${\mathbf{c}}_i^{\ast }=\underset{{\mathbf{c}}_i}{\mathrm{argmin}}\left\{\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}{\left[{\kappa }_{ij}({\varphi }_j-{\varphi }_i)-{\mathbf{P}}_i^{⊺}({\mathbf{x}}_{ij}){\mathbf{c}}_i\right]}^2{W}_{ij}\right\},$$

$$\begin{array}{rl}& {\mathbf{c}}_i^{\ast }={\mathbf{M}}_i^{-1}\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}{\mathbf{P}}_i({\mathbf{x}}_{ij})W_{ij}{\kappa }_{ij}({\varphi }_j-{\varphi }_i),\text{with}{\mathbf{M}}_i=\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}{\mathbf{P}}_i({\mathbf{x}}_{ij})W_{ij}{\mathbf{P}}_i^{⊺}({\mathbf{x}}_{ij}),\\ & {\nabla }_h{\varphi }_h({\mathbf{x}}_i)=\frac{1}{2}(\nabla {\mathbf{P}}_i{)}^{⊺}({\mathbf{x}}_i){\mathbf{c}}_i^{\ast }{\nabla }_h^2{\varphi }_h({\mathbf{x}}_i)=\frac{1}{4}({\nabla }^2{\mathbf{P}}_i{)}^{⊺}({\mathbf{x}}_i){\mathbf{c}}_i^{\ast }.\end{array}$$

$${\mathbf{n}}_i\cdot [\frac{1}{2}(\nabla {\mathbf{P}}_i{)}^{⊺}({\mathbf{x}}_i){\mathbf{c}}_i^{\ast }]=g_i,$$

$${\int }_{{\Gamma }_n}\mathbf{\sigma }\cdot d\mathbf{A}=\sum _{i\in {\Gamma }_n}({\mathbf{\sigma }}_i\cdot {\mathbf{n}}_i){\mathcal{A}}_i.$$

$${\int }_{{\Gamma }_n}\mathbf{\sigma }\cdot d\mathbf{A}=\sum _{i\in {\Gamma }_n}\left(p_i\mathbf{I}+\nu \left[{\nabla }_h{\mathbf{u}}_{h,i}+({\nabla }_h{\mathbf{u}}_{h,i}{)}^{⊺}\right]\right)\cdot {\mathbf{n}}_i{\mathcal{A}}_i,$$

$$\Vert {\mathbf{e}}_i^t{\Vert }^2=\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}\Vert \nabla {\mathbf{u}}_j-{\nabla }_h{\mathbf{u}}_{h,i\to j}{\Vert }^2.$$

$${\nabla }_h{\mathbf{u}}_{h,i\to j}=(\nabla {\mathbf{P}}_i^{div}{)}^{⊺}({\mathbf{x}}_j){\mathbf{c}}_i^{\ast }.$$

$$\Vert {\mathbf{e}}_i^r{\Vert }^2=\sum _{j\in {\mathcal{N}}_{{ϵ}_i}}\Vert \mathbf{R}[\nabla \mathbf{u}{]}_j-{\nabla }_h{\mathbf{u}}_{h,i\to j}{\Vert }^2.$$

$$\mathbf{R}[\nabla \mathbf{u}{]}_j=\frac{1}{N_j}\sum _{k\in {\mathcal{N}}_{{ϵ}_j}}{\nabla }_h{\mathbf{u}}_{h,k\to j}=\frac{1}{N_j}\sum _{k\in {\mathcal{N}}_{{ϵ}_j}}(\nabla {\mathbf{P}}_k^{div}{)}^{⊺}({\mathbf{x}}_j){\mathbf{c}}_k^{\ast },$$

$$\text{SOLVE}\to \text{ESTIMATE}\to \text{MARK}\to \text{REFINE}$$

$$\left[\begin{array}{cc}\mathbf{K}& \mathbf{G}\\ \mathbf{B}& \mathbf{L}\end{array}\right]\left[\begin{array}{c}\tilde{\mathbf{u}}\\ p\end{array}\right]=\left[\begin{array}{c}{\mathbf{f}}_{tot}\\ g\end{array}\right].$$

$${\left[\begin{array}{cc}\mathbf{S}& \mathbf{G}\\ 0& \mathbf{L}\end{array}\right]}^{-1},$$

$$E_{tot}\le TOL,$$

$$\kappa =\{\begin{array}{ll}30& \text{for}\text{x}\in {\Omega }_{f,1}\\ 1& \text{for}\text{x}\in {\Omega }_{f,2}\end{array}.$$

Full text

1 Introduction

Fluid-structure interactions (FSI) problems are ubiquitous in science and engineering applications. Examples include, but are not limited to, complex fluids [1, 2, 3, 4, 5, 6], additive manufacturing [7, 8, 9], and energy harvesting [10, 11, 12, 13]. Two traits shared by these problems are: i) moving fluid-solid interfaces of complex geometries; and ii) hydrodynamic interactions between many bodies. Thus, numerical methods for solving such problems must be able to efficiently track the moving interfaces and to accurately capture bidirectional fluid-solid hydrodynamic couplings, which can be challenging particularly when the solids display large displacements and rotations, and/or the moving solids have complex and feature-rich geometries. Moreover, in the zero Reynolds-number limit, solutions to the Stokes flow are well-known to exhibit singular pressures in the vicinity of boundaries with sharp corners and in narrow lubrication gaps occurring between solid bodies. In the flows of colloid suspensions, colloid kinematics of dense suspensions at zero-Reynolds number is dominated by the so-called lubrication effects. Thus, the accurate prediction of these flows depends on the accurate resolution of the singularities over a length scale potentially orders of magnitude smaller than the characteristic length of flow. Such suspension flows therefore provide a particularly interesting FSI application to apply adaptive schemes.

Meshless methods provide attractive and flexible discretization for FSI problems, due to their ability to circumvent costly management of deforming meshes. This is important not only in handling large deformations relevant to FSI, but also in handling the insertion and deletion of nodes in $h$-adaptive schemes without the need to preserve mesh quality. In this work, we consider the generalized moving least squares (GMLS), a recently developed meshless method that has been shown to provide stable high-order solution to the div-grad [14] and steady Stokes [15] problems. GMLS provides a flexible optimization framework to obtain optimal recovery of linear target functionals over a desired reproduction space directly from scattered samples (GMLS particles), and is built upon a rigorous approximation theory [16, 17]. In the div-grad case, stability can be obtained by constructing a staggered discretization analogous to mesh-based primal/dual discretizations; the flexibility of GMLS enables the $\epsilon$-ball graph of neighbor connectivity to be used as a surrogate for primal/dual meshes in, e.g., staggered finite volume methods. For the Stokes problem, GMLS enables the reconstruction of velocity over a space of divergence-free polynomials, which can be efficiently performed locally in a least-squares sense, contrary to mesh-based methods necessitating costly global construction of divergence-free shape functions. The spatial compatibility of these approaches in GMLS stands in contrast to our previous work using smoothed particle hydrodynamics (SPH) [18, 19, 20, 21, 22, 23], a different meshless method in which to date researchers have been unable to simultaneously achieve spatial compatibility and consistency.

As with traditional mesh-based methods utilizing quasi-uniform triangulations, GMLS suffers from deteriorating accuracy in the vicinity of corner singularities and discontinuities in material properties. Additionally, successfully resolving lubrication gaps in Stokes problems while maintaining a uniform resolution imposes a restrictive computational cost. To this end, we propose the introduction of spatial adaptivity into the GMLS discretization.

An appropriate adaptivity criterion must be selected to ensure recovery of optimal convergence in the presence of singularities. Here, we chose to consider an error estimator. In the literature of adaptive finite element methods (AFEM) for solving second-order elliptic partial differential equations (PDEs), there are mainly two types of a posteriori error estimators – the residual-based [24, 25] and recovery-based [26, 27]. A residual-based error estimator is derived from the variational formulation, and hence, strongly depends on the problem being solved [25] while remaining inappropriate for the strong-form discretization used in GMLS. The recovery-based error estimator is solely based upon the computed solution and does not depend on the specific governing PDEs [27]. Thus, a recovery-based error estimator is more applicable to the current scheme than the residual type error estimators. Moreover, although the theoretical justification of the recovery-based error estimator is still under investigation [28, 29], in practice, it has been shown to provide more accurate estimation of true errors when compared to the residual-based error estimators  [30, 28]. Therefore, we pursued a recovery-based error estimator as the criterion to direct spatial adaptivity in the GMLS discretization.

Given the adaptive criterion, we proposed a four-step adaptive algorithm: $\text{SOLVE}\rightarrow\text{ESTIMATE}\rightarrow\text{MARK}\rightarrow\text{REFINE}$, to adaptively refine the spatial resolutions. In the SOLVE step, the governing equations are numerically solved using the GMLS discretization, and all of the variables are updated. In the ESTIMATE step, the recovery-based a posteriori error estimator is evaluated at all GMLS particles, which measures the difference between the direct and recovered gradient. Here, the recovered gradient at a GMLS particle is defined as the average of the direct solutions of the gradient (e.g., velocity gradient) obtained on its neighbor particles. In the MARK step, the GMLS particles with larger recovery-based errors are marked for refining, following a strategy similar to the Dörfler’s strategy used in AFEM [31, 32, 25]. Finally, in the REFINE step, the GMLS particles marked are refined by splitting them into smaller particles. At each time step, these four steps are iterated until the total estimated error of all GMLS particles is no more than a preset tolerance.

The proposed spatially adaptive GMLS method was first validated and assessed by solving the div-grad problems with corner singularity or subject to discontinuous coefficients. By comparing with the analytical solutions, we examined whether the adaptive GMLS solutions recover the optimal, theoretical convergence orders and how accurate the proposed recovery-based error provides estimations of the true errors. For GMLS, only a truncation error analysis for differential operators exists [17], and the lack of formal error analysis renders the notion of optimal convergence ill-defined. In this work, we informally refer to optimal convergence as matching the convergence observed in previous work for non-singular solutions and regular domains [14, 15]. Next, we employed the adaptive GMLS for solving Stokes flows subject to bidirectional fluid-solid couplings with the FSI applications of colloid suspensions. By comparing the numerical results with the analytical solutions and experimental data, we examined the accuracy of the proposed adaptive GMLS method in solving the FSI problems. By comparing with the GMLS solutions with uniform resolutions, we demonstrated the efficiency gains of the adaptive GMLS to achieve the same accuracy.

This paper is organized as follows. Section 2 presents the governing equations to be solved for the div-grad and steady Stokes problems, respectively. We provide in Section 3 a concise statement of the GMLS discretization employed in this work. Section 4 explains the proposed recovery-based error estimator and algorithm for adaptive refinement. In Section 5, we validated the adaptive GMLS through numerical tests. We solved two div-grad problems and four Stokes flows with applications to colloid suspensions – the div-grad problem in an L-shaped domain, div-grad problem with rough coefficients, Wannier flow, dynamics of an L-shaped active colloid, dynamics of two square colloids under a shear flow, and collective dynamics of active asymmetric colloidal dimers. In this exercise, the numerical results were compared with the analytical solutions, GMLS solutions with uniform resolutions, or experimental data. Finally, we conclude and suggest directions of future work in Section 6.

2 Governing equations

The proposed method was employed for solving both diffusion (div-grad) and Stokes problems. The computational domain can be represented as $\Omega=\Omega_{f}\cup\Gamma$, where $\Omega_{f}$ is the sub-domain occupied by the fluid, and $\Gamma$ denotes the boundary.

2.1 Div-grad problem

The div-grad boundary value problem is governed by:

[TABLE]

where $\kappa\in\mathbf{L^{2}}(\Omega)$ describes a potentially discontinuous diffusion coefficient; all $\kappa$, $f$, and $\phi_{\Gamma}$ are given data.

2.2 Steady Stokes problem

For the Stokes problem, we consider the coupled dynamics of freely suspended colloids and steady Stokes flow of incompressible fluid. Assume the domain $\Omega$ contains $N_{c}$ colloids, and each of them with a boundary $\Gamma_{n}$, $n=1,\dots,N_{c}$. Assuming finite separation between any two colloids, the boundary of the domain may be partitioned into the disjoint union $\Gamma=\Gamma_{w}\cup\Gamma_{1}\cdots\cup\Gamma_{N_{c}}$, where $\Gamma_{w}$ denotes the wall boundaries. Each colloid has a position $\mathbf{X}_{n}$, an orientation $\mathbf{\Theta}_{n}$, and undergoes rigid-body kinematics prescribed by a translational velocity $\dot{\mathbf{X}}_{n}$ and rotational velocity $\dot{\mathbf{\Theta}}_{n}$. We refer to the collection of all colloids as a suspension. For a given $\left\{\mathbf{X}_{n},\mathbf{\Theta}_{n}\right\}_{n=1,\dots,N_{c}}$, the geometrical configuration of the suspension (i.e., the set $\left\{\Gamma_{n}\right\}_{n=1,\dots,N_{c}}$) is prescribed, and the fluid flow is coupled to the colloid motion via the following steady Stokes problem, providing a boundary value problem for the fluid velocity and pressure:

[TABLE]

where $\nu$ is the kinematic viscosity of fluid; $\mathbf{f}$ denotes a body force such as gravity exerted on the fluid; $\mathbf{w}$ is the velocity of the wall boundaries, and $\mathbf{w}=\mathbf{0}$ for stationary walls.

For a divergence-free velocity field, the vector identity $\nabla^{2}\mathbf{u}=-\nabla\times\nabla\times\mathbf{u}$ holds. Thus, Eq. (2) may be recast in the following equivalent form:

[TABLE]

Here, $\mathbf{n}$ is the unit normal vector outward facing at boundary $\Gamma$. By recasting the Stokes equation (2) in this equivalent form, the saddle-point structure of the Stokes operator can be avoided, and after discretization one obtains instead a system with elliptic matrices along the diagonal blocks that is more amenable to standard preconditioning techniques. This necessitates a discretization faithful to the divergence-free constraint. In the GMLS discretization, the velocity is directly constructed from an appropriate space of divergence-free vectors, as explained in Section 3.

In general, the translational and angular dynamics of each colloid is governed by:

[TABLE]

where $M_{n}$ and $I_{n}$ denote the mass and moment of inertia of each colloid, respectively; $\bm{\sigma}=-p\mathbf{I}+{\nu}\left[\nabla\mathbf{u}+(\nabla\mathbf{u})^{\intercal}\right]$ is the stress exerted by the fluid on colloids; $\mathbf{F}_{e,n}$ and $\mathbf{T}_{e,n}$ represent the external force and torque applied on each colloid, respectively. Assuming the inertia effect is negligible, which is valid for small colloids in a viscous fluid, the above dynamics of each colloid is subject to the force- and torque-free constraint; i.e.,

[TABLE]

Therefore, by solving Eqs. (3) and (5) concurrently as a monolithic system, we obtained a boundary value problem for the colloids’ translational and angular velocities $\left\{\dot{\mathbf{X}}_{n},\dot{\mathbf{\Theta}}_{n}\right\}_{n=1,\dots,N_{c}}$ as well as the fluid’s velocity and pressure fields. Given $\left\{\dot{\mathbf{X}}_{n},\dot{\mathbf{\Theta}}_{n}\right\}_{n=1,\dots,N_{c}}$, the colloids’ positions and orientations $\left\{\mathbf{X}_{n},\mathbf{\Theta}_{n}\right\}_{n=1,\dots,N_{c}}$ were numerically evolved by the 4th-order Runge-Kutta time-integration scheme [33].

3 GMLS discretization

To numerically solve the governing equations given in Section 2, separate treatment was applied to the velocity $\mathbf{u}$ and to the scalar fields $\phi$ and $p$ appearing in Eqs. (1) and (3), respectively. While the velocity $\mathbf{u}$ in the Stokes equations was handled via a GMLS approximation reconstructed over the space of divergence-free vector polynomials, $\phi$ in the diffusion equation and the pressure $p$ in the Stokes equations were approximated via a staggered GMLS discretization. In this section, we briefly review the basics of both the divergence-free [15] and staggered [14] approaches employed in this work. Regarding their derivation and analysis, we refer to the previous work [14, 15] for details.

3.1 GMLS approximation

The fluid domain $\Omega$ is discretized by a set of collocation points (called GMLS particles); another set of GMLS particles are distributed along the colloid and wall boundaries $\Gamma_{w}\cup\Gamma_{1}\cdots\cup\Gamma_{N_{c}}$ for imposing the boundary conditions (BCs). We refer to these as interior and boundary particles, respectively. For a GMLS particle at $\mathbf{x}_{i}$ and a given scalar function $\psi$ evaluated at its neighbor locations: $\psi_{j}=\psi(\mathbf{x}_{j})$, a polynomial $\psi_{h}(\mathbf{x})$ of order $m$ is sought to approximate $\psi$ and its derivatives $D^{\alpha}\psi$ at $\mathbf{x}_{i}$. To this end, $\psi_{h}(\mathbf{x})=\mathbf{P}^{\intercal}(\mathbf{x})\mathbf{c}^{*}$ with a polynomial basis $\mathbf{P}(\mathbf{x})$ and coefficient vector $\mathbf{c}^{*}$ such that the following weighted residual functional is minimized:

[TABLE]

As such, for $\psi\in span(\mathbf{P}(\mathbf{x}))$, $\psi$ can be exactly reconstructed. From this polynomial reconstruction property, it follows that high-order accuracy can be achieved via the GMLS approximation by taking large $m$, e.g., $m=4$ [16]. Following standard arguments for the minimization of a symmetric positive definite quadratic form, the solution of Eq. (6) is given by:

[TABLE]

To approximate derivatives of the underlying function, an arbitrary $\alpha^{th}$ order differential operator can be written as [17]:

[TABLE]

In Eqs. (6) and (7), the weight function $W_{ij}=W(r_{ij})$, where $r_{ij}=\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|$; $W(r)=1-(\frac{r}{\epsilon})^{4}$ for $r<\epsilon$, otherwise, $W(r)\equiv 0$, with $\epsilon$ the compact support. Hence, it is only necessary to include GMLS particles within an $\epsilon$-neighborhood of the $i$th GMLS particle, i.e., $j\in{\cal{N}}_{\epsilon_{i}}=\left\{\textbf{x}_{j}~{}\text{s.t.}~{}\left\|\textbf{x}_{i}-\textbf{x}_{j}\right\|<\epsilon_{i}\right\}$. As a result, the approximation error in both the interpolant and derivatives of $\psi$ is governed by the choice of $\epsilon$. While a smaller $\epsilon$ yields higher accuracy and sparser linear operators, $\epsilon$ must be chosen sufficiently large to ensure unisolvency over the reconstruction space, and thus, a solution to the quadratic program [16]. In this work, $\epsilon_{i}=4.5\Delta{x}_{i}$ with $\Delta{x}_{i}=V_{i}^{1/d}$, where $V_{i}$ is the volume occupied by particle $i$, and $d$ denotes the dimension, e.g., $d=2$ in 2D. For the adaptive strategy considered in this work, particles can have different $\Delta{x}$, and in turn, different $\epsilon$. We therefore require a definition of $W_{ij}$ whose support adapts to the local particle density in the vicinity of $\mathbf{x}_{i}$. To ensure a symmetric weighting ($W_{ij}=W_{ji}$), we define the weight function as:

[TABLE]

3.2 Divergence-free GMLS reconstruction for the velocity $\mathbf{u}$

Based on the GMLS approximation discussed above, we seek a polynomial reconstruction $\mathbf{u}_{h}(\mathbf{x})=(\mathbf{P}^{div})^{\intercal}(\mathbf{x})\mathbf{c}^{*}$ to approximate $\mathbf{u}$ and to discretize its gradient and curl-curl operators. To enforce compatibility with the divergence-free constraint for the velocity, we chose the polynomial basis $\mathbf{P}^{div}(\mathbf{x})$ from the space of $m^{th}$ order divergence-free vector polynomials. For example, in 2D and for $m=2$, a divergence-free vector polynomial basis $\mathbf{P}^{div}(\mathbf{x})$ at particle $i$ is given by:

[TABLE]

The vectorial extension of Eqs. (7) and (8), providing expressions for the reconstruction, gradient, and curl-curl operators is given by [15]:

[TABLE]

The Dirichlet BC in Eq. (3) for the velocity was then enforced by assigning the values of $\mathbf{u}$ as specified in Eq. (3) over the set of boundary GMLS particles.

3.3 Staggered discretization of div-grad operator

To ensure a compatible discretization of the div-grad operator for $\phi$ in the diffusion equation and for the pressure $p$ in the Stokes equations, a staggered GMLS discretization [14] was employed. Motivated by compatible mesh-based discretization on primal-dual grids [34, 35], a local primal-dual complex was constructed for each GMLS particle via an $\epsilon$-neighborhood graph (virtual cell). In each virtual cell, a set of primal edges can be built as: $\mathbf{E}_{i}=\{\mathbf{x}_{i}-\mathbf{x}_{j}~{}|~{}\mathbf{x}_{j}\in{\cal{N}}_{\epsilon_{i}}\}$; and, each edge is associated with a midpoint $\mathbf{x}_{ij}=\frac{\mathbf{x}_{i}+\mathbf{x}_{j}}{2}$ and a virtual dual face at the midpoint and normal to the edge. The staggered GMLS discretization treats the div-grad operator as a composition of a GMLS divergence recovered from local virtual dual faces together with a topological gradient over primal edges, providing a generalization of familiar staggered finite difference methods to unstructured point clouds [14]. We succinctly describe below the staggered GMLS discretization for $\phi$, but it also applies to $p$. A detailed derivation and analysis of the discretization may be found in [14].

Select as reconstruction basis $\mathbf{P}$ the $\epsilon_{i}$-scaled Taylor monomials, which may be expressed in 2D as:

[TABLE]

Define the polynomial $q_{i}(\mathbf{x})=\mathbf{P}^{\intercal}_{i}(\mathbf{x})\mathbf{c}^{*}_{i}$, with coefficients determined by solving

[TABLE]

where $\kappa_{ij}=\frac{\kappa(\mathbf{x}_{i})+\kappa(\mathbf{x}_{j})}{2}$. The optimal coefficient vector, gradient, and Laplacian are then provided by:

[TABLE]

Evaluation of the gradient and Laplacian operators at interior GMLS particles follows from a direct evaluation of Eq. (14). For boundary GMLS particles, either a Dirichlet BC must be imposed in Eq. (1), or an inhomogeneous Neumann BC must be handled in Eq. (3). The Dirichlet BC was imposed by simply setting: $\phi_{i}=\phi_{\Gamma}(\mathbf{x}_{i})$. To impose an inhomogeneous Neumann BC of the form $\partial_{\mathbf{n}}|_{\Gamma}=g$, we added a single equality constraint to the quadratic program Eq. (13) to enforce

[TABLE]

where $\mathbf{n}_{i}$ is the unit outward facing normal at boundary particle $i$. Further details regarding the imposition of Neumann BCs may be found in [14].

3.4 Approximation of integrals at colloid boundaries

The integral of the stress tensor $\bm{\sigma}$ in Eq. (5) requires the construction of a meshless quadrature rule. Consider the collection of boundary particles $\mathbf{x}_{i}\in\Gamma_{n}$. We assumed the boundary may be partitioned into a disjoint collection of planar/curved faces $\Gamma_{ni}$. For each $\Gamma_{ni}$, we associated a unit outward facing normal $\mathbf{n}_{i}$ and measure of face area $\mathcal{A}_{i}$. A composite quadrature rule for the integral may then be prescribed as:

[TABLE]

The stress tensor $\bm{\sigma}=-p\mathbf{I}+{\nu}\left[\nabla\mathbf{u}+(\nabla\mathbf{u})^{\intercal}\right]$ must be reconstructed at $\mathbf{x}_{i}$, where $p$ is already available, and it only remains to approximate the gradient of the velocity $\nabla\mathbf{u}$. Following the formulations in Section 3.2, we approximated $\nabla\mathbf{u}$ at each boundary particle. The integral was then approximated as:

[TABLE]

where $\nabla_{h}\mathbf{u}_{h,i}=\nabla_{h}\mathbf{u}_{h}(\mathbf{x}_{i})=(\nabla\mathbf{P}^{div}_{i})^{\intercal}(\mathbf{x}_{i})\mathbf{c}_{i}^{*}$. For the purposes of this work, this quadrature rule sufficed to obtain the reported high-order convergence results, while the interested reader may consider referring to [15] for a discussion of higher-order quadrature rules.

4 Algorithm for adaptive refinement

In this section, we introduce the algorithm for achieving adaptive refinement in the GMLS discretization. The algorithm relies on a recovery-based a posteriori error estimator as the adaptive criterion.

4.1 Recovery-based error estimator

For the div-grad problem, the error estimator was evaluated for $\phi$; for the Stokes problem, the errors were estimated based on the velocity $\mathbf{u}$. We describe below the proposed recovery-based error estimator for $\mathbf{u}$, but it also applies to $\phi$.

Assuming the exact solution of the velocity gradient $\nabla\mathbf{u}$ and the direct solution $\nabla_{h}\mathbf{u}_{h}$, the true error of the GMLS approximation at particle $i$ is:

[TABLE]

Here, the energy norm $\|\cdot\|^{2}$ is defined as the element-wise inner product; i.e., $\|\mathbf{e}\|^{2}=\mathbf{e}:\mathbf{e}=e^{mn}e^{mn}$, where the Einstein summation convention is adopted.

The direct solution $\nabla_{h}\mathbf{u}_{h,i\rightarrow j}$ in Eq. (18) was computed by reconstructing the velocity gradient at particle $i$ and evaluating it at particle $j$; i.e.,

[TABLE]

Thus, if $j$ was considered the neighbor of a different particle $i^{\prime}$, although still evaluated at particle $j$, $\nabla_{h}\mathbf{u}_{h,i^{\prime}\rightarrow j}$ would be different. The definition of $\nabla_{h}\mathbf{u}_{h,i\rightarrow j}$ is further illustrated in Figure 1, where the curved surfaces with different colors represent the velocity gradient reconstructed at different particles: $\nabla_{h}\mathbf{u}_{h}=(\nabla\mathbf{P}^{div}_{i})^{\intercal}(\mathbf{x})\mathbf{c}_{i}^{*}$, $i=$ 1, 2, 3, and 4. Particle 5 is a neighbor particle shared by particles 1, 2, 3, and 4. Evaluating $\nabla_{h}\mathbf{u}_{h,i\rightarrow j}$ at particle 5 ($j=5$) reads four different values.

The true error as defined in Eq. (18) cannot always be evaluated due to a lack of the exact solution in practical applications. Instead, the recovered error, which measures the difference between the direct and recovered velocity gradient [36], can be practically determined in any simulation. In the GMLS, this recovered error can be evaluated locally on each GMLS particle at every time step as:

[TABLE]

Here, $\mathbf{R}[\nabla\mathbf{u}]$ represents the recovered velocity gradient, and we proposed the following definition for it:

[TABLE]

where $k$ and $N_{j}$ denote the index and total number of the neighbor particles, respectively, of particle $j$. Thus, the recovered velocity gradient at particle 5 in Figure 1 can be defined as: $\mathbf{R}[\nabla\mathbf{u}]_{5}=\frac{1}{4}(\nabla_{h}\mathbf{u}_{h,1\rightarrow 5}+\nabla_{h}\mathbf{u}_{h,2\rightarrow 5}+\nabla_{h}\mathbf{u}_{h,3\rightarrow 5}+\nabla_{h}\mathbf{u}_{h,4\rightarrow 5})$, assuming it has only four neighbors.

Following the AFEM [30, 29], being an appropriate adaptive criterion, the recovered error should have the same convergence order as the true error; i.e., the a posteriori error estimator is reliable and robust. In Section 5, we numerically assessed and validated the proposed recovery-based error estimator given by Eqs. (20) and (21).

4.2 Adaptive algorithm

Given the recovery-based a posteriori error estimator as the adaptive criterion, we introduce the following adaptive algorithm. At each time step, starting from a uniform, coarse resolution, the spatial resolution was adaptively refined following a four-step iteration procedure:

[TABLE]

In the SOLVE step, the governing equations, Eq. (1) for the div-grad problem or Eqs. (3) and (5) for the Stokes problem, were numerically solved by GMLS, and all of the variables were updated. In particular, to efficiently solve the resulting linear system after the GMLS discretization, we employed different linear solvers for different equations. Specifically, for Eq. (1), since it is a Poisson problem, we used the algebraic multigrid (AMG) method. For Eqs. (3) and (5), the resulting linear system has the following block structure:

[TABLE]

Here, we combined the velocities at all interior and boundary GMLS particles and colloids’ translational and angular velocities into the vector $\tilde{\mathbf{u}}$. $\mathbf{K}$ contains contributions from the curl-curl ($\nabla\times\nabla\times$) operator in Eq. (3) and the velocity gradient in the stress $\bm{\sigma}$ in the force- and torque-free constraint (5). $\mathbf{G}$ encompasses the $\frac{1}{\rho}\nabla$ operator in Eq. (3) and the contribution from the pressure in the stress $\bm{\sigma}$ in the force- and torque-free constraint (5). $\mathbf{L}$ corresponds to the Laplacian operator ($\nabla^{2}$) in Eq. (3). $\mathbf{B}$ is non-zero only at the boundaries and contains the contribution from the $\nu\mathbf{n}\cdot\nabla\times\nabla\times$ operator in the inhomogeneous Neumann BC in Eq. (3). On the right-hand side of Eq. (23), $\mathbf{f}_{tot}$ combines the body force $\mathbf{f}$, $\mathbf{F}_{e,n}$ and $\mathbf{T}_{e,n}$ applied on each colloid, and the velocity $\mathbf{w}$ of the wall boundaries; $g$ contains $\nabla\cdot\mathbf{f}$ and $\mathbf{n}\cdot\mathbf{f}$ in Eq. (3). Following the previous work [15], we used the preconditioned general minimal residual (PGMRes) method to solve the above linear system with a block triangular preconditioner given as follows:

[TABLE]

where $\mathbf{S}=\mathbf{K}-\mathbf{G}\,\text{diag}(\mathbf{L})^{-1}\,\mathbf{B}$ is an approximate Schur complement. To employ this preconditioner, both $\mathbf{S}$ and $\mathbf{L}$ needed to be inverted. In our implementation, we used the AMG method for the inversion. Due to the singularity of the solutions and adaptive refinement, standard relaxation schemes, such as Gauss-Seidel smoother, were not effective in the AMG. Thus, we applied a drop-tolerance type incomplete LU (ILU) smoother [37] on the finest level of the AMG to improve the robustness of AMG both for solving the linear system of the div-grad problem and inverting the matrices in the block preconditioner (24) for the Stokes problem.

In the ESTIMATE step, the a posteriori recovered error was firstly estimated locally at each GMLS particle based on Eqs. (20) and (21). Next, the particles were ranked according to the magnitude of the estimated recovered error. Finally, the total error of all GMLS particles: $E_{tot}=\sum\limits_{i}\|\mathbf{e}^{r}_{i}\|^{2}V_{i}$, was evaluated.

In the MARK step, a portion of GMLS particles were selected and marked for refining according to the local, recovered error estimated in the ESTIMATE step. The particles with larger estimated errors are the potential candidates for refinement. In practice, the top $N_{t}$ particles ranked in the ESTIMATE step with a total error of $\alpha_{s}E_{tot}$ were marked. This MARK strategy is similar to the Dörfler’s strategy used in the AFEM [31, 32, 25]. Here, $\alpha_{s}\in(0,1)\times 100\%$ is the preset marking percentage. The value of $\alpha_{s}$ should be chosen depending on the specific problem and according to the distribution of the estimated recovered errors. A higher marking percentage is recommended if the estimated errors vary by orders across different regions of the domain. Otherwise, a smaller $\alpha_{s}$ is sufficient. Note that a larger $\alpha_{s}$ calls for more particles marked and refined, thereby a more expensive simulation.

In the REFINE step, the GMLS particles marked in the MARK step were refined. In particular, a marked interior GMLS particle was split into 4 small particles (in 2D); a marked boundary particle was split into 2 small particles, as illustrated in Figure 2. The spacing length $\Delta X_{H}$ of the new, small particles is half of $\Delta X_{L}$ of the original, big particles. Hence, the compact support $\epsilon$ of the new, small particles is also half of that of the original, big particles; their volume $V=\Delta x^{2}$ is $\frac{1}{4}$ of the original, big particles. There are two caveats herein. First, if a boundary has a sharp corner, we always put a particle (black in Figure 2) at the corner to ensure accurate representation of the boundary. This corner particle maintained its compact support $\epsilon$ and volume $V$ the same as those of its two neighbor boundary particles. That says, if its neighbor boundary particles were refined, this corner particle synchronized its $\epsilon$ and $V$ with them. Second, due to the meshless nature of the GMLS discretization, the interior GMLS particles are allowed not perfectly aligned with the boundaries. After several times of refinement, there could occur near a boundary voids of size larger than the average spacing of surrounding interior particles, which would result in larger local numerical errors. To remedy it, if any interior particle near a boundary needed to be refined, we identified the images of this particle and its neighbor interior particles on the opposite side of the boundary and refined those image (imaginary) particles as well. The resulting new, small particles within the solution domain were kept, which could effectively fill the voids.

The above four steps (22) were repeated iteratively at each time step until the total error satisfied:

[TABLE]

where $TOL$ is the preset tolerance.

5 Simulation results

Two div-grad problems and four Stokes flow problems were studied to assess the accuracy and convergence of the proposed spatially adaptive GMLS method.

5.1 Constant coefficient div-grad problem in an L-shaped domain

We first considered the div-grad problem with $\kappa\equiv 1$ in Eq. (1) on an L-shaped domain. This benchmark problem was also studied in a finite element context [38, 25] and isogeometric analysis [39]. The geometry of the solution domain is described in Figure 3. The source term in Eq. (1) was taken as $f=1$. Under polar coordinates, the analytical solution of this problem is [25]: $\phi(r,\theta)=r^{\frac{2}{3}}\sin(\frac{2\theta}{3})-\frac{r^{2}}{4}$.

Because the L-shaped geometry has a corner singularity, the numerical solution of this problem with uniform discretization cannot achieve the theoretical order of convergence, as demonstrated in [25]. Indeed, we identified the same behavior in our simulations performed with a uniform distribution of GMLS particles. As shown in Figure 4, regardless of the polynomial order employed in the GMLS approximation, the convergence order of the numerical solutions remains around 0.67, consistent with reported FEM results [25].

Instead of refining the discretization uniformly, we employed adaptive refinement as proposed in Section 4.2. Figure 5 presents a snapshot of particle distribution under adaptive refinement. As anticipated, the maximum errors occurred near the singular corner of the L-shaped domain. As a result, the GMLS particles around that corner were refined for multiple times to reduce the errors. During these iterations of refinement, the generated adaptive resolution entailed much faster convergence of the numerical solutions, as depicted in Figure 4. Notably, e.g., from Figure 4 (b), to achieve the same accuracy, employing adaptive refinement can significantly reduce the number of GMLS particles needed for solving this problem. Furthermore, the adaptive GMLS solution exhibits second or fourth order convergence for second or forth order polynomials used in the GMLS reconstruction, matching the convergence rates for smooth solutions in the literature [14].

In addition, we note that the estimated errors by Eq. (20) based on the recovered gradient of $\phi$ display the same convergence order as the true errors, in both uniform-resolution and adaptive-resolution scenarios.

5.2 Div-grad problem with rough coefficients

We next solved the div-grad problem with a discontinuous coefficient $\kappa$ in Eq. (1), which describes various physical phenomena, e.g., transport in porous media with discontinuous permeability [40, 41, 42, 43] and electrostatics with discontinuous dielectric permittivity [44, 45, 46]. The solution domain of this problem is a square of dimension $10\times 10$. A circle of radius $2$ positioned at the center divides the square into two sub-domains: inside the circle ($\Omega_{f,1}$) and outside the circle ($\Omega_{f,2}$). The coefficient $\kappa$ in Eq. (1) is discontinuous across the two sub-domains and given as:

[TABLE]

The source term in Eq. (1) was taken as $f=0$. To compare with the analytical solution reported in [45] obtained in an infinite domain, the value of $\phi_{\Gamma}$ in Eq. (1) was imposed on the boundary of the truncated square domain as a Dirichlet BC [45].

Provided such a large coefficient ratio between the two adjacent sub-domains as in Eq. (26), the numerical solution near the interface tends to develop oscillations, and hence, it’s challenging to solve this problem stably and accurately. Trask et al. has previously demonstrated that the compatible, staggered GMLS discretization can effectively eliminate numerical oscillations and produce stable and accurate solutions for this problem [14]. However, we find that due to the discontinuous coefficient, the convergence of the numerical solution is diminished when using uniformly distributed GMLS particles. As depicted in Figure 6, the true errors of the numerical solutions with uniform refinement exhibit less than 1st-order convergence. Here, to avoid the Runge phenomenon [47, 48], the 1st-order polynomial basis was employed in the reconstruction of the gradient, which theoretically entails 1st-order convergence.

Thus, we demonstrate herein that the optimal convergence order can be achieved by introducing adaptive refinement in the staggered GMLS discretization. Figure 7 illustrates the distribution of GMLS particles after several iterations of refinement. The proposed error estimator directed the adaptive refinement with the finest resolution near the interface of the two sub-domains.

Figure 8 presents the computed $\frac{\partial\phi}{\partial x}$ (flux along $x$) along the line ($y=0$) passing through the center of the square domain. The numerical solution agrees well with the analytical solution [45]. The convergence study shown in Figure 6 confirms that the adaptive GMLS solution can achieve the theoretical 1st-order convergence.

The estimated errors based on the recovered gradient of $\phi$ display the same convergence order as the true errors, in both uniform-resolution and adaptive-resolution scenarios.

5.3 Wannier flow

We now move forward to the Stokes problems. In this section, we examined Wannier flow [49], a Stokes flow confined between two rotating cylinders. Wannier flow has its analytic solutions for both pressure and velocity, providing a good benchmark for assessing accuracy and convergence of the proposed method. GMLS was previously employed for solving Wannier flow with the convergence assessed for uniform resolution [15].

The specific setup of the problem is sketched in Figure 9, where the radii of the two cylinders are $r_{1}$ and $r_{2}$, respectively; the distance between the two cylinders’ centers is defined as $\Delta S$; and the two cylinders are rotating in the same direction with the rotational speeds $\omega_{1}$ and $\omega_{2}$, respectively. The fluid was assumed to have the density $\rho=1.0~{}$\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3} and kinematic viscosity $\nu=1.0~{}$\mathrm{m}^{2}\mathrm{/}\mathrm{s}.

Starting with a coarse, uniform distribution of GMLS particles, we again followed the algorithm described in Section 4 to apply adaptive refinement. During the process, the marked GMLS particles were split, and hence, the total number of GMLS particles ($N^{2}$) increased. In each iteration of refinement, the true and estimated recovery-based errors were evaluated, as depicted in Figure 10. It can be seen that the optimal convergence order is obtained for the gradient of velocity, regardless of 2nd or 4th-order polynomial reconstruction employed in GMLS. For comparison, we also employed uniform refinement and evaluated the true and estimated errors, also depicted in Figure 10. Notably, using the same number of GMLS particles, the adaptive refinement yield more accurate solutions than uniform refinement. Given the same acceptable error, employing adaptive refinement would require many fewer degrees of freedom in the simulation. Also, we note that the estimated errors by Eq. (20) based on the recovered gradient of velocity closely match the true errors, in both uniform and adaptive refinement. In this specific test, we set $\Delta S=\frac{\pi}{5}~{}$\mathrm{m}$$.

By adjusting the center position of the inner cylinder, the minimum gap width between the two cylinders ($r_{2}-r_{1}-\Delta S$) varies. As the two cylinders draw close, a numerical solution could fail due to lack of sufficient discretization points within the narrow gap between the two cylinders. In this case, adaptive refinement is necessary to obtain an accurate and computationally tractable numerical solution. Figure 11 depicts the resulting distribution of GMLS particles after 6 iterations of refinement. Here, the minimum gap width between the two cylinders ($r_{2}-r_{1}-\Delta S$) is as small as $\frac{r_{1}}{20}$. From Figure 11, we see that adaptation allows for sufficient discretization points within the narrow gap, thereby enabling a stable and accurate numerical solution in such an extreme case. The computed pressure distribution shown in Figure 11 agrees with the analytical solution [49].

Further, we computed the drag force exerted on the inner cylinder by the fluid at varying gap widths. As shown in Figure 12, the numerical predictions show good agreement with the analytical solution [49] for all cases considered. For the results presented in both Figure 11 and Figure 12, the 4th-order GMLS polynomial reconstruction was employed for both the velocity and pressure.

5.4 Dynamics of an L-shaped active colloid

In this section, we reproduced an interesting phenomenon concerning the motion of an asymmetric self-propelled colloid in fluid, which was first observed experimentally [5]. The colloid is L-shaped with two asymmetric arms, and its motion is constrained in 2D, as in the experiment [5]. In our simulation, the L-shaped colloid was suspended in a square domain chosen large enough that the boundaries had insignificant effects on the colloid dynamics. As the colloid moved, we always ensured that the colloid’s center of mass (COM) coincided with the center of the square domain, as depicted in Figure 13. The BCs specified in Eq. (3) were enforced on the boundaries of the square domain and L-shaped colloid. The colloid was subject to the gravitational force $\mathbf{g}$ with $\|\mathbf{g}\|=0.5\times 10^{-9}~{}$\mathrm{N}. To mimic the colloid’s self-propelling mechanism as in the experiment, the equivalent force $\textbf{F}_{e}$ and torque $\textbf{T}_{e}$ were applied at the colloid’s COM. During the motion of the colloid, $\textbf{F}_{e}$ was exerted always along the longer arm of the colloid. The torque $\textbf{T}_{e}$ was related to $\textbf{F}_{e}$ by $\textbf{T}_{e}=(l-x^{\prime}_{p})\textbf{e}_{\textbf{x}^{\prime}}\times\textbf{F}_{e}$, where $x^{\prime}_{p}$ and $\textbf{e}_{\textbf{x}^{\prime}}$ were defined in the local coordinate system on the colloid with the origin $O^{\prime}$ set at the bottom left corner and $\textbf{e}_{\textbf{x}^{\prime}}$ along the shorter arm. The coordinates of the colloid’s COM $(x^{\prime}_{p},~{}y^{\prime}_{p})=(1.98\times 10^{-6}~{}$\mathrm{m}$,~{}2.02\times 10^{-6}~{}$\mathrm{m}$)$. The fluid assumed the properties of water: $\rho=1.0\times 10^{3}~{}$\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3} and $\nu=1.0\times 10^{-6}~{}$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$$, consistent with the experiment [5].

As found in the experiment [5], by varying the magnitude of $\textbf{F}_{e}$ (and thereby the magnitude of $\textbf{T}_{e}$ as well), the L-shaped active colloid can display rich dynamic behaviors owing to its nonlinear hydrodynamic interaction with the surrounding fluid. As a result, the colloid can follow different types of trajectories with different magnitudes of $\textbf{F}_{e}$. More specifically, for a small $\textbf{F}_{e}$, since the colloid’s self-propelling motion is slow, and the impact of gravitational force is dominant, the colloid exhibits straight downward swimming (SDS). As $\textbf{F}_{e}$ gradually increases, the colloid’s active motion becomes more pronounced, and its hydrodynamic interaction with the surrounding fluid starts to compete with the gravitational force, leading to straight upward swimming (SUS) or even trochoid-like motion (TLM). Because the associated fluid dynamics falls into the low Reynolds number regime [5], it may be treated as a Stokes problem.

To solve the Stokes equations (3) and (5), we employed the proposed adaptive GMLS method. The numerical solutions captured all three typical types of motions, as illustrated in Figure 14. In Figure 15, we further demonstrate two snapshots of the adaptive GMLS simulation, where the distribution of GMLS particles with adaptive refinement can be seen. The refinement directed by the error estimator concentrated around the corners of the L-shape, which should be anticipated. The two snapshots correspond to the two points ($P_{1}$ and $P_{2}$) marked on the TLM trajectory shown in Figure 14, respectively. (Refer to the animation in the Supporting Information for the entire dynamics of TLM.) In Figure 15, we also present the computed pressure field surrounding the colloid. To verify the adaptive GMLS solution, we compared the computed trajectory of the colloid with that predicted by GMLS with a uniform, high resolution. In particular, the most complicated TLM was chosen for verification. As depicted in Figure 14, the two solutions overlap with each other. Notably in this test, while an average of 13800 GMLS particles were used in the adaptive-resolution GMLS, 58810 GMLS particles were needed in the simulation with the uniform, high resolution. Obviously, the adaptive GMLS required much fewer degrees of freedom for the same accuracy.

5.5 Dynamics of two square colloids under a shear flow

In this section, we solved the Stokes problem with two colloids immersed in fluid. In particular, the two colloids were square-shaped and subject to a shear flow. Figure 16 illustrates the setup of this problem. Two colloids of equal size were suspended in a square domain large enough to neglect wall effects. The fluid had the same density and kinematic viscosity as in Section 5.3. To generate a plane shear flow, the top and bottom walls were assigned velocities consistent with the target shear rate. Following the shear flow, the two colloids rotate while approaching each other and then moving apart. When the two colloids get very close, it can be challenging for a numerical solver to be stable and still provides accurate solutions without invoking a subgrid-scale lubrication model [50, 51, 52].

In literature, a similar problem but with two circular colloids has been solved both analytically [53] and numerically [54, 15, 23]. Solving the problem with square colloids is more challenging due to the presence of sharp corners causing pressure singularities. To this end, the GMLS discretization enables a sharp-representation of arbitrary shapes, and the adaptive GMLS allows for resolving the singularity with sufficient resolution to ensure stable and accurate solutions. To validate, we also solved the problem with two circular colloids and compare with its analytical solution [53]. The 4th-order GMLS polynomial reconstruction was employed for both the velocity and pressure.

Figure 17 depicts the trajectories computed by adaptive GMLS. Due to the symmetry of the two colloids, we only tracked the trajectory of the COM of the colloid initially on the left. As shown in Figure 17, our numerical prediction well agrees with the analytical solution [53] for the case of circular colloids, which validated the adaptive GMLS solution. We also note that the shape of colloids plays a significant role on their dynamics and trajectories. While the trajectory for the circular colloids was symmetric, the square colloids followed an asymmetric path and did not return to their initial positions. For convenience of comparison, the diameter of the circular colloids was set equal to the square colloids’ edge length.

In Figure 18, we present the distributions of GMLS particles at three different times corresponding to $P_{0}$, $P_{1}$, and $P_{2}$, respectively, marked on the trajectory shown in Figure 17. (Refer to the animation in the Supporting Information for more snapshots.) Since the simulation domain is large compared with the size of the colloids, if we started from a uniform, coarse distribution of particles throughout the entire domain, the computation would be still too expensive even with adaptive refinement. Thus, we instead started from an initial multi-resolution scenario. In particular, the simulation domain was divided into four sub-domains discretized with different resolutions depending on their distances to the center of the domain. The furthest subdomain was discretized with the coarsest resolution, and for each two adjacent sub-domains, the resolution ratio was set to 2. With this multi-resolution scenario as the initial distribution of GMLS particles, we applied the adaptive refinement. The error estimator directed refinement near the sharp corners of the square colloids (see Figure 18). The computed pressure field is also presented in Figure 18, from which we can see that the 4th-order GMLS solution with adaptive refinement is able to stably and correctly capture the short-range lubrication interaction between the two square colloids even when they get very close. Here, no subgrid-scale lubrication model or correction was used.

5.6 Collective dynamics of active asymmetric colloidal dimers

Lastly, we solved a more complicated but interesting problem, regarding the collective dynamics of multiple active asymmetric colloidal dimers. This study was motivated by an experimental work [3], where each dimer was composed of two overlapped spheres of different radii. As the dynamics of interest for these dimers occurred only along the substrate plane in the experiment, it may be posed as a 2D problem. Thus, each dimer was represented by two overlapped circular colloids in our simulations. While the larger lobe of each dimer had a radius of $R_{L}=4\times 10^{-6}~{}$\mathrm{m}, the radius of the smaller lobe was: $R_{S}=2.8\times 10^{-6}~{}$\mathrm{m}. The center-to-center distance of each dimer’s two lobes was equal to $R_{L}$. A collection of such colloidal dimers were suspended in a square domain with each edge length of $L=100\times 10^{-6}~{}$\mathrm{m}. The snapshots from the simulations shown in Figures [19](#S5.F19) and [20](#S5.F20) also illustrate the configuration of this problem. The fluid was assumed to be water with $\rho=1.0\times 10^{3}~{}$\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3} and $\nu=1.0\times 10^{-6}~{}$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$$. Due to the small length scale and slow dynamics [3], this may be treated as a Stokes problem.

According to the experimental work [3], each colloidal dimer can lie on the substrate surface and be actively moving in the fluid environment, propelled by an external field. The active motion of each dimer is along the center line of the two lobes toward the smaller lobe. As any two dimers come sufficiently close, a repulsive force acts between the two larger lobes and also between the two smaller lobes of the two dimers. Due to the presence of a standing, stationary dimer, the active, lying dimers nearby can be attracted toward it and form clusters of 2 to 6 dimers. It is the interplay and competition between the aforementioned active motion, repulsive force, and attraction that determine the collective dynamics and clustering of the dimers. In our simulations, the active motion of each dimer was driven by a constant force $F_{e}=1\times 10^{-9}~{}$\mathrm{N}$$ exerted along the center line of its two lobes toward the smaller lobe. The repulsive force between lobes was modeled by: $F_{r}=R^{2}(\frac{2R}{r})^{4}$, where $R=R_{L}$ or $R_{S}$, and $r$ is the center-to-center distance of the two interacting lobes. This repulsive force was applied along the center line between the two larger lobes or between the two smaller lobes of two dimers. Note that a torque could be generated for each dimer by this repulsive force. The attraction between dimers was represented by a force $F_{a}=R^{2}(\frac{R}{r})^{2}$, where $R=R_{L}$ or $R_{S}$, and $r$ is the distance from the center of the larger/smaller lobe to the standing dimer. The force was exerted at the centers of the two lobes of each dimer and toward the position of the standing dimer. A torque could also be generated for each dimer by this attractive force.

Using the adaptive GMLS, we simulated two groups of dimers: the first with three dimers and the second with six dimers. In Figure 19, three dimers initially lying apart may form a chiral cluster. The chirality (left-handedness or right-handedness) depends on the initial relative orientations of the three dimers. The equilibrated chiral cluster of three dimers rotates clockwise or counter-clockwise, depending on its chirality, with a steady rotational speed. (Refer to the animation in the Supporting Information.) In Figure 20, six dimers initially set apart finally assemble into an achiral cluster, regardless of their initial relative orientations. The equilibrated achiral clusters of six dimers remain stationary without rotating or moving. (Refer to the animation in the Supporting Information.) These different collective and clustering behaviors are determined by the competition between the active motions of dimers and the repulsive and attractive forces between the dimers. Our predictions for both groups of dimers are consistent with the experimental observations [3]. We hence demonstrate that the proposed adaptive GMLS is applicable to address more complicated but practical applications.

6 Conclusion

We have introduced a spatially adaptive high-order meshless method. It employs the GMLS discretization for solving the governing PDEs, and hence inherits the GMLS’ meshless nature and attributes of consistency, compatibility, and high-order accuracy. By applying the recovery-based a posteriori error estimator as the adaptive criterion, the method introduces spatial adaptivity in a rigorous and optimal way. Following the four-step adaptive algorithm, the spatial resolution can be adaptively refined near point singularities and fluid-solid interfaces. The method is found flexible and robust in solving various div-grad and Stokes problems with corner singularity, discontinuous diffusion coefficient, and moving boundaries of non-trivial geometries. The proposed recovery-based error displayed the same convergence order as the true errors, providing accurate estimations of the true errors. The numerical solutions produced by the adaptive GMLS method could achieve optimal convergence orders even in the presence of singularities, matching the convergence for non-singular solutions and regular domains in the literature [14, 15]. Compared with the GMLS solutions with uniform resolution, the adaptive GMLS solutions yield the same accuracy but used much fewer degrees of freedom. In the simulations of colloidal suspensions, by employing adaptive refinement in lubrication gaps and sharp representation of fluid-solid interfaces of arbitrary geometries, the short-range lubrication interaction between colloids was resolved accurately and efficiently without invoking any artificial subgrid-scale lubrication model. Thus, it is anticipated that the proposed adaptive GMLS method can be applied for solving inherently multiscale problems.

Looking ahead, the adaptive GMLS can be extended for solving FSI in a wide range of Reynolds numbers, beyond the Stokes problems. Its meshless nature will allow for advecting the GMLS particles with the flow and hence solving the Navier-Stokes equation in a Lagrangian sense. In addition, although the simulations herein were conducted in 2D, the proposed method is directly applicable to a 3D implementation, which insofar efficiency gains are concerned, stands to benefit even more from the proposed adaptive refinement strategy. With these extensions, we expect that the proposed spatially adaptive high-order meshless method can be applied to address a multitude of FSI applications.