Fictitious Domain Method Research Articles

Improving accuracy of volume penalised fluid-solid interactions

We analyse and improve the volume-penalty method, a simple and versatile way to model objects in fluid flows. The volume-penalty method is a kind of fictitious-domain method that approximates no-slip boundary conditions with rapid linear damping inside the object. The method can then simulate complex, moving objects in general numerical solvers without specialised algorithms or boundary-conforming grids. Volume penalisation pays for this simplicity by introducing an equation-level error, the model error, that is related to the damping time η≪1. While the model error has been proven to vanish as the damping time tends to zero, previous work suggests convergence at a slow rate of O(η1/2). The stiffness of the damping implies conventional volume penalisation only achieves first order numerical accuracy. We analyse the volume-penalty method using multiple-scales matched-asymptotics with a signed-distance coordinate system valid for arbitrary smooth geometries. We show the dominant model error stems from a displacement length that is proportional to a Reynolds number Re dependent boundary layer of size O(η1/2Re−1/2). The relative size of the displacement length and damping time leads to multiple error regimes. Our key finding derives a simple smoothing prescription for the damping that eliminates the displacement length and reduces the model error to O(η) in all regimes. This translates to second order numerical accuracy. We validate our findings in several comprehensive benchmark problems and finally combine Richardson extrapolation of the model error with our correction to further improve convergence to O(η2).

A Numerical Model for Steel Continuous Casting Problem in a Time-variable Domain

A mathematical model and numerical method for simulation of the continuous casting process in a variable in time domain are presented. The variable geometry of the slab is caused by the change in time of the width of the mould. The mathematical model of the process is a Stefan problem with prescribed convection and non-linear Robin boundary condition. Considered differential equation is approximated by a finite difference scheme, which is constructed in several steps. First, a semi-discrete problem is constructed using the method of characteristics with respect to the time variable. Then, at each time level, the current elliptic problem in a curvilinear domain is replaced by a problem in the parallelepiped domain using the fictitious domain method. Finally, the boundary-value problem in the parallelepiped domain is approximated by a finite-difference scheme. The constructed fully discrete problem in algebraic form is a system of nonlinear equations containing a diagonal monotone operator and a linear part with a symmetric and positive definite $$M$$ -matrix. To solve the resulting system of nonlinear algebraic equations, well-known iterative solution methods can be applied.

A higher‐order Trace finite element method for shells

A higher-order fictitious domain method (FDM) for Reissner-Mindlin shells is proposed which uses a three-dimensional background mesh for the discretization. The midsurface of the shell is immersed into the higher-order background mesh and the geometry is implied by level-set functions. The mechanical model is based on the Tangential Differential Calculus (TDC) which extends the classical models based on curvilinear coordinates to implicit geometries. The shell model is described by PDEs on manifolds and the resulting FDM may typically be called Trace FEM. The three standard key aspects of FDMs have to be addressed in the Trace FEM as well to allow for a higher-order accurate method: (i) numerical integration in the cut background elements, (ii) stabilization of awkward cut situations and elimination of linear dependencies, and (iii) enforcement of boundary conditions using Nitsche's method. The numerical results confirm that higher-order accurate results are enabled by the proposed method provided that the solutions are sufficiently smooth.

A fictitious domain method with distributed Lagrange multipliers on adaptive quad/octrees for the direct numerical simulation of particle-laden flows

In this article, we extend our Distributed Lagrange Multiplier/Fictitious Domain method previously implemented on simple regular Cartesian grids to quadtree/octree adaptive grids. The objective is to improve both the accuracy and efficiency of our DLM/FD particle-resolved simulation method by extending its computing capabilities through dynamic local mesh refinement. The main features of our numerical method, such as a first-order operator splitting time algorithm and a second-order reconstruction of the velocity field close to the boundary of the immersed rigid bodies (of arbitrary shape), are unchanged. We implemented our adaptive DLM/FD algorithm within Basilisk, a parallel platform to solve partial differential equations on dynamic quadtree/octree grids. The quadtree/octree structure of the grid and specific design rules of Basilisk impose a special treatment of some of the operations performed on the grid in the DLM/FD-Uzawa algorithm. The new computational method is then tested and validated on a set of flow configurations including the challenging problem of accurately computing lubrication interaction forces without resorting to using any ad hoc correction. Finally, we illustrate the potential of our code to compute complex particle-laden flow configurations that were not attainable in the past with a DLM/FD algorithm implemented on a simple regular Cartesian grid.

Combining boundary‐conforming finite element meshes on moving domains using a sliding mesh approach

AbstractFor most finite element simulations, boundary‐conforming meshes have significant advantages in terms of accuracy or efficiency. This is particularly true for complex domains. However, with increased complexity of the domain, generating a boundary‐conforming mesh becomes more difficult and time consuming. One might therefore decide to resort to an approach where individual boundary‐conforming meshes are pieced together in a modular fashion to form a larger domain. This article presents a stabilized finite element formulation for fluid and temperature equations on sliding meshes. It couples the solution fields of multiple subdomains whose boundaries slide along each other on common interfaces. Thus, the method allows to use highly tuned boundary‐conforming meshes for each subdomain that are only coupled at the overlapping boundary interfaces. In contrast to standard overlapping or fictitious domain methods the coupling is broken down to few interfaces with reduced geometric dimension. The formulation consists of the following key ingredients: the coupling of the solution fields on the overlapping surfaces is imposed weakly using a stabilized version of Nitsche's method. It ensures mass and energy conservation at the common interfaces. Additionally, we allow to impose weak Dirichlet boundary conditions at the nonoverlapping parts of the interfaces. We present a detailed numerical study for the resulting stabilized formulation. It shows optimal convergence behavior of the interface coupling for both Newtonian and generalized Newtonian material models. Simulations of flow of plastic melt inside single‐screw as well as twin‐screw extruders demonstrate the applicability of the method to complex and relevant industrial applications.

Dual fat boundary method: the fat boundary method in elasticity with an extension of the application scope

The fat boundary method (FBM) is a fictitious domain method, proposed to solve Poisson problems in a domain with small perforations. It can achieve higher accuracy around holes, which makes it very suitable to solve elasticity problems because stress concentrations often appear around holes. However, there are some strict restrictions of the FBM limiting the wide range of applications. For example, the original FBM deals with perforated rectangular domain with only Dirichlet boundary conditions. Furthermore, because the global domain is extended to the holes, analytical solutions in holes corresponding to the Dirichlet boundary conditions around holes are required. This limits both the boundary conditions around holes and the shape of holes, because for arbitrary holes it is difficult to get the analytical solutions. This article makes an attempt to break these limitations and apply the FBM to elasticity. Firstly, we review the FBM and introduce Neumann boundary conditions to the rectangular domain. A mathematical proof of the conditional convergence of the algorithm is presented. Furthermore, the FBM is compared with the Lagrange multiplier method to clarify that the FBM is one kind of weak imposition method. Then we apply the FBM to linear elasticity and the dual fat boundary method is proposed to solve problems without analytical solutions in holes. Some numerical examples are presented to verify the method proposed here.

Effects of the collision model in interface-resolved simulations of particle-laden turbulent channel flows

The effects of the particle collision model in a direct-forcing fictitious domain method on the fluid and particle statistics of a fully developed turbulent channel flow laden with finite-size neutrally buoyant particles are numerically investigated. The particle collisions are described by a combination of the discrete element method and the lubrication force correction model. We first validate our code via several benchmark tests, including the normal particle–wall collisions at different impact Stokes numbers and the oblique collisions with varied incidence angles. Subsequently, the effects of the lubrication correction and the particle stiffness on the fluid and particle statistics of the particle-laden turbulent flows are examined. The results show that the lubrication force correction has an important effect on the particle pair statistics at the near-contact regime. Both the lubrication force between the particles and the decrease in the particle stiffness result in the decrease in the flow friction mainly due to the increase in the fluid Reynolds stress. The flow friction is always larger for smaller particles at the same particle volume fraction irrespective of the lubrication correction. The particle–particle lubrication force decreases the near-wall particle concentration, whereas the particle–wall lubrication force has the opposite effect.

Fractional centered difference scheme for high-dimensional integral fractional Laplacian

In this work we study the finite difference method for the fractional diffusion equation with high-dimensional hyper-singular integral fractional Laplacian. We first propose a simple and easy-to-implement discrete approximation, i.e., fractional centered difference scheme with γth-order (γ≤2) convergence for the fractional operator. Based on the established approximation, we then construct a finite difference scheme to solve fractional diffusion equations and analyze the stability and convergence in discrete energy norm (0<α≤2) and in discrete maximum norm (1<α≤2). We further present a fast solver for the linear system which is obtained by discretization on rectangular domain and use the fictitious domain method to extend the fast solver to the non-rectangular one. Several numerical results are provided to support our theoretical results.

Enhanced numerical integration scheme based on image compression techniques: Application to rational polygonal interpolants

Polygonal finite elements offer an increased freedom in terms of mesh generation at the price of more complex, often rational, shape functions. Thus, the numerical integration of rational interpolants over polygonal domains is one of the challenges that needs to be solved. If, additionally, strong discontinuities are present in the integrand, e.g., when employing fictitious domain methods, special integration procedures must be developed. Therefore, we propose to extend the conventional quadtree-decomposition-based integration approach by image compression techniques. In this context, our focus is on unfitted polygonal elements using Wachspress shape functions. In order to assess the performance of the novel integration scheme, we investigate the integration error and the compression rate being related to the reduction in integration points. To this end, the area and the stiffness matrix of a single element are computed using different formulations of the shape functions, i.e., global and local, and partitioning schemes. Finally, the performance of the proposed integration scheme is evaluated by investigating two problems of linear elasticity.

Simultaneous material, shape and topology optimization

Using three design fields we develop an optimization environment that can simultaneously optimize material, shape and topology. We use the implicit representation of the boundaries with level-set functions that define the shape and topology. Differentiable R-functions allow us to combine these shapes and topology descriptions with Boolean operations. Additionally, we incorporate design dependent-stiffness materials with another design field. Notably, this framework accommodates design dependent loads, has the ability to introduce holes, and ensures the satisfaction of optimality criteria. It builds upon the fictitious domain, ersatz material, material interpolation and level-set methods. It also borrows from parameterized density-based topology optimization methods. Since analytical sensitivities can be computed, we use efficient nonlinear programming algorithms to update the design instead of the Hamilton–Jacobi’s scheme of level-set methods. We illustrate the features of our framework by designing a cantilever beam with octet truss microlattice, a dam with design-dependent loads, and a composite clevis plate.

Numerical investigation on effects of alkali metals release during combustion of char particle groups

Alkali metals, especially Na, are generally accepted as the primary cause of serious ash-related problems in the industrial applications of pulverized coal. This paper investigates the effects of different variables, including particle size, ambient temperature, O2 mass fraction, and atmospheres: O2/N2 and O2/CO2 at particle level. The fictitious domain method (FDM) provides a true direct numerical simulation (TDNS) for the investigation. Computational results show the correlation between alkali metals release and carbon consumption during char particle combustion. Here, high ambient temperature is more positive to alkali metals release than carbon consumption. In an oxygen-dominated atmosphere, a small amount of CO2 tends to form a higher temperature field than N2, which promotes a greater release of alkali metals. This paper also investigates the specific effects of volume fraction and non-uniform distribution of particles for the combustion of particle groups. The computational results show that both temperature field and alkali metals release are significantly affected by the volume fraction of particles. The effects of the non-uniform distribution of particles, e.g., uniform distribution, dense distribution on the windward side, dense distribution on the leeward side, and asymmetric distribution, on temperature field and alkali metals release vary with different volume fractions of particles. As such, the non-uniform distribution of particles plays a role in changing the mass transfer of the oxidant and the heat transfer.

Consistent Velocity–Pressure Coupling for Second-Order L2-Penalty and Direct-Forcing Methods

The present work studies the interactions between fictitious-domain methods on structured grids and velocity–pressure coupling for the resolution of the Navier–Stokes equations. The pressure-correction approaches are widely used in this context but the corrector step is generally not modified consistently to take into account the fictitious domain. A consistent modification of the pressure-projection for a high-order penalty (or penalization) method close to the Ikeno–Kajishima modification for the Immersed Boundary Method is presented here. Compared to the first-order correction required for the L 2 -penalty methods, the small values of the penalty parameters do not lead to numerical instabilities in solving the Poisson equation. A comparison of the corrected rotational pressure-correction method with the augmented Lagrangian approach which does not require a correction is carried out.

Consistent Velocity-Pressure Coupling for Second-Order L 2 -Penalty and Direct-Forcing Methods

The present work studies the interactions between fictitious-domain methods on structured grids and velocity-pressure coupling for the resolution of the Navier-Stokes equations. The pressure-correction approaches are widely used in this context but the corrector step is generally not modified consistently to take into account the fictitious domain. A consistent modification of the pressure-projection for a high-order penalty (or penalization) method close to the Ikeno-Kajishima modification for the Immersed Boundary Method is presented here. Compared to the first-order correction required for the L 2-penalty methods, the small values of the penalty parameters do not lead to numerical instabilities in solving the Poisson equation. A comparison of the corrected rotational pressure-correction method with the augmented Lagrangian approach which does not require a correction is carried out.

Enhanced numerical integration scheme based on image-compression techniques: application to fictitious domain methods

In the present work, we propose a new approach, the so-called compressed adaptive integration scheme (C-AIS), for the computation of the stiffness and mass matrices in fictitious domain methods requiring the integration of discontinuous functions. The novel approach extends the conventional quadtree-decomposition-based adaptive integration scheme (AIS) by an additional step, in which established image-compression techniques are exploited to decrease the number of integration sub-cells. The benefits of the C-AIS are manifold: First, the compression of the sub-cells inevitably leads to significant savings in terms of computational time required by the numerical integration. Second, the compression procedure, which is executed directly after the quadtree-decomposition algorithm, can be easily included in existing codes. Third, if applied to polynomial integrands, the C-AIS yields exactly the same accuracy as the conventional AIS. Finally, the fourth advantage is seen in the fact that the C-AIS can readily be combined with other approaches seeking a reduction of the number of integration points such as the Boolean-FCM. The efficiency of the C-AIS approach is presented in the context of the FCM based on Cartesian meshes applied to problems of linear elastostatics and modal analysis, while it is also suitable for the quadrature in other fictitious domain approaches, e.g., CutFEM and cgFEM.

A unified finite strain theory for membranes and ropes

The finite strain theory is reformulated in the frame of the Tangential Differential Calculus (TDC) resulting in a unification in a threefold sense. Firstly, ropes, membranes and three-dimensional continua are treated with one set of governing equations. Secondly, the reformulated boundary value problem applies to parametrized and implicit geometries. Therefore, the formulation is more general than classical ones as it does not rely on parametrizations implying curvilinear coordinate systems and the concept of co- and contravariant base vectors. This leads to the third unification: TDC-based models are applicable to two fundamentally different numerical approaches. On the one hand, one may use the classical Surface FEM where the geometry is discretized by curved one-dimensional elements for ropes and two-dimensional surface elements for membranes. On the other hand, it also applies to recent Trace FEM approaches where the geometry is immersed in a higher-dimensional background mesh. Then, the shape functions of the background mesh are evaluated on the trace of the immersed geometry and used for the approximation. As such, the Trace FEM is a fictitious domain method for partial differential equations on manifolds. The numerical results show that the proposed finite strain theory yields higher-order convergence rates independent of the numerical methodology, the dimension of the manifold, and the geometric representation type.

Direct numerical simulation of particle pore-scale transport through three-dimensional porous media with arbitrarily polyhedral mesh

A new direct numerical simulation algorithm is developed for particle pore-scale transport through the porous media with arbitrarily polyhedral mesh. In the algorithm, the Navier-Stokes Equation is used to describe the continuous phase motion in the Eulerian framework; Newton's Second Law is used to describe the particle dynamics in the Lagrangian framework; Discrete element method is used to describe the particle-particle interactions and particle-wall interactions; RIGID is used to detect the contact state between particles with arbitrarily shaped pore walls. To suppress the spurious force oscillations (SFO) and improve the numerical accuracy of the evaluation of fluid-particle interaction, a novel consistent fictitious domain method (CFDM) in the arbitrarily collocated polyhedral mesh is developed. Numerical results of six test cases show that CFDM is accurate and second order in space, and no obvious SFO is found. Finally, the new direct numerical simulation algorithm is used to simulate the particle transport through three-dimensional porous media reconstructed from micro-CT scans from a real rock. The numerical results of a serial of tests with different particle sizes reveal several distinct microscopic flow mechanisms and the corresponding macroscopic characteristics. The change of channel resistance leads to the formation of particle motion paths in succession; Along a certain motion path, the particle moving velocity can be different at different sites; With the increase of particle size, the particle average retention time and particle average transit time increase; Particle velocity presents lognormal distribution, which becomes wider with the increase of particle size. The newly developed algorithm can be adopted as a direct numerical simulation tool to simulate particle motion in arbitrarily complex pore space.

Interface-resolved direct numerical simulations of the interactions between spheroidal particles and upward vertical turbulent channel flows

The interactions between finite-size spheroidal particles and upward turbulent flows in a vertical channel are numerically simulated with a direct-forcing fictitious domain method at two particle settling coefficients $u_{s}$ (the ratio of the particle Stokes free-fall velocity to the bulk velocity) of 0.1 and 0.3, a bulk Reynolds number of 2873, a ratio of the particle equivalent diameter to the channel width of 0.05, a particle volume fraction of 2.36 % and particle aspect ratios of $1/3$ , 1 and 2. Our results show that the flow friction is largest for the case of a sphere, and smallest for the oblate case when the particle sedimentation effect is weak ( $u_{s}=0.1$ ), whereas the flow friction is smallest for the case of a sphere, and largest for the oblate case when the particle sedimentation effect is moderately strong ( $u_{s}=0.3$ ). The reason for the lower flow friction of the spherical particles is that the large-scale vortices are more strongly attenuated by the spherical particles than by the non-spherical particles in the case of $u_{s}=0.3$ . The settling particles tend to migrate towards the channel centre due to the Saffman effect, and the migration is strongest for the spherical particles. The non-spherical particles tend to align their long axes with the streamwise direction in the near-wall region, and perpendicular to the streamwise direction in the bulk region due to the significant settling effect.

Modelling of hydro-mechanical processes in heterogeneous fracture intersections using a fictitious domain method with variational transfer operators

Fluid flow in rough fractures and the coupling with the mechanical behaviour of the fractures pose great difficulties for numerical modeling approaches, due to complex fracture surface topographies, the non-linearity of hydromechanical processes and their tightly coupled nature. To this end, we have adapted a fictitious domain method to enable the simulation of hydromechanical processes in fracture-intersections. The main characteristic of the method is the immersion of the fracture, modelled as a linear elastic solid, in the surrounding computational fluid domain, modelled with the incompressible Navier Stokes equations. The fluid and the solid problems are coupled with variational transfer operators. Variational transfer operators are also used to solve contact within the fracture using a dual mortar approach and to generate problem specific fluid meshes. With respect to our applications, the key features of the method are the usage of different finite element discretizations for the solid and the fluid problem and the automatically generated representation of the fluid-solid boundary. We demonstrate that the presented methodology resolves small-scale roughness on the fracture surface, while capturing fluid flow field changes during mechanical loading. Starting with 2D/3D benchmark simulations of intersected fractures, we end with an intersected fracture composed of complex fracture surface topographies, which are in contact under increasing loads. The contributions of this article are: (1) the application of the fictitious domain method to study flow in fractures with intersections, (2) a mortar based contact solver for the solid problem, (3) generation of problem specific grids using the geometry information from the variational transfer operators.

Interface-resolved numerical simulations of particle-laden turbulent channel flows with spanwise rotation

Interface-resolved simulations of particle-laden turbulent channel flows with spanwise rotation at a Reynolds number of 180 and different rotation numbers ranging from 0.1 to 1.0 are performed with a fictitious domain method. The difficulty of the centrifugal force on the particles not satisfying the periodic boundary condition is circumvented by the feature of the fictitious domain formulation for the neutrally buoyant case, where the centrifugal force in the particle motion equation vanishes, and by only considering a low rotation number of 0.1 and setting the rotation center to be far away from the channel for the non-unity density ratio case. Our results show that the heavy particles (i.e., the particle density being larger than the fluid density) migrate towards the pressure wall, whereas the light particles migrate towards the suction wall. For the density ratio being unity, the particle concentration is higher near the pressure wall than near the suction wall, and we attribute the reason to the effects of the mean secondary flow structure (i.e., the Taylor–Görtler vortices), since similar particle concentration distribution and secondary flow structure are observed in a rotating laminar channel flow. The mean velocities of heavy particles are smaller in the pressure-side half channel except the near-wall region, and larger in the suction-side half channel, compared to the fluid mean velocity; the opposite occurs for the light particle case. The addition of the finite-size particles increases the flow drag. The flow drag is not sensitive to the density ratio for the light particles and increases with increasing density ratio for the heavy particles. The effects of the particles on the fluid root-mean-square velocities of the rotating turbulent channel flow are generally similar to the non-rotating channel case, but become more complicated because of the asymmetric turbulence intensity and particle concentration distribution near two walls caused by the channel rotation.

The Algebraic Immersed Interface and Boundary Method for Elliptic Equations with Jump Conditions

A new fictitious domain method, the algebraic immersed interface and boundary (AIIB) method, is presented for elliptic equations with immersed interface conditions. This method allows jump conditions on immersed interfaces to be discretized accurately. The main idea is to create auxiliary unknowns at existing grid locations which increases the degrees of freedom of the initial problem. These auxiliary unknowns allow to impose various constraints to the system on interfaces of complex shapes. For instance, the method is able to deal with immersed interfaces for elliptic equations with jump conditions on the solution or discontinuous coefficients with a second order of spatial accuracy. As the AIIB method acts on an algebraic level and only changes the problem matrix, no particular attention to the initial discretization is required. The method can be easily implemented in any structured grid code and can deal with immersed boundary problems too. Several validation problems are presented to demonstrate the interest and accuracy of the method.