Code verification of non‐linear immersed boundary simulations using the method of manufactured solutions

Non‐standard finite element technologies, such as immersed boundary approaches, are typically based on novel algorithms and advanced methods, which require reliable testing of the implemented code. For this purpose, the method of manufactured solutions (MoMS) offers a great framework, enabling an easy and straightforward derivation of closed‐form reference solutions. In this contribution, the focus is kept on non‐linear analysis via the finite cell method (FCM), which is typically based on an unfitted geometry discretization and higher‐order shape functions. The code verification via MoMS generally requires the application of boundary conditions to all boundaries of the simulation domain, which need to be enforced in a weak sense on the immersed boundaries. To avoid this, we propose a novel way of deriving manufactured solutions, for which the necessary constraints on the embedded boundaries are directly fulfilled. Thus, weak boundary conditions can be eliminated from the FCM simulation, and the simulation complexity is reduced when testing other relevant features of the immersed code. In particular, we focus on finite strain analysis of 3D structures with a Neo‐Hookean material model, and show that the proposed technique enables a reliable code verification approach for all load steps throughout the deformation process.


INTRODUCTION
In immersed boundary methods, such as the finite cell method (FCM) [1], the considered physical domain Ω phys generally has a complicated shape, such that a geometry-conforming spatial discretization is often cumbersome to achieve.To overcome this, immersed approaches are based on embedding the region of interest into a larger domain of simple shape Ω e , which is easy to discretize by Cartesian meshes.However, this inevitably gives rise to a fictitious domain Ω f ict = Ω e ⧵Ω phys of theoretically zero stiffness and cells intersected by domain boundaries.The unfitted discretization approach of the FCM is depicted in Figure 1 in the context of large deformations (see Section 2 for more context).The cut cells pose several challenges when it comes to accurate numerical integration, enforcement of boundary conditions, and conditioning of the system.These issues become even more significant in non-linear analyses, where stability and computation time play a crucial role [2][3][4].
F I G U R E 1 Simulation of large deformations using the finite cell method (reproduced from Ref. [2]).
For verification of the advanced methods needed to overcome numerical challenges posed by the cut cells in the FCM, we see a great benefit using the method of manufactured solutions (MoMS), whose key idea can be summarized as follows [5,6]: Instead of seeking a solution  to a partial differential equation (PDE), the problem is reversed, and one manufactures (assumes)  with known spatial and temporal derivatives.This can be easily inserted into the PDE, which then typically yields a source term (in our context a body load ).Then, if  is applied to the numerical framework, while also prescribing the manufactured solution as boundary conditions (either Dirichlet or Neumann), the numerical solutions  h should yield an approximate reproduction of the manufactured one, provided that the code is implemented correctly.Note that typically, all fields derived from  are obtained analytically1 , which makes global and local error estimation exceptionally accurate and easy without the need for overkill FEM solutions, or rarely available analytical reference solutions.

Kinematics and constitutive relations
In this contribution, both geometrical and physical non-linearities are considered in the context of finite strain analysis with a hyperelastic material based on a Neo-Hookean material model.Let  and  = () denote the position vectors of a material point in the initial and current configurations, respectively, where  is the deformation map based on the displacement field () =  + ().Without derivation, all the relevant quantities regarding the kinematics and constitutive relations are listed in Equations ( 1)-(7) [7].Here, Grad(⋅) and Div(⋅) stand for the gradient and divergence operators w.r.t., and the material parameters are represented by the Lamé constants  and .
Deformation gradient ∶  =  + Grad() Right Cauchy-Green tensor ∶  =  T Green-Lagrange strain tensor ∶  = 1 2 ( − ) Strain energy functional (Neo-Hooke mat.) ∶  =  2 (tr() − 3) +  4 ( 2 − 1) − Second Piola-Kirchoff stress tensor First Piola-Kirchoff stress tensor ∶  = In the static case, the equilibrium in the initial configuration is defined by Equation (8), where () denotes the body loads.Furthermore, as depicted in Figure 1, the boundary is composed by the Dirichlet and Neumann boundaries, on which prescribed displacements ū and tractions t are defined, according to Equations ( 9) and (10).Here,  0 denotes the outward pointing surface normal in the initial configuration, while Γ D 0 and Γ N 0 are the associated boundary pieces in the initial configuration.
In this contribution, for simulating finite strains, we follow the total Lagrangian formulation.Without derivation, the non-linear weak form of equilibrium in the initial configuration reads [3,7] which in the context of FCM, is formulated over Ω e .Thus, the standard weak form is extended by the indicator function () defined as  = 1 in Ω phys and  << 1 in Ω f ict , to suppress the energy contribution of the fictitious domain.The first term in Equation (11) considers the internal virtual work, while the second term corresponds to the virtual work of the external forces.Furthermore,  refers to the variation of displacements, and  is the variation of the Green-Lagrange strain tensor.Due to the highly non-linear nature of Equation (11), it is typically solved iteratively with the help of a Newton-Raphson scheme.After taking the directional derivative of   e , the linearized weak form reads where Δ and Δ define the increment of the displacement vector and the Green-Lagrange strain tensor, respectively.

Discretization of the weak form in the initial configuration
The linearized weak form over Ω e in Equation ( 12) can now be discretized in the initial configuration using a Cartesian mesh consisting of  c finite cells, as depicted in Figure 1.The discretized problem results in where   T refers to the global tangent stiffness matrix, furthermore,   int and   ext correspond to the global internal and external load vectors, respectively.Note that the former two depend on the current value of the displacement field.The parameter Λ  ∈ [0, 1] defines a load factor at the load increment .The equation system in Equation ( 13) is solved at every Newton-Raphson iteration  for obtaining Δ +1 .The global quantities in Equation ( 13) are obtained by the assembly of the cell-specific local quantities following the standard Bubnov-Galerkin approach.The tangential stiffness matrix  , T , the internal load vector  , int , and the external load vectors  , body and  , surf of the cell  at the Newton iteration step  are defined as Here, ℂ V and  V are the elasticity and second Piola-Kirchoff stress tensors in Voigt notation, respectively.Furthermore,  is a matrix containing the shape functions,  corresponds to the discrete gradient operator and  is the strain operator.

Problem statement
Let us now consider a rectangular plate of size 4 × 4 × 1 with a cylindrical hole Ω void of radius  = 2 in its lower left corner as depicted in Figure 2A.We formulate the manufactured solution in cylindrical coordinates, where a radial displacement field is assumed Here, the framework of the MoMS allows for an arbitrary choice for  r .For the current manufactured solution, we define  r = Λ  ( − ) 2⏟⎴⏟⎴⏟ ûr , where ûr provides the qualitative shape of the displacement field dependent on .Additionally, we also include the load increment Λ  as an additional parameter.The displacement field is depicted in Figures 2A and 2B 2 for Λ  = 1.Note that  r is constructed, such that  r ( = ) = 0 and d r d for all values of Λ  , which will be relevant later for Equation (21).Below, two important facts regarding the current manufactured solution are highlighted: 1. Following from the radial displacement field, when expressed in cylindrical coordinates, all the derived tensors are of the simple forms given below3 2. Furthermore, due to including Λ  in  r , all the derived quantities are also functions of Λ  , that is, the analytical state of the manufactured problem is known for every load step.Note, however, that due to the non-linear expressions in Equations ( 2)-( 7), the dependency is no longer linear.This is the case for the body load field as well which contains higher powers of Λ  .Due to the analytically derived expression being very lengthy for the body load, we do not include it in its full detail in this contribution, but we can provide the Mathematica script upon request.
Our aim is to manufacture the displacement field , such that tractions on the hole's boundary Ω void are vanishing.Following from Equation (18), and from the fact that  0 = −[1, 0, 0] T for Ω void , Equation (10) reduces to the simple 1D condition   () =  rr () = 0 for  =  . (20) The analytical expression for  rr () based on a displacement field of the form  r () = Λ  ûr () term reads which, based on Equation (17), indeed vanishes at  =  due to the second bracket-term evaluating to zero.In fact, for the current choice of  r ,  rr ( = ) = 0 for all values of Λ  , that is, Equation ( 20) is fulfilled throughout the entire deformation process.Thus, the appropriate boundary conditions with a combination of  ensure that manufactured solution can be reproduced by the FCM approach for all load steps without the need for weak boundary conditions.The values of  rr and  r over  are depicted in Figure 2B.Note that field values are not smooth for  < , however, we are not interested in the results within the fictitious domain.

Results
Let us now reproduce the manufactured solution of Figure 2A using the FCM framework, as depicted in Figure 3A.On the two faces highlighted in blue, the manufactured displacement field is prescribed to the numerical solution  h = .
The cylindrical void region's boundary is treated as a free boundary due to Equation (20), while on all the remaining faces,  h is restricted in the normal direction.Furthermore, the material properties are chosen as  = 50 and  = 0.34 .Note that the simulation is conducted in Cartesian coordinates, thus, when applying the prescribed displacements and body loads, an appropriate transformation of the coordinates is required.For the numerical integration of Equation ( 14) in the cut cells (red cells in Figure 3A), several approaches could be applied that were shown to yield reliable results even for non-linear analysis [2,3].In this example, however, due to the cylindrical shape of the void region, we use the smart octree numerical integration scheme based on integration sub-cells, that are perfectly aligned to the curved boundary via blending functions [8].Thus, the integration error is basically eliminated from the current problem.For the integration points within the fictitious domain,  = 10 −12 is used.We perform convergence studies for  = 1 and 2 (Serendipity family), while using a uniform discretization by 3 2 , 5 2 , 7 2 , 9 2 and 11 2 finite cells in the  1 - 2 -plane, and 2 finite cells in the  3 -direction 5 .In all cases, 20 load steps are used.Figure 3B depicts the highly deformed deformation state at the last load step using the finest discretization and  = 2. Already here, basically vanishing radial displacements at the hole's boundary can be observed 6 , indicating a highly accurate solution.In Figure 4A, the strain energy is depicted over the different load steps for  = 1 and 2. In both cases, the numerical results correspond to a discretization by 3 × 3 × 2 finite cells.Figure 4A reveals that already for  = 2, the manufactured solution is very well reproduced, not only at the end of the simulation, but throughout the entire deformation process.This is due to the fact, that owing to Equation ( 19), the hole's boundary is traction-free for all values of Λ  .For further validation of the approach, the relative error of the displacement field in the  2 -norm is evaluated for the last deformation step.On the right, the  2 -norm for an arbitrary vector field  is defined.Furthermore, || −  h ||  2 is computed on the FCM mesh, while the reference value ||||  2 is evaluated analytically using Mathematica.
The corresponding results are depicted in Figure 4B, where the theoretical algebraic convergence rates are obtained, as indicated by the black lines. 7,8 Finally, in Figure 5, the local errors are investigated by evaluating the absolute deviations | r −  h r | and | rr −  h rr | for the last solution step computed with the finest discretization and  = 2. Figure 5A reveals, that while the biggest deviations occur at the immersed boundary, even these are negligibly small, which is in good agreement with Figure 3B.Similarly, in Figure 5B, very low errors at the free boundary can be observed.

CONCLUSION
In this contribution, the MoMS was presented as a fast and reliable code verification approach for immersed boundary methods, such as the FCM.We have shown, that with some additional care, analytical manufactured solutions can be derived even for non-linear analyses, such that boundary conditions on the immersed boundaries are readily fulfilled throughout the entire deformation process.Thus, a weak imposition of Dirichlet boundary conditions can be completely avoided in the simulation.Furthermore, the derived manufactured solution is available in a closed-form for load steps, enabling a highly accurate and simple error estimation of the entire non-linear FCM simulation.Finally, we have shown that with the carefully derived body loads, the reference solution can be indeed reproduced, and theoretical convergence √  2 1 +  2 2 , the square root term cannot be fully eliminated from the displacement field.Thus ûr () cannot be exactly approximated by the polynomial Ansatz space of the FCM simulation.
rates when performing various ℎ-refinements can be obtained.This indicates not only an error-free implementation of the non-linear FCM code, but also the suitability of the proposed code verification framework itself.

2
Various fields of the manufactured solution depicted in a problem domain and along a cut line.

F I G U R E 3 4
Simulation of the manufactured solution in the context of the FCM.The red region indicates the cut cells intersected by Ω void .FCM, finite cell method.Global accuracy of the when reproducing the manufactured solution.FCM, finite cell method.

F I G U R E 5
Contour plots for evaluating the local deviation of the numerical solution from the manufactured one.
Open access funding enabled and organized by Projekt DEAL.O R C I D Márton Petö https://orcid.org/0000-0002-4373-2848R E F E R E N C E S