Isogeometric Analysis of Soil Plasticity

In this paper we present numerical simulations of soil plasticity using isogeometric analysis comparing the results to the solutions from conventional finite element method. Isogeometric analysis is a numerical method that uses nonuniform rational B-splines (NURBS) as basis functions instead of the Lagrangian polynomials often used in the finite element method. These functions have a higher-order of continuity, making it possible to represent complex geometries exactly. After a brief outline of the theory behind the isogeometric concept, we give a presentation of the constitutive equations, used to simulate the soil behavior in this work. The paper concludes with numerical examples in twoand three-dimensions, which assess the accuracy of isogeometric analysis for simulations of soil behavior. The numerical examples presented show, that for drained soils, the results from isogeometric analysis are overall in good agreement with the conventional finite element method in twoand three-dimensions. Thus isogeometric analysis is a good alternative to conventional finite element analysis for simulations of soil behavior.


Introduction
In the design of foundations and geotechnical structures it is essential to predict soil behavior under different loading conditions. In the last decade finite element analysis (FEA) has become a widely spread tool for predicting soil behavior.
Much research has been carried out to improve the ability of simulating the behavior of different soils and a number of new constitutive models have been developed. Another important aspect when modeling geotechnical problems is the interaction between soil and structure, which can have a large influence on the structural design [1].
Since first introduced by Hughes et al. in [2], isogeometric analysis (IGA) has been used for a broad number of engineering problems, ranging from analysis of fluid-structure interaction [3] to analysis of shells [4]. Some of the areas where isogeometric analysis has shown to be advantageous compared to the standard finite element analysis are also of interest in geotechnics. One example is the simulation of fluid flow through porous media, where the higher-order basis functions of IGA lead to continuous pressure gradients over element boundaries ensuring local mass conservation, [5] and [6]. The basic idea behind isogeometric analysis is to use splines as basis functions for computational analysis [2]. Originally the purpose was to use the geometry from the Computer Aided Geometric Design (CAGD) software's also for analysis and consequently reduce the time required to recreate models for analysis. However, the use of spline basis functions proved to have numerical benefits compared to the standard Lagrangian basis-functions used in FEA [7] [8]. Initially the IGA framework suffered from the lack of a method for local refinement, frequently sought during analysis. To overcome this initial drawback, a number of methods have been developed to introduce local refinements within the IGA framework, where T-splines [9] or Hierarchical NURBS [7] [10] are the most widely used methods. The higherorder smoothness of the basis functions in IGA also opens up the possibilities of using quadrature rules evaluated on element boundaries reducing the total number of quadrature points [2] [7]. Another interesting possibility within the IGA framework is the use of isogeometric collocation methods, a one-point quadrature rule that reduces the computational cost of analysis [11]. Since collocation methods are based on the discretization of the strong form of the partial differential equations, it would be suitable for the isogeometric concept. The use of isogeometric collocation has shown potential to increase the computational efficiency of the isogeometric framework, outperforming both the standard Galerkin form of isogeometric analysis and standard C 0 finite element analysis with regard to computational costs [7]. As pointed out by Nguyen et al. [12] for nonlinear analysis the numerical error strongly depends on the integration scheme and the order of the Non-Uniform Rational B-splines (NURBS) basis functions.
To evaluate the isogeometric framework for geotechnical applications, it is interesting to assess how the isogeometric framework performs for soil plasticity.
For simplicity and for setting up the numerical framework we restrict our self to quadratic NURBS basis functions in this work. As soil typically show strong plastic flow under hydrostatic pressure, special B-bar methods are not necessary in our work, see [13] for a discussion of these methods in the case of nearly incompressible material behavior.
In this work we have evaluated how isogeometric analysis performs compared to the conventional finite element method for soil plasticity in two and three dimensions. The focus of the work is to evaluate the convergence behavior and mesh size dependencies for isogeometric analysis and to compare them with results from conventional FEA. To provide a background, we start with outlining the basic concepts of isogeometric analysis, introducing the definitions of B-splines and NURBS that are used as basis or shape functions in this work. We later continue with the basic governing equations for elasto-plasticity followed by a brief recapitulation of the Drucker-Prager constitutive equations. The paper is concluded with two-and three-dimensional numerical examples, comparing the results to conventional finite element analysis.

Isogeometric Analysis
The fundamental idea behind IGA is to employ the same basis functions for both geometrical discretization and analysis. Herein lays also the most profound difference between IGA and standard FEA, where isogeometric analysis utilizes the basis functions from CAGD, capable of representing the exact geometry also for analysis. Whereas, in conventional finite element analysis, the piecewise polynomials chosen to approximate the solution fields are also used to approximate the geometry [2]. The original reason to use the isogeometric analysis was to decrease the overall computational cost by removing the need to recreate the geometry and constructing a mesh for analysis. In this work we have analyzed the performance of IGA for soil plasticity using a NURBS-based isogeometric formulation. To give a brief introduction to the concept of IGA and to elucidate some of the differences between conventional FEA and IGA this section will review the basic concepts of isogeometric analysis. For a more extensive description the reader is referred to [2].

Basic Concept of B-Splines
To get an overview of the concept of NURBS-based isogeometric framework which is referred to as the Cox-de Boor formula, first presented in [14] and [15].
One example of a second order B-spline with three control points is illustrated in Figure 1. For polynomial degrees p < 2 the basis functions of IGA and conventional FEA take the same shape, however for p ≥ 2 they deviate in shape and support. By plotting the B-spline basis functions and Lagrangian basis functions over the parameter space some of the differences between IGA and conventional FEA are made visible, see Figure 2 and Figure 3. Figure 3 shows that for p ≥ 2, Lagrangian basis functions, the continuity varies between internal and corner-or end-nodes, whereas B-spline basis functions show a continuous and homogeneous pattern only shifted relative to each other and with N i,p (ξ) ≥ 0. Another important aspect of spline basis functions, that distinguishes them from conventional FEA basis functions, is that each p th order function has p − 1 continuous derivatives across element boundaries.

NURBS Representation
One issue that arises when using B-spline, is that not all types of geometries can be represented exactly using polynomial functions. To overcome this problem rational B-splines were introduced by Versprille in [16]. This generalization of Bsplines is constructed by introducing a non-negative weight, w i , to each control point and making use of the definition of rational functions as the ratio of two polynomials [17].
where W(ξ) is called a weighting function, defined as where all the weights take the same value. Using the NURBS basis functions and their corresponding control points P i , a piecewise NURBS curve, is constructed from To be able to perform two-and three-dimensional analysis, NURBS surfaces and bodies needs to be defined. This is done in a similar manner. A NURBS body is given by where the NURBS basis function,

( )
, , , , , , p q r i j k R ξ η ζ , is constructed from three sets of knot vectors, and their respective weighting functions In the isogeometric framework the concept of elements is represented by the non-zero valued knot spans. This is illustrated in Figure 4, where a two dimensional plate is constructed from two knot vectors, a set of control points P i,j and their corresponding weights w i,j . The geometry is constructed of one knot vector,

Strong and Weak Form
In this section we give a brief recapitulation of the strong and weak form of the equilibrium equations of the quasi-static balance of linear momentum. Let u i denote the displacement vector, then the infinitesimal strain tensor, ε ij , is defined as the symmetric part of the displacement gradient The governing strong form of the equilibrium equation is given as where f i is the body force. The strong form of the equilibrium equation is complemented by the essential and natural boundary conditions.
where g i and h i are known quantities on ∂Ω g and ∂Ω h . The variational or weak form of Equation (11) is formed by multiplying the governing equation with a test function v i and performing an integration by parts over the domain Ω The weak form of the problem is complemented by a constitutive relation with { } κ being a set of internal variables to describe the plastic behavior of the material. For a more detailed derivation of the weak form, the reader is tentatively referred to any of [18] or [19].

Discretization with Isogeometric Analysis
The main difference between conventional FEA and IGA are the basis functions used for discretization. In isogeometric analysis the same basis functions that are used to discretize the geometry are also used to solve for the ap- In a similar manner we rewrite Equation (17) where the matrix B e is an operator mapping the element discrete displacements to the local strains.
The test function v i can be discretized in the same manner as the displacement field. Introducing the relations above it is possible to rewrite Equation . Afterwards, we update the displacement vector by a ← a + ∆a and stop the iteration if the residual r(a) < TOL with TOL a user given tolerance, e.g. TOL = 10 −6 .

Constitutive Model (Heading 4)
In this section we will give a brief recapitulation of the Drucker-Prager criterion and the theory of plasticity used to evaluate the influence of the NURBS-basis functions on the plastic strains for granular materials in this work.

Drucker-Prager Criterion
The Drucker-Prager criterion can be seen as a smooth approximation to the

Elasto-Plastic Constitutive Model
We assume a small strain setting and thus additively split the strain tensor ε into an elastic part e ε and a plastic part p ε as For geo-materials in general, associative flow rules for the plastic strain p ε displays an excessive dilatant behavior. This can be avoided by using a non- associative flow rule, that is, the potential function  is not equal to the yield function  . In this work the potential function  is defined by replacing the angle of internal friction, φ′ , with the angle of dilation, ψ, in Equation (25) forming sin and . sin sin tan cos 3 The potential function is then written as ( ) , see also Equation (15). For the numerical integration of the constitutive evolution equations we use a standard Euler backward scheme. The reader is referred to ([20], Sec. 7.2) which gives a detailed description of the numerical implementation procedure and also the format of the algorithmic tangent stiffness D ats . We like to mention that at the apex of the yield surface cone, the return vector is contained by a complementary cone, illustrated in Figure 7, see [20] for further details. The IGA and FEA as well as the elasto-plastic Drucker-Prager model discussed above have been implemented in our in-house Fortran code.

Numerical Studies
To validate the performance of isogeometric analysis for soil-plasticity, three numerical benchmark models have been established. The models have all been simulated for soils in saturated conditions using the Drucker-Prager criterion.
The models evaluated consist of one two-dimensional model of a strip footing and two three-dimensional cylindrical soil profiles subjected to a prescribed force and prescribed displacement respectively.

Two-Dimensional Study
A two-dimensional model of a strip footing on sandy silt has been analyzed. Due to the symmetric properties of the problem the analysis contains only half of the footing. The geometry and boundary conditions are shown in Figure 8. To model the load from a flexible footing a vertical pressure of 180 kPa is applied in 50 equal steps.
The problem is solved for plane strain conditions using quadratic NURBSbased IGA and conventional FEA with 5 different meshes. In order to compare the two methods the element meshes have been constructed using quadratic isoparametric elements for both IGA and conventional FEA. The element size ranges from a coarse mesh with elements of 2 × 2 meters down to a mesh with elements of 0.25 × 0.25 meters. The mesh data and degrees of freedom for each mesh are shown in Table 1. The material parameters of the soil in the two dimensional model are shown in Table 2.     Figure 9. Ground surface displacements underneath along the loading surface.
that there is a minor difference between the displacements from the isogeometric-and conventional finite element analysis at the edge of the footing (x = 2). Figure 10 displays the load/displacement response in point A and B for the finest mesh. By studying the load/displacement response in Figure 10, the differences between the IGA and FEA results in point B become more apparent.
The results for point B showed in Figure 10 are of special interest. The NURBS basis functions used in the isogeometric analysis provides continuous stress fields at the edge of the footing, where discontinuous stresses can be expected.
The continuous stress will in turn affect the development of plastic strains at the edge of the footing.
The effects of the continuous stress fields that result from the NURBS basis functions can be seen in Figure 11, where the evolution of the plastic zone under the footing is displayed. From load steps 20 -22 in Figure 11 it can be seen that the initial plastic strains are slightly different at the edge of the footing in the conventional finite element analysis compared to the isogeometric analysis. The difference in plastic strains at the edge of the footing corresponds to the diverging displacements seen in Figure 9 and Figure    results from the finest mesh used. The graphs show that the convergence for the isogeometric analysis and the conventional finite element analysis are comparable.

Three-Dimensional Studies
The two three-dimensional studies in this work are composed of a cylindrical soil profile with a 2:1 height/diameter proportion. In the first example, the cylindrical soil profile is subjected to a prescribed displacement at the top of the cylinder, in the second study the soil profile is subjected to a confining pressure and an increased vertical load at the top of the cylinder. The geometry and boundary conditions of the three-dimensional examples are presented in Figure 13.
In both examples, a perfect elasto-plastic Drucker-Prager criterion, matched to the compressive meridian of the Mohr-Coulomb criterion has been used to model soil plasticity. The material properties for the soil in both studies are given in Table 3. For both IGA and conventional FEA 27 node isoparametric elements have been used to model the cylinder. The mesh data used are presented in Figure 14.
Using NURBS parameterization an exact cylinder can be modeled. Whereas, for conventional FEA the geometry depends on the density of the mesh.  Figure 13. Boundary conditions for the three-dimensional models. The displacement-controlled analysis is fixed in the bottom of the cylinder and has no horizontal constraints at the top (a). In the load controlled analysis the bottom of the cylinder is assumed to be fixed (b). Figure 14. The left hand side shows a horizontal view of the finite element mesh and on the right hand side shows the conventional FEA (top) and IGA (bottom) meshes viewed from the top. The figure also includes mesh data used in the simulations.

Displacement Controlled Analysis
In the displacement controlled analysis a total strain of 7% is applied over 100 equal time by prescribing a displacement at the top surface. To compare the results from the isogeometric analysis with the finite element analysis, Figure 15 shows the applied pressure normalized to the cohesion against the axial strains.
To study the effects of the mesh density, Figure 16 shows the average vertical  stress components over the top surface of the soil profile for each mesh normalized to the vertical stress from the finest mesh. From Figure 16 it is clear that the deviation of the vertical stresses is slightly less pronounced in the results from the isogeometric analysis compared to the conventional finite element me-thod. Comparing the results from the two dimensional example, the authors judge that the difference arises from the beneficial geometrical parameterization using NURBS in the isogeometric analysis.

Force Controlled Analysis
In the force-controlled example the soil profile is subjected to a confining pressure, σ c , and an axial load, σ v , acting on the top surface as illustrated in Figure   13(b). The analysis is run over 100 load steps. The load is applied by first increasing the confinement pressure to 3 kPa after which the vertical load is increased to 14 kPa. To compare the results from the isogeometric analysis to the results from the conventional finite element analysis, the applied force, normalized to the cohesion, is plotted against the axial strains in Figure 17. The results show that the solutions are in good agreement but that there is a difference between the displacements from the isogeometric analysis compared to the finite element results. Comparing the deviation of the displacement at the top surface for the load controlled example, the influence of using the exact geometry in the isogeometric analysis becomes more evident than for the displacement-controlled analysis, see Figure 18 and Figure 16. This would be expected as the confining pressure σ c will be more accurate using a soil profile made out of a perfect cylinder, than one discretized with Lagrangian polynomials. This can be seen in Figure 18 where the results from the coarser meshes clearly are in error.

Summary and Conclusion
The aim of this work is to evaluate the isogeometric framework for numerical