1 Introduction

The finite cell method (FCM) ^[1,2] can be regarded as an embedding or fictitious domain method ^[3,4,5] combined with the p-version finite element method (p-FEM) ^[6] . The efficiency of the p-FEM for problems of elasto-plasticity has been thoroughly investigated ^{[7,8,9,10,11,12]} . The p-version yields clearly superior accuracy in comparison with its h-version counterparts. Although the FCM has been demonstrated to be an efficient method for linear problems, little work has been done to explore its merits in the elasto-plastic structural analysis ^[13,14,15] . In the FCM,the polynomial degree of the finite element approximation is increased (p-version) to improve the accuracy of the solution,rather than following the traditional mesh refinement approach (h-version). Furthermore,different local refinement strategies have been developed ^[16,17] for different types of discontinuities and singularities to improve the quality of the FCM approximation.

The aim of this paper is to numerically compare the accuracy of the FCM for the PrandtlReuss flow theory of plasticity with that of the h-version FEM. The computations are based on the elasto-plastic J₂ flow theory,using the von Mises yield criterion. Further assumptions include small displacements and small strains. The results show that the FCM makes the analysis of complex geometry problems very convenient and easy as well as reliable. The accuracy of stresses for nearly incompressible materials is always of interest. Incompressibility effects may be introduced by the elasto-plastic material law and may slow down the convergence. The numerical examples presented in this paper demonstrate that a good convergence of stress components can be obtained using the FCM,even if the point of interest is located within the plastic region.

The paper is organized as follows. Section 2 explains the FCM for the elasto-plastic formulation and a modified quadtree integration scheme. Section 3 presents numerical examples to evaluate efficiency of the FCM compared to the h-FEM for the J₂ plasticity,which is followed by the summary and conclusions. 2 FCM 2.1 Elasto-plastic formulation

The FCM is a combination of high-order FEM and fictitious domain methods. Therefore, we use the well-established terminology of the FEM in this paper. To solve the weak form of the equilibrium equation using the high-order FCM,the domain of interest Ω,as shown in Fig. 1, is extended into a domain Ω_e,i.e.,,only with the advantage that Ω_e can be discretized with a Cartesian structured grid of cells.

The weak formulation of the equilibrium conditions of a solid body reads as finding the displacement field u which satisfies the (homogeneous) Dirichlet boundary conditions such that

in which b is the vector of volume loads,t is the vector of prescribed tractions on the Neumann boundary ΓN,the operator ‘:’ denotes the double-contracted product,and v is the test function. α is defined as In each cell,the displacement field u is approximated using hierarchical shape functions based on the integrated Legendre polynomials ^[18] . Hierarchic bases have been used to implement two different types of Ansatz spaces in two dimensions ^[6,19] : the trunk space and the tensor product space . Three different types of Ansatz spaces have been im plemented in three dimensions ^[6,19] : the trunk space ,the tensor product space ,and the space The physically nonlinear problem under contemplation is the J₂ flow theory for small strains with nonlinear isotropic hardening ^[20] . The strains are small and can be decomposed into an elastic and a plastic component. The stresses satisfy both the constitutive equations and the equilibrium conditions. Since the nonlinear stress(σ)-strain(ε) relationship depends on the load history,the weak formulation is a nonlinear functional to be solved incrementally. Using the Newton-Raphson method for each step [t_(n),t_(n+1)] to linearize the weak formulation,it is possible to compute a sequence of equilibrated load steps, where denotes a variable at the ith Newton-Raphson iteration during the load step in [t_(n),t_(n+1)]. The actual stress state as well as the algorithmic tangent modulus have to be recomputed at each step i ^[20] . In linear elasticity problems,the stress(σ)-strain(ε) relationship is linear. The tangent modulus tensor in Eq. (3) denotes the elasticity tensor,and the equation is solved in one step. The spatial discretization of Eq. (3) is performed using the FCM. 2.2 Integration scheme

In the embedded FEM,the boundary of the physical domain is not captured by the mesh. Rather,accurate integration schemes may capture the boundary. One of such schemes used for the FCM is a modified quadtree integration scheme. Quadtree is a tree data structure used for partitioning a domain by recursively subdividing it into four quadrants. The extension of a quadtree for 3D space is called octree,in which each internal node or domain has eight children ^[21] .

With the help of a quadtree or an octree,it is easy to perform the appropriate numerical integration for Eq. (3) using the composed integration over subcells (see Fig. 2).

A mapping has to be used to establish a relationship between the coordinates of the subcell (r,s) and the cell (ξ,η) ^[2] . In this investigation,we use a modification of the standard quadtree integration scheme which has two specific features (see Fig. 3).

(i) Since the subcells are much smaller than the cell,the function over the subcells can be regarded as a lower-order function. It is accordingly possible to use a low-order Gauss quadrature. In the standard p-FEM,p + 1 Gauss points are used in each direction for the elements of order p ^[8] . We therefore begin with p + 1 Gauss points in each direction. Then,we reduce the number of Gauss points incrementally down to 1 following the refining of the quadtree. Any integration point outside the boundary of the physical domain defined by the quadtree is ignored. Figure 3(a) shows the position of active Gauss points of subcells mapped from the local coordinate (r,s) to the parent cell.

(ii) The second feature involves avoiding the ill-conditioning inherent in fictitious methods by assigning a small α to the Gauss points outside the physical domain ^[1] . Therefore,the cell is reconsidered with p + 1 Gauss points in each direction,and the points outside the physical domain are only taken into account to avoid ill-conditioning of the problem (see Fig. 3(b)). As an example,144 integration points are used for integrating over a cell of order p = 4 with 3 levels of quadtree refinement (R = 3),as shown in Fig. 3. 3 Numerical examples

In each example,the unit of dimensions and displacements is mm,and the unit of stresses is MPa. The material is assumed to be elastic/perfectly-plastic with a shear modulus of µ = 80.193 8 MPa,the bulk modulus κ = 164.206 MPa,and the yield stress σ yield = 450 MPa ^[22] .

The plane strain situation and a unit thickness are assumed for 2D problems. The computations for FCM are conducted using the p-FEM code AdhoC ^[23] ,and the direct solver SPOOLES ^[24] is used to solve the overall equation system. The results are visualised by PARAVIEW ^[25] using a fine post-processing mesh. In this paper,we restrict our investigations to a quadtree/octree refinement level of = 3,and the trunk space and The FCM leaves the mesh unchanged and progressively increases the degree p of the piecewise polynomial approximation,until some desired level of precision is reached. In the standard h-FEM,however,the mesh is refined and the polynomial degree of the shape functions remains unchanged.

The h-version mesh generation and computations were carried out with the commercial finite element package ABAQUS ^[26] which uses special pressure-displacement elements for incompressible materials,based on a mixed finite element formulation. These elements also include an additional and internal mode trial function ^[27] . Since it is difficult to create a structured mesh when faced with complex geometries,we use the unstructured mesh generating scheme. We choose the 2D and 3D enhanced low-order approximation elements CPE4R (p = 1) and C3D8R (p = 1) in the h-FEM. The reduced integration scheme and a dense mesh are used to avoid locking and hour-glassing phenomena ^[27] . 3.1 Perforated plate

The first numerical example to be considered is a square plate with a central perforation under cyclic loading. A quarter of the square plate is modelled,and the symmetry conditions are imposed on the lower and right sides of the plate (see Fig. 4). The uniform traction t_n = 100 is scaled with a factor −4.5 6 λ 6 4.5 in 205 load steps. As a reference solution,we use the results provided by Wieners ^[22] based on a very fine mesh with Q2P1 elements. Results of interest are the stress σ_yy at point A = (9,0),the stress σ_xx at point B = (6.195 3,3.804 7), the displacement u_x at point C = (0,10),and u_y at point D = (10,10).

The discretization of h-FEM and the p-version of the FCM are adjusted to achieve almost the same degrees of freedom for comparison. Figure 5 shows 100 cells used with p = 10,which corresponds to 9 600 degrees of freedom,and 4 686 elements of h-FEM,which corresponds to 9622 degrees of freedom.

As shown in Figs. 6−9,the displacements are approximated accurately by both methods, but the stresses are more accurate using the FCM. Let us consider this problem again with a load factor λ = 4.5 and study the convergence behaviour in terms of degrees of freedom. The first and last meshes used for the h-FEM are shown in Fig. 10. For FCM calculations,the polynomial degree varies from 4 to 14 using 100 cells (see Fig. 5). The results are plotted in the form of relative errors in a double logarithmic scale in Figs. 11 and Fig. 12. As these figures indicate,the FCM yields smoother and more accurate results particularly at point A,where the most critical stress concentration occurs. Using the FCM,only 4 740 degrees of freedom with a corresponding polynomial degree p = 7 are needed to compute the stress components with a relative error less than 5%. Considering these numerical investigations,it is apparent that it is easy to obtain a degree of accuracy suitable for practical implementation by means of the FCM.

3.2 Perforated thick plate t_n =100

For the second example,we consider a thick perforated plate with a unit thickness under the uniform traction of t_n = 100 that is raised monotonously to a factor of λ = 4 in 49 load steps (see Fig. 13). The symmetrical boundary conditions are used at the planes of y = 0,z = 0,and x = 10. By a reference solution,we use the p-version FEM with 199 hexahedral elements of p = 10 (see Fig. 14),whereas the curved boundary has been represented precisely with blending functions ^[6] . The result of interest is the stress σ_yy at the critical point E = (9,0,0),where the stress concentration occurs. 100 cells with p = 8 are used for the FCM calculations,which in this case correspond to 30 468 degrees of freedom. For the h-version calculations,12 445 elements give 31 092 degrees of freedom (see Fig. 15).

Considering the results in Fig. 16,the FCM approach for p = 8 conforms better to the reference solution than the h-FEM. It is evident that the FCM provides a more accurate approximation using higher order p. This shows that the most accurate results can be achieved by the FCM.

A three-dimensional view of the von Mises stress for the FCM with p = 8 and the p-version reference solution at the last load step λ = 4 is shown in Fig. 17.

Since two different codes,ABAQUS and AdhoC,are used in the present study,it is difficult to achieve a fair comparison for the CPU computational time. It is,however,safe to conclude that it is faster to pursue solution convergence in the FCM,because it has a fixed mesh and the polynomial order is only increased for the sake of a more accurate approximate solution. It is also possible to use a fixed dense pattern of integration points for a p-extension,allowing the entire pre-processing procedure including generating the Cartesian mesh and setting up the integration points to be performed at the beginning of the computation. Besides,the FCM has all the advantages of the p-version FEM,such as condensing interior degrees of freedom which are confined to the element. The elimination of interior degrees of freedom can be performed in parallel on element level without any communication. This increases the efficiency of the overall parallel approach,which takes advantage of the fact that the computation of element matrices of high order constitutes a significant part of the overall numerical effort. Furthermore,it can be shown that the elimination of the interior degrees of freedom improves also the conditioning of the resulting global equation system which is important when using iterative solvers ^{[28,29,30,31]} . 4 Conclusions

This paper compares the FCM numerically with the h-version FEM for the Prandtl-Reuss flow theory of plasticity. The results demonstrate the efficiency of the FCM to solve materially nonlinear problems. Numerical examples show that the FCM provides efficient and accurate approximations to this class of physically nonlinear problems. It is easier to control the discretization error with the FCM than that with the h-FEM. Therefore,the degree of error can be reduced by increasing the order of elements and not necessarily by a new mesh. Error monitoring and control are achieved automatically,with minimum user interaction and without substantial loss of computational efficiency. For two- and three-dimensional benchmark problems,we find that the models based on the FCM deliver much more accurate results compared with their h-version counterparts. To summarise,the application of the FCM in elasto-plastic structural analysis can be recommended. These observations suggest that it would be beneficial to extend the use of the FCM to problems of elastoplasticity with finite strains. Results in this context will be presented in forthcoming papers.

Acknowledgements   This work is the result of an institutional partnership being supported by the Alexander von Humboldt Foundation. This support is gratefully acknowledged.