An eigenvalue stabilization technique to increase the robustness of the finite cell method for finite strain problems

Broken cells in the finite cell method—especially those with a small volume fraction—lead to a high condition number of the global system of equations. To overcome this problem, in this paper, we apply and adapt an eigenvalue stabilization technique to improve the ill-conditioned matrices of the finite cells and to enhance the robustness for large deformation analysis. In this approach, the modes causing high condition numbers are identified for each cell, based on the eigenvalues of the cell stiffness matrix. Then, those modes are supported directly by adding extra stiffness to the cell stiffness matrix in order to improve the condition number. Furthermore, the same extra stiffness is considered on the right-hand side of the system—which leads to a stabilization scheme that does not modify the solution. The performance of the eigenvalue stabilization technique is demonstrated using different numerical examples.


Introduction
In fictitious domain or immersed boundary methods, such as the finite cell method (FCM) [13,32], CutFEM [6][7][8], or CutIGA [15], the cells/elements generally do not conform to the boundary of the geometry of interest. As a result, the geometry can be easily discretized using, for instance, Cartesian grids. Nevertheless, this simplification in the mesh generation leads to broken cells/elements that are cut by the boundary of the domain. This can result in an ill-conditioned global system matrix, especially if badly broken cells/elements are filled only with a small fraction of material. In the extended finite element method (XFEM), conditioning problems can also arise. An overview of different approaches to overcome the conditioning problem is given in [29]. In nonlinear finite strain analysis, badly broken B Wadhah Garhuom wadhah.garhuom@tuhh.de Khuldoon Usman khuldoon.usman@tuhh.de Alexander Düster alexander.duester@tuhh.de 1 Numerical Structural Analysis with Application in Ship Technology (M-10), Hamburg University of Technology, Am Schwarzenberg-Campus 4 (C), Hamburg 21073, Germany cells/elements tend to deform significantly compared to nonbroken cells which leads to convergence problems during the iterative/incremental procedure.
To this end, different stabilization techniques have been developed for fictitious domain methods to overcome the problem of ill-conditioning. In the context of the CutFEM, Burman et al. [6][7][8] proposed a stabilization approach called the ghost penalty method, to obtain a better condition number for elements that are cut by the boundary. To this end, an extra term is added to the weak form so as to penalize jumps of the normal derivatives of the shape functions of the neighboring cut and non-cut elements. Consequently, the cut elements are supported by the interior non-cut elements. Here, the penalty term is referred to as ghost since it is related to the fictitious domain.
In the context of the FCM, the fictitious material or αmethod is suggested in [13,32], where a very soft material is introduced in the fictitious domain of every broken cell. This is done by scaling the material parameters with the indicator function α, which has a value of one for points inside the physical domain and a value of α = 10 −q for points outside the physical domain. For practical problems, we set q = 5, . . . , 9 to achieve reasonable results. To improve the condition number without changing the problem under consideration substantially, it is important not to add too much stiffness in the fictitious domain [10]. This approach also helps to control the solution within the fictitious domain while at the same time improving the condition number. This is especially necessary when simulating nonlinear finite strain problems to prevent convergence issues. Generally, the approach based on the fictitious material works well for linear and nonlinear problems under small strains [10,13,22,40]. One drawback of this approach is that the same artificial stiffness is added to all points in the fictitious domain, which can result in a significant modification of the solution.
Additionally, for the analysis using high-order shape functions, degrees of freedom that have a small support in the physical domain tend to increase the condition number and decrease the robustness of the FCM. To this end, [15,17,42] proposed a basis function removal (BFR) strategy to improve the condition number by removing shape functions that do not contribute much to the global stiffness matrix. Thereby, a global criterion was developed based on the energy contribution of each degree of freedom -allowing to determine which shape function needs to be removed. In doing so, the high-order modes of the broken cells that cause the problems are removed and consequently the condition number is improved significantly.
Moreover, the severe distortion of the mesh in nonlinear finite strain problems can cause the analysis to terminate before reaching the desired deformation state. To this end, a remeshing strategy proposed in [17,18] can improve the robustness of the FCM by creating a new mesh whenever the old mesh cannot take any more deformations. Starting with an initial mesh, the structure is deformed until the remeshing criteria indicate that a remeshing step is required because the mesh is strongly distorted. Then, a new mesh is created that covers the deformed geometry. Thanks to the fictitious domain approach, creating the mesh is straightforward. Subsequently, important quantities between the old and new mesh are interpolated. Finally, the analysis is continued until the desired load step is reached. It was shown in [17] that a combination of the remeshing strategy, the basis function removal, and the α-method can improve the robustness even for very large deformations.
In [34], a preconditioning technique was developed to prevent poor conditioning of the FCM caused by broken cells. This approach is referred to as the Symmetric Incomplete Permuted Inverse Cholesky preconditioner. To this end, the preconditioning involves a diagonal scaling of the shape functions as well as an orthonormalization process utilizing the Gram-Schmidt procedure.
In this paper, we present an eigenvalue stabilization technique, based on the approach of Loehnert [29], to improve the robustness of the FCM for nonlinear problems by reducing the condition number of broken cells without affecting the solution significantly. To this end, an eigenvalue decomposition of the cell stiffness matrix is computed for every broken cell. Afterwards, the modes that have small eigenvalues or even zero eigenvalues can be grouped together. Those modes cause the FCM to be less robust in nonlinear computations. They appear especially when badly broken cells are present together with high-order shape functions. The modes with zero eigenvalues that belong to rigid body motions are not considered for stabilization-in order to avoid non-physical behavior. Finally, the bad modes are supported by adding extra stiffness to the cell stiffness matrix of the corresponding mode and to the load vector, to ensure that the solution is not modified significantly. Furthermore, we introduce a criterion that allows to adaptively add stiffness to the modes depending on their corresponding eigenvalue. This paper is organized as follows: Sect. 2 gives a brief summary of the basic formulation of the FCM. In Sect. 3, the stabilization of the FCM based on the α-method is explained. In Sect. 4, an eigenvalue stabilization technique for the FCM is discussed in detail. In Sect. 5, we study the accuracy of the proposed method on a linear benchmark example and on more complex nonlinear finite strain problems. Finally, the paper is summarized in Sect. 6.

The finite cell method
This section offers a brief overview of the basics of the finite cell method (FCM)-for more details see [13,14,32,36]. The FCM is a combination of high-order finite elements with a fictitious domain approach. To this end, hierarchic shape functions based on integrated Legendre polynomials are utilized [38,39]. The FCM was successfully applied to a variety of problems in solid mechanics such as thermoelasticity [44], geometrical nonlinearities [18,37], explicit and implicit elastodynamics [16,24], biomechanics [35,42], elastoplasticity [3,26,40], and microstructured materials [20,21,27]. Figure 1 illustrates the basic concept of the method using a twodimensional geometry of solid mechanics. Here, the physical domain is embedded into a fictitious domain e \ which results in an extended domain e of a simple shape. The extended domain e can then be easily discretized using, for example, Cartesian grids. Dirichlet boundary conditions d as well as Neumann boundary conditionst are applied on the boundary of the physical domain . For the sake of simplicity, we neglect body forces in this context. As a result of the fictitious domain approach, the approximation of the displacements is decoupled from the approximation of the geometry-in contrast to the standard finite element method (FEM). Therefore, the term cells is used to differentiate them from the geometry conforming finite elements. Next, to account for the boundary of the geometry, the indi- is introduced, which has a value of one inside the physical domain and a value of zero elsewhere. The nonlinear weak form G α e reads [43] G where P denotes the first Piola-Kirchhoff stress tensor and η corresponds to the test function. Since our focus later on will be on large deformation analysis, Eq. (2) is highly nonlinear. Therefore, we utilize the Newton-Raphson method, which is based on the linearization of the weak form. Assuming that the external loads are independent of the deformation, the directional derivative of G α e reads where and d represents the displacement increment. Accordingly, the linearized weak form results in Next, the extended domain is discretized using finite cells based on hierarchic shape functions in order to compute an approximate solution of the problem at hand, as illustrated in Fig. 1. Inserting the discretization of the displacement and the test function into the weak form results in the following linear equation system which needs to be solved within each Newton-Raphson iteration for the unknown displacement increment. Here, K i T (d i ) refers to the global tangent stiffness matrix. Furthermore, G i (d i ) defines the global out of balance vector. During the assembly process of the local stiffness matrices and vectors of each finite cell, the global quantities are acquired as follows The simplification in the mesh generation using the FCM leads to broken cells which are intersected by the boundary of the domain. This results in discontinuous integrands, see Eq. (5). Consequently, standard Gauss quadratures do not perform well for those cells. Therefore, special integration schemes need to be used-such as, for instance, the adaptive integration scheme based on spacetree decomposition [2,13,28,33], the moment fitting quadrature [11,12,19,25], or the approach based on equivalent polynomials [1,41].
Furthermore, fictitious domain methods such as the FCM often suffer from ill-conditioning of the resulting system matrices caused by the broken cells. Figure 2 shows the condition number for a plate with a hole, discretized using four finite cells applying different values of the radius r for the circle. The condition number is plotted against the polynomial degree p of the hierarchic shape functions applied in the finite cells. It can be seen that the condition number increases drastically

Material stabilization (˛-method)
The commonly used method to improve the condition number of the FCM is the so-called material stabilization technique or the α-method [10,13]. The main idea is to introduce a soft material into the fictitious domain. To this end, the indicator function α is set to a small positive value for integration points that are located in the fictitious domain In doing so, the condition number can be improved by introducing an artificial stiffness, as will be shown in Sect. 5. By setting the indicator function α to zero in the fictitious domain, the geometry is exactly preserved. To achieve reasonable results without losing too much accuracy, the q values should be chosen in a range of q = 5, . . . , 9. One advantage of this approach is the simple implementation, since a new set of Gaussian points can be introduced in the fictitious domain to reduce the ill-conditioning of the system. Figure 3 shows the number of integration points needed utilizing a shape function order of p = 4. Here, the red dots refer to the points required to integrate the physical domain and the blue dots represent the stabilization points. In Fig. 3a, the moment fitting scheme is applied with (2 p + 1) points in each direction to integrate the physical domain. Additionally, a set of ( p + 1) points in each direction is added to stabilize the fictitious domain, excluding the points that are located in The integration points used for stabilization (blue dots) and the points used for integrating the physical domain (red dots) applied on one broken cell with an ansatz order of p = 4. a Moment fitting with (2 p + 1) 2 physical points plus ( p + 1) 2 stabilization points. b Quadtree with tree depth k = 2 and ( p + 1) 2 physical points in each sub-cell plus ( p + 1) 2 stabilization points. c Quadtree with tree depth k = 2 and ( p + 1) 2 physical points in each sub-cell plus ( p − 1) 2 stabilization points. (Color figure online) the physical domain. This scheme will be used in Sect. 5.1 for the linear computations. In the nonlinear analysis (Sect. 5.2 and 5.3), the adaptive octree is applied to integrate the physical domain with ( p + 1) points in each direction of every sub-cell, excluding the points located in the fictitious domain. Additionally, a number of ( p +1) points in each direction are applied to stabilize the fictitious domain, excluding the points located in the physical domain, as can be seen in Fig. 3b for a tree depth of k = 2.
One drawback of the material stabilization approach is the loss of accuracy especially, if large α values are to be used. In the nonlinear analysis, the α-method will be combined with an eigenvalue stabilization technique (see Sect. 4) to increase the robustness of the FCM when undergoing large deformations. To this end, a set of ( p − 1) points will be added in each direction of the fictitious domain, excluding the points located in the physical domain and setting a small value of α = 1e−7 for those points, as can be seen in Fig. 3c.

Eigenvalue stabilization technique
In this section, following Loehnert [29], an eigenvalue stabilization technique will be applied to improve the conditioning and robustness of the FCM. This approach is based on the eigenvalues of the cell's stiffness matrix. Usually, cells that have a small volume fraction and are intersected by the geometry will exhibit eigenvalues that are very close to zero-in addition to the zero eigenvalues that correspond to rigid body modes (RBM). This means that some of the modes will tend to have linear dependencies which result in an increasing condition number. To overcome this problem, following the ideas of Loehnert [29,30] and Loehnert and Beese [31], we can compute the eigenvalues and eigenmodes of each broken cell by factorizing the cell stiffness matrix k c using the eigenvalue decomposition Here, k c is a symmetric positive semi-definite matrix which can be factorized into the diagonal matrix that contains the eigenvalues with the corresponding eigenmodes V . The eigenvalues consist of non-zero eigenvalues¯ , with the corresponding eigenvectorsV , and of eigenvalues that are equal to zero with the related eigenvectors V 0 . Next, the eigenmodes that correspond to a zero eigenvalue or to a small eigenvalue are grouped into the eigenspaceṼ with their corresponding eigenvalues˜ . The modes that are related to rigid body translations and rotations are removed from the eigenspaceṼ . To this end, all modes remaining within the eigenspaceṼ need to be supported. To do so, we construct the stiffness matrix Here, j defines a fixed stabilization factor for all eigenvalues, set to a small positive value, while n refers to the total number of modes that need to be supported. Theṽ j corresponds to the jth eigenmode that belongs toṼ . Alternatively, one can define a nonlinear stabilization factor γ j that depends on j as well as the corresponding eigenvalue˜ j . To this end, the additional stiffness γ j can be adaptively defined depending on the size of the eigenvalue, see with The condition of the system is improved by adding k to the cell stiffness matrix k c , which results in a modified stiffness matrix as well as to the out of balance vector g c (right-hand side-RHS) of the nonlinear computation-for the sake of consistency-as follows Finally, the Newton-Raphson method is applied to solve the modified system whereK i T (d i ) denotes the modified global tangent stiffness matrix andG i denotes the modified global residual vector, which are the result of the assembly process (7). Having solved the linear system (16), the solution is updated: The procedure for the eigenvalue stabilization is summarized in Algorithm 1. In Sect. 5, we will investigate the influence of the proposed method on the robustness and the accuracy of the FCM. Identify the modes to be supportedṼ 8: Remove the modes related to RBMs 9: Build the stabilization matrix k (12) 10: end if

Detection of the rigid body modes
The eigenvectors with zero eigenvalues that are related to the rigid body modes need to be extracted from the set of modes to be stabilizedṼ . Otherwise, stabilizing those modes would lead to a non-physical behavior of the system. To this end, one way to detect the rigid body modes of a cell is to first compute the RBMs analytically. Afterwards, the RBMs can be easily identified and removed from the spaceṼ by means of a Gram-Schmidt orthogonalization procedure, as proposed by Loehnert [29]. The RBMs can be computed cheaply based on the nodal coordinates of the cells. Generally, for instance, a hexahedral element in a three-dimensional space has six RBMs: three translation modes-which can be calculated by applying a small displacement at every node in the corresponding direction of the mode-and three rotational modes that can be computed using a spherical coordinate system, as suggested by [4,5,23]. The high-order modes of the hierarchic basis can be ignored, since the rigid body motions are represented with the nodal modes.

Numerical investigations
Before applying the proposed eigenvalue stabilization technique for nonlinear applications, we first investigate it for a linear problem. Thus, in the first example (Sect. 5.1), a linear elastic isotropic material model is applied with a Young's modulus of E = 50 MPa and a Poisson's ratio of ν = 0.3 assuming small strains. In the second and third application, (Sects. 5.2 and 5.3), a hyperelastic material model is utilized in a finite strain analysis. The strain energy density function [9] is defined as follows where C = F T F refers to the right Cauchy-Green tensor and J = det (F). The material parameters read where λ and μ denote the Lamé first parameter and the shear modulus, respectively. The first Piola-Kirchhoff stress tensor is computed as (21) and the elasticity tensor reads where I refers to the fourth order identity tensor (I i jkl = δ ik δ jl ), and M is computed in index notation as M i jkl = −F jk F li . For all numerical examples, Eq. (12) is used to compute the additional stiffness matrix of the eigenvalue stabilization.

Plate with a cylindrical hole
In the first example, we consider a plate with a cylindrical hole. The goal is to compare the performance of the eigenvalue stabilization technique with the α-method in terms of accuracy and improvement of the condition number. The results will be compared to a reference solution generated from the standard p-FEM. The geometry of the plate has a side length of a = 100 mm and a radius of r = 60 mm. Symmetry boundary conditions are applied together with a prescribed pressure on the top surface oft y = 0.1 MPa, as can be seen in Fig. 5. The reference solution is obtained with 3200 curved p-version hexahedral elements utilizing a polynomial degree of p = 5. The FCM computations are The moment fitting quadrature is used to integrate the broken cells of the FCM with (2 p + 1) 3 integration points and an octree depth of k = 8 is applied to compute the moments. In order to stabilize the cells (partially) located in fictitious domain using the α-method, the standard Gauss-quadrature is utilized with ( p + 1) 3 integration points, and the points located in the physical domain are removed, as illustrated in Fig. 3a. Applying the eigenvalue stabilization, no points are added in the fictitious domain.

Eigenvalue stabilization without correcting the RHS
First, we investigate the eigenvalue stabilization technique without modifying the load vector as is the case in Eq. (15). This is because, in a linear analysis, it is not possible to compute the correction term (k d c ) directly without knowing the solution vector d c . We study the condition number of the global stiffness matrix, plotted against the number of degrees of freedom, see Fig. 7. Here, different factors are used for (top of Fig. 7) and for α (bottom of Fig. 7). It can be seen The results for the eigenvalue stabilization are depicted at the top, the results for the α-method at the bottom. Using the α-method, we do not observe any change in the accuracy for α = 0 to α = 1e−7. However, by increasing the values from α = 1e−5 to α = 1e−3, we start to loose accuracy quite significantly as compared to the results of the eigenvalue stabilization. There, we still have a good accuracy until = 1e−2 before we start losing accuracy using = 9e−2.
For better visualization, in Fig. 9, the relative error of the von Mises stress is plotted versus the condition number for an increase of the polynomial degree of the shape functions. It can be seen that the condition number is reduced significantly when using the eigenvalue stabilization (top), without affecting the solution to a great extent as compared to the α-method (bottom). This is because, in the eigenvalue stabilization, only modes that have a small support (small eigenvalue) are stabilized-whereas a larger loss of accuracy is observed for the α-method, where all points in the fictitious domain are stabilized with the same amount.
Furthermore, we also investigate the accuracy of the strain energy. In Fig. 10, the relative error in energy norm is plotted versus the number of degrees of freedom for an increase of the polynomial degree of the shape functions. To this end, using the α-method (bottom of Fig. 10), we do not see any change in the accuracy for α = 0 to α = 1e−7. However, increasing the values from α = 1e−5 to α = 1e−3, we start to loose accuracy quite significantly as compared to the results of the eigenvalue stabilization at the top of Fig. 10.
For a better visualization, in Fig. 11, the relative error in energy norm is plotted versus the condition number. It can be seen that the condition number is reduced significantly when

Eigenvalue stabilization with correcting the RHS
In this section, we apply the eigenvalue stabilization by also modifying the load vector on the RHS. This is done in two steps. First, we compute the solution d 1 by modifying only the stiffness matrix K = K + K (23) and solving the system where F ext defines the global external load vector. In the next step, we solve the system again by correcting the load vector iteratively as follows In the first iteration (i = 1), we apply the solution d 1 , which is obtained from Eq. (24). Please note that an exact correction of the RHS would call for a repeated solution of the linear system (25) until the displacement does not change anymore. However, we found that three correction steps are already able to provide a good accuracy. By using direct solvers, it is furthermore possible to store the factorization of the global stiffness matrixK and re-use it to solve the system for multiple RHSs. In Fig. 12, the relative error of the von Mises stress (top) as well as the relative error in energy norm (bottom) are plotted against the condition number, applying three correction steps. It can be seen that the condition number is improved significantly without affecting the solution at all, provided that the load vector is corrected accordingly. In the nonlinear computations (Sects. 5.2 and 5.3), we correct the RHS during the iterative/incremental procedure using the last computed solution, as illustrated in Algorithm 1.

Cost of the stabilization
In this section, we investigate the cost of the stabilization utilizing the α-method as well as the eigenvalue approach. We define the cost of the stabilization by measuring the difference between the computation time of the whole simulation with stabilization (T s ) and without stabilization (T ) as follows The same problem setup as described in (Sect. 5.1) is used here with α = 1e−5 and = 9e−4. Here, a computer with 2 CPUs each of with 10 cores (40 threads) and 2.4 GHz is used for the simulations. In Fig. 13, the stabilization cost is plotted versus the number of degrees of freedom. It can be seen that the cost of the eigenvalue stabilization is comparable to the α-method. This shows that the cost of computing

Single cube connector
In the following application, we consider a single cube connector [17]. The motivation is to investigate the robustness of the eigenvalue stabilization technique for nonlinear computations undergoing large deformations-and to compare it to the α-method. The geometry of the cube connector is defined using the following level set function with an inner radius of r = 11.25 mm, an outer radius of R = 15 mm, and a design parameter d = 4.6×10 4 mm 4 . The Here, every point x within the reference domain (x ∈ ) corresponds to a value of the level set function that is greater or equal to zero (φ(x) ≥ 0). The strain energy function of the hyperelastic material model is given in Eq. (18) and the corresponding material parameters are presented in Eqs. (19) and (20). Symmetry boundary conditions are applied together with a prescribed displacement at the top surface (compression). The geometry is discretized with 129 cells (6 × 6 × 6 subdivisions), as illustrated in Fig. 14. For integrating the broken cells, the adaptive octree with a tree depth of k = 3 is utilized and ( p + 1) 3 integration points are distributed in each sub-cell. For stabilizing the fictitious domain using the α-method, the standard Gauss-quadrature is utilized with ( p + 1) 3 integration points, as can be seen in Fig. 3b. Applying the eigenvalue stabilization, a set of ( p − 1) 3 integration points is added in the fictitious domain with α = 1e−7, as can be seen in Fig. 3c. We apply displacement increments of 0.05 mm in each load step, trying to compress the single cube connector as much as possible.
To this end, we start by setting the ansatz order to p = 2 and apply different stabilization factors for the α-method, It can be observed that using the α-method, the structure can be deformed up to a displacement ofū z = 4.0 mm utilizing q = 0, as can be seen in Fig. 15 (bottom). However, it is clear that the structure becomes too stiff because of the large stiffness added to the system. For q values between 3 and 5, the α-method produces reasonable results with a maximum deformation ofū z = 2.85 mm. On the other hand, utilizing the eigenvalue stabilization, we can reach the final deformation state with a displacement ofū z = 7.00 mm using q = 4, which is about 2.5 times more as compared to the α-method. Additionally, it can be seen that the solution is not modified when using higher values. Furthermore, the Green-Lagrange strain component E zz is plotted for point A (x = 0.73 mm, y = 0.73 mm, z = 0.9 mm) of a broken cell at every load step. It can be observed that using the α-method with q = 0 and q = 1 the solution is affected largely for such a local quantity, see Fig. 16 (bottom). On the other hand, utilizing the eigenvalue stabilization we can reach a much further deformation state without affecting the solution largely, as can be seen in Fig. 16 (top).
Next, we also investigate the strain energy using an ansatz order of p = 3, as illustrated in Fig. 17. It can be observed that using the α-method with q = 0 and q = 1 results in a high loss of accuracy. Also, the structure becomes too stiff. However, for q values between 3 and 5, the α-method produces reasonable results with a maximum deformation ofū z = 1.95 mm using q = 3, as can be seen in Fig. 17 (bottom). Utilizing the eigenvalue stabilization, we can reach a further deformation state with a displacement ofū z = 5.75 mm with q = 4, which is about 2.95 times more as compared to the α-method, as can be seen in Fig. 17 (top).
These results show that the strategy of stabilizing the FCM for hyperelastic problems using the technique based on eigenvalues is more robust and accurate as compared to the α-method. In Fig. 18, the von Mises stress of the single cube connector is plotted at different load steps using the ansatz order of p = 3 and = 1e−4.

Single pore of a foam
In the last application, we analyze the performance of the proposed eigenvalue stabilization for a more complex structure such as a pore of a foam [21]. The geometry of the foam is obtained from a CT-scan and then converted into a triangulated surface, as depicted in Fig. 19. For the computation, the bottom surface is fixed in all directions, while the top surface is fixed in x and z-directions with a prescribed displacement in the y-direction to compress the foam. The geometry is discretized with 1815 cells (20 × 20 × 20 subdivisions). The numerical integration is done the same way as in the previous example (Sect. 5.2). We apply displacement increments of 0.02 mm at each load step and try to compress the geometry as much as possible. To this end, we set the ansatz order to p = 2 and apply different stabilization factors for the αmethod with q = {0, 1, 2, 4, 5}, as well as for the eigenvalue stabilization with q = {4, 6, 8}.
To this end, the strain energy function is plotted for the α-method at different load steps, as shown in Fig. 20 (bot- . It can be observed that using q = 0 and q = 1 results in a drastic change in the solution. For q = 4 and q = 5, the accuracy obtained is reasonable. However, the structure could only be deformed up to a displacement ofū y = 1.2 mm using q = 4. On the other hand, utilizing the eigenvalue stabilization, we can reach a much more pronounced deformation state with a displacement ofū y = 3.82 mm with q = 4, which is about 3.2 times more as compared to the α-method. Furthermore, we observe no high loss of accuracy when the values increased from q = 8 to q = 4, as can be seen in Fig. 20 (top).
Furthermore, the Green-Lagrange strain component E yy is plotted for point A (x = 1.9 mm, y = 3.45 mm, z = 4.9 mm) of a broken cell at every load step. It can be observed that using the α-method with q = 4 and q = 5 does not affect the solution largely. However the simulations failed to converge much earlier as compared to the eigenvalue stabilization, as can be seen in Fig. 21. Finally, the von Mises stress for the single pore of a foam is plotted in Fig. 22 at different load steps for = 1e−4. Fig. 19 Single pore of a foam. Geometry, boundary conditions, and mesh

Conclusions
In this paper, we proposed an eigenvalue stabilization technique for the finite cell method in order to improve its robustness for nonlinear analysis at finite strains. The approach, originally proposed by Loehnert [29], is based on the eigenvalue decomposition of the stiffness matrix for cells that are cut by the boundary of the geometry. To this end, eigenvalues that have a smaller value than a certain tolerance are grouped together. Due to the linear dependencies between the modes, the eigenvalues of some of the modes can even be close to zero-which is why they need additional support. The zero modes that are related to rigid body translations and rotations are not considered for stabilization. Afterwards, the bad modes are supported by adding extra stiffness to the corresponding cell stiffness matrix. Additionally, a correction term is added to the load vector so as to avoid a large modification of the system. Furthermore, we proposed an adaptive stabilization criterion so that the modes with the worst eigenvalues receive more support than modes with better support. We investigated the approach using different numerical examples, comparing the results to the commonly used αmethod. To this end, we showed for the linear analysis that the condition number can be improved significantly without affecting the solution at all, provided that the load vector is corrected accordingly. Moreover, we showed that the method is suitable to improve the robustness of the FCM significantly for nonlinear analysis.
There are several aspects which are of interest for future work. One important point is to extend the proposed stabilization scheme to distinguish between numerical and geometrical or physical instabilities. The stabilization should be applied only to instabilities arising from badly broken cells. Geometrical instabilities related to local buckling or physical instabilities caused by the material behavior should not be affected by the proposed scheme. Furthermore, the stabilization scheme could be combined with iterative solvers acting as a kind of preconditioner for nonlinear problems to be solved with the FCM.
Acknowledgements The authors gratefully acknowledge support by the Deutsche Forschungsgemeinschaft in the Priority Program 1748 "Reliable simulation techniques in solid mechanics. Development of non-standard discretization methods, mechanical and mathematical analysis" under the project DU 405/8-2.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.