S.Morganti,F.Fahrendorf,L.De Lorenzis,J.A.Evans,T.J.R.Hughesand A.Reali

1Department of Electrical,Computer and Biomedical Engineering,Università degli Studi di Pavia,Pavia,27100,Italy

2Institute of Applied Mechanics,Technische Universität Braunschweig,Braunschweig,38106,Germany

3Department of Mechanical and Process Engineering,Zurich,8092,Switzerland

4Department of Aerospace Engineering Sciences,University of Colorado Boulder,Colorado,80309,USA

5Institute for Computational Engineering and Sciences,The University of Texas at Austin,Texas,78712,USA

6Department of Civil Engineering and Architecture,Università degli Studi di Pavia,Pavia,27100,Italy

ABSTRACT We investigate primal and mixed u-p isogeometric collocation methods for application to nearly-incompressible isotropic elasticity.The primal method employs Navier’s equations in terms of the displacement unknowns,and the mixed method employs both displacement and pressure unknowns.As benchmarks for what might be considered acceptable accuracy,we employ constant-pressure Abaqus finite elements that are widely used in engineering applications.As a basis of comparisons,we present results for compressible elasticity.All the methods were completely satisfactory for the compressible case.However,results for low-degree primal methods exhibited displacement locking and in general deteriorated in the nearly-incompressible case.The results for the mixed methods behaved very well for two of the problems we studied,achieving levels of accuracy very similar to those for the compressible case.The third problem,which we consider a“torture test”presented a more complex story for the mixed methods in the nearly-incompressible case.

KEYWORDS Isogeometric analysis;isogeometric collocation;nearly-incompressible elasticity

1 Introduction

The problem of nearly-incompressible linear elasticity and the limit problem of incompressibility,governed by the Stokes equations,play important roles in engineering analysis.The equations of nearly-incompressible elasticity are pertinent to the analysis of rubber materials,and the Stokes problem is the standard model for slow viscous fluid flow and represents an important modeling step toward the development of computational formulations of the full Navier-Stokes equations.Near-incompressibility also arises in metal plasticity applications in which the elastic deformations are compressible,but plastic flow is modeled as incompressible.The limit problem of incompressibility is appropriate for the elastic-plastic analysis of undrained soils.Given the physical importance of the nearly-incompressible model,it has been widely studied in the finite element literature and a variety of approaches have been developed to achieve numerical stability and overcome “mesh” locking,see,e.g.,[1].

Recent works of Elguedj et al.[2,3] have investigated NURBS-based isogeometric Galerkin formulations based on-bar and-bar projection methods of nearly-incompressible small-and large-deformation elasticity and plasticity problems.Success was achieved by employing displacement fields one polynomial degree higher and one degree smoother than the volumetric strain field.The family of elements studied were thus maximally smooth isogeometric “QpQp-1elements”forp=1,2,3,4.Thus the classicalQ1Q0element,first introduced in [4],often referred to as the“constant pressure”element and denotedq1p0,was the lowest-order member of the family studied.Due to the widespread use of this element in engineering applications,it proved very useful as a benchmark against which the behavior of the higher-order approaches could be measured.

In addition to Galerkin formulations,Isogeometric Analysis is also amenable to discretization by collocation methods due to the increased smoothness of basis functions.The primary motivation for developing collocation methods is computational efficiency,as they minimize the number of quadrature points to one per node,independent of polynomial order.The study of isogeometric collocation methods was initiated in [5] and further developed in [6–11].The efficiency of isogeometric collocation methods,compared to finite element and isogeometric Galerkin methods,was studied in Schillinger et al.[12] and in De Lorenzis et al.[13].An observation made in this work was that collocation has advantages that increase significantly with the degree of the basis.It may also be noted that the cost of formation of element arrays for the standard Gaussian quadrature rules,involvingO(p3) quadrature points per element in three dimensions,isO(p9),which is prohibitively expensive for higher-degree elements.For collocation it isO(p3),the optimal result,which becomes decisive at higher degree.However,the Galerkin method has always been the gold standard as far as accuracy is concerned and it is difficult to make sweeping generalizations about it compared with collocation.This is especially true in light of recent developments that considerably speed up the isogeometric Galerkin method by utilizing weighted quadrature,sum factorization and row or column formation and assembly attainingO(p4) cost;see Hiemstra et al.[14],or reduced quadrature schemes which require only two points per parametric direction regardless of the discretization order by exploiting the concept of variational collocation [15,16].

Mixed methods constitute perhaps the most important finite element technology for addressing constrained media problems,such as the equations of nearly-incompressible elasticity.However,there have been very few studies of mixed isogeometric collocation methods for problems of this type.In a previous work,the authors considered a mixed isogeometric collocation formulation of nearly-incompressible elasticity and plasticity in which the entire stress field is approximated in addition to the displacement field [17].Models of this type also have relevance to the analysis of viscoelastic fluids in the Eulerian formulation.However,perhaps the most studied mixed formulation of nearly-incompressible elasticity is the one in which pressure is approximated in addition to displacements.This is the primary focus of this paper,which we believe is the first time this model has been investigated from the standpoint of isogeometric collocation.We also consider a so-called primal formulation in which we employ the standard displacement equations unaltered.For the primal formulation we consider maximally smooth NURBS displacement fields of polynomial degreep=2,3,...,7.For the mixed elements,we consider displacement fields of degreep=3,4,...,7 and pressures one degree lower.Herein,we refer to these discretizations asQpQp-1.As a benchmark for acceptable accuracy of a mixed formulation in compressible and nearly-incompressible cases,we presentq1p0-results from the commercial finite element code Abaqus.For sufficiently high-degree collocation results we obviously expect faster convergence results,so this is utilized to provide a context,but not a comparison.

To gauge the performance of the methods,we present results for three linear,isotropic,elasticity problems;two two-dimensional problems and one three-dimensional problem.In all cases the NURBS geometry map is non-affine.We do not present mathematical error estimates.To the best of our knowledge,the only rigorous proofs of stability,convergence and error estimates for collocation methods are those presented in [5] and they are only valid in one dimension.Our computational results in multidimensions were consistent with the mathematical results for one dimension,but the generalization of the mathematical theory to multiple dimensions,to the best of our knowledge,remains an open problem.

The first two-dimensional problem utilizes a plane strain manufactured solution that is divergence-free with zero displacement boundary conditions,i.e.,homogeneous Dirichlet boundary conditions.The body force derived from the manufactured solution is proportional to the shear modulus,μ,and is independent of the other Lamé parameterλ.This problem may be thought of as providing insight into cases that are dominated by shear deformations,such as occur in plasticity.

The second two-dimensional problem is the classical,internally pressurized,plane strain,thick-walled cylinder problem.The exact solution to this problem involves both distortional and dilatational deformations.All boundary conditions are of Neumann type.The applied internal pressure is constant and independent of bothμandλ.

The third problem may be thought as a “torture test.” We employ a three-dimensional domain with homogeneous Dirichlet boundary conditions and assume a manufactured solution that has both distortional and dilatational components of the same order,resulting in body force loading with bothμ-andλ-proportional components.In the nearly-incompressible case,taken asλ/μ=104herein,this results in a very large and dominant pressure loading.The physical relevance of this problem may be questioned,but it does let us explore another aspect of element behavior besides locking in nearly-incompressible cases,namely,“pressure robustness,” a concept introduced by Alexander Linke and investigated by him and his collaborators in a number of recent works;see,e.g.,[18] and references therein.

An outline of the remainder of the paper follows:In Section 2 we describe the collocation methods for linear,isotropic elasticity.In Section 2.1 we present the primal formulation and in Section 2.2 we present the mixed formulation.Numerical results for the two-dimensional problems and the three-dimensional problem are presented in Sections 3 and 4,respectively.Conclusions and a summary of results are presented in Section 5.

2 Isogeometric Collocation for Linear Elasticity

In this section,following the presentation of [19],we apply the ideas of isogeometric collocation to linear elasticity problems.Primal and mixedu-pformulations are described.The problem we consider is represented in Fig.1 and consists of an elastic bodyΩ⊂Rdsubjected to body forces f,to prescribed displacements g on a portion of the boundaryΓD,and to prescribed tractions h on the remaining portionΓN,withΓ=ΓD∪ΓNthe boundary of the domain,and withΓD∩ΓN=Ø.Suitable regularity requirements are assumed to hold for f,g,and h.

Figure 1:Sketch of a generic elastic body Ω subjected to volume forces f,to prescribed displacements g on a portion of the boundary ΓD,and to prescribed tractions h on the remaining portion ΓN

2.1 Primal Formulation

The primal (i.e.,displacement-based) formulation of the small-strain linear elastostatic prob

lem in strong form (Navier’s equations) is given by

complemented by the Dirichlet boundary conditions

and by the Neumann boundary conditions

where u(x) is the unknown displacement field (x being the position vector),∇Sis the symmetric part of the gradient operator ∇,n is the unit outward normal to the boundary of the domain,and C is the fourth-order elasticity tensor defined as follows:

where I and I are the fourth-and second-order identity tensors,respectively,andλandμare the Lamé constants.

In this paper,the collocation approach is applied in the context of Isogeometric Analysis(IGA),and B-splines or NURBS [20] are used to represent both geometry and problem variables in an isoparametric fashion [21,22].

The basic ingredient for the construction of B-spline and NURBS basis functions is theknot vector,i.e.,a set of non-decreasing coordinates in the parameter space:Ξ={ξ1=0,...,ξn+p+1=1},wherepis the degree of the B-spline andnis the number of basis functions.In the present work,we always employ so-calledopenknot vectors,where the first and the last knots have multiplicityp+1.Basis functions formed from open knot vectors are interpolatory at the ends of the parametric interval [0,1] but are not,in general,interpolatory at interior knots.

The construction of the IGA collocation method is obtained following [6] by seeking an approximation uMfor the unknown displacement field u of the elastic problem in the form

to be a set of basis functions defined on:=[0,1]d,dbeing the number of dimensions,such that the physical domainωin (1) can be described by a global geometry function

whereci∈Rdare the control points.We assume that the map G is a bijection.The unknown vectors∈Rdare referred to as displacement coefficients,or control variables.Expression (5) is substituted into Eqs.(1)–(3).

Let us assume for simplicity thatd=2.We denote bym1andm2the number of basis functions in the two parametric directions.ThenM=m1m2is the total number of unknown coefficients per displacement component.We chooseMcollocation pointsτkl,k={1,...,m1},l={1,...,m2}located at the images of the tensor-product Greville abscissae of the knot vectors.1In the three-dimensional setting (d=3),we proceed in a completely analogous way,considering the third parametric direction(e.g.,the M=m1m2m3 collocation points will be indicated as’τklm,with k={1,...,m1},l={1,...,m2},m={1,...,m3}).

The Greville abscissae [23] related to a univariate spline space of degreepand knot vectorΞ={ξ1=0,...,ξn+p+1}are points of the parametric space defined as

The Greville abscissae are simple to compute and have proven effective in many applications.However,there are other possibilities with interesting properties (e.g.,[16,24]).In the two-dimensional primal case,2Mscalar equations are needed to determine the displacement coefficients.In the patch interiorω,we obtain 2(m1-2)(m2-2)scalar equations by collocating Eq.(1) at the pointsτkl,k={2,...,m1-1},l={2,...,m2-1}:

At the Dirichlet boundaryΓDwe impose

To enforce Neumann boundary conditions,Eq.(3) is collocated at the pointsτkl∈ΓNaccording to the following strategy,see [7]:

where nLand nRare the unit outward normals of the two edges meeting at the corner,and hLand hRare the respective imposed tractions.In addition,we refer the reader to [7] for a detailed discussion on the conditions to be imposed in more complicated situations,like at the interfaces of multi-patch geometries.

We note that,as has been shown in [13],the above “basic” approach to imposing Neumann boundary conditions may lead to difficulties in situations when non-uniform meshes are adopted.In such cases,alternative methods for imposing Neumann boundary conditions should be adopted,and in [13] two effective strategies are described.

2.2 Mixed u-p Formulation

A mixedu-p2Throughout,the most popular symbols for both spline degree and pressure are used,namely p;The differentiation between the two is made clear by the context.formulation is readily obtained starting from the equilibrium equations in differential form (1)–(3),written in terms of displacements,and introducing the “pressure-like”variablep=-λ∇· u,yielding:

complemented by the Dirichlet boundary conditions

and the Neumann boundary conditions

We note that in the incompressible limit (λ→+∞)pcorresponds to thehydrostatic pressure.Otherwise,pis simply a scalar field defined by Eq.(13).See [1],Chapter 4,for elaboration.Despite this,in the remainder of the paper,we simply refer topas the “pressure.”

In the mixed formulation we need to represent the pressure field in a similar,but not,in general,identical way as the displacement field.So we assume the pressure field takes the form

Thepiare called the pressure coefficients,or pressure control variables,and theare the corresponding B-spline or NURBS basis functions.Note,as per our assumption,in general,are different than thebasis functions of the displacement field.Expression (16) is now to be substituted into Eqs.(13) and (15).

Denoting byn1andn2the number of pressure basis functions in the two parametric directions,N=n1n2is the total number of unknown coefficients for the pressure.We chooseNcollocation points,j={1,...,n1},h={1,...,n2} located at the images of the tensor-product Greville abscissae of the knot vectors for the pressure field.In this case,2M+Nscalar equations are needed to determine the displacement and pressure coefficients.

In the patch interiorΩ,we obtain 2(m1-2)(m2-2)+n1n2scalar equations by collocating Eq.(13) at the points,and,k={2,...,m1-1},l={2,...,m2-1}andj={1,...,n1},h={1,...,n2}as follows:

The Dirichlet and the Neumann boundaries are treated analogously to the primal case previously described.In particular,to enforce Neumann boundary conditions,Eq.(15) is collocated at the points∈ΓNaccording to the following strategy:

3 Numerical Results in 2D

The isogeometric collocation method is tested on two classical plane-strain elasticity problems:The first one entails homogeneous Dirichlet boundary conditions on the whole boundary,while the second one presents mixed Dirichlet and Neumann boundary conditions.Both benchmarks are implemented for an annular shape to include a non-trivial geometry map between the parametric and the physical domain.

For each problem we consider both compressible (λ/μ=1) and nearly incompressible (λ/μ=104) situations and we report convergence plots for all the unknown fields of theu-pformulation.After presenting the results obtained with the primal collocation formulation,clearly affected by locking in the nearly incompressible regime,the improvements produced by using theu-pcollocation approach are shown.Unequal order approximations for the displacement and pressure fields are investigated.The polynomial order of the shape functions used for the displacement field is taken one order greater than that of the shape functions of the pressure field,i.e.,we useQpQp-1elements.The idea is similar to the Taylor-Hood element of the Galerkin method.In that case both the displacement and pressure areC0-continuous.Here the displacement and pressure areCp-1-continuous andCp-2-continuous,respectively.Formulations of this type have been used successfully in isogeometric Galerkin methods,see [2,3].

In Fig.2 the positions of the collocation points of the displacement and pressure fields are represented for theQ3Q2case (m1=m2=10,n1=n2=9) for the quarter of annulus geometry considered in the following examples.

Additionally,we compare the isogeometric collocation results obtained with finite element solutions obtained using a hybridu-pformulation in Abaqus (Simulia,Dassault Systémes,Providence,USA),namely CPE4H elements,adopting piecewise bilinear displacement and constant pressure approximation.This element is a variant of the classicalq1p0element [4].

Figure 2:Position of the collocation points for the displacement field (blue stars) and pressure field (red circles) for the case of same meshes,Q3Q2 elements;number of displacement collocation points:m1=m2=10,n1=n2=9

3.1 Quarter of an Annulus with Non-Uniform Body Load and Homogeneous Dirichlet Boundary Conditions

For the first benchmark,we refer to [6] and we consider a quarter of an annulus as sketched in Fig.3,with an external radiusRe=4 and an internal radiusRi=1.

The domain is exactly represented by a single NURBS patch and is assumed to be clamped.Following [25],we assign a divergence-free manufactured solution in terms of displacement components,which satisfies the prescribed boundary conditions:

The load is then calculated starting from the manufactured solution.

Figure 3:Quarter of an annulus with non-uniform body load and homogeneous boundary conditions

3.1.1 Primal Formulation

In Fig.4,we report forλ/μ=1 andλ/μ=104the convergence plots of the L2-norm of the displacement erroru*defined as

where the superscriptestands for the exact displacement components,while the superscripthis used to indicate the components of the approximated displacement field.Primal pressures are derived fromp=-λ(∇· u).The exact pressures in this example are zero.Consequently,the pressure errors are reported as the absolute values.

Figure 4:Primal formulation.Convergence plot of the L2-norm of the displacement and pressure error for the quarter of an annulus with non-uniform body load and homogeneous boundary conditions problem for λ/μ=1 and λ/μ=104

In Fig.4,we also report the convergence plots of the L2-norm of the pressure errorp*.

Remark.Throughout,convergence rates are determined from the two points of the convergence plots corresponding to the two finest meshes.

It is clear that,in the nearly incompressible regime,i.e.,whenλ/μ=104,severe volumetric locking appears when using quadratic and cubic basis functions.For higher degrees (i.e.,p＞3),the rates of convergence are the same as those obtained forλ/μ=1,while the absolute errors are two to four orders of magnitude higher.The same observations can be made from Table 1 summarizing displacement and pressure error convergence rates and values.

3.1.2 Mixed u-p Formulation

The results obtained with the mixedu-pcollocation method are reported in terms of convergence plots of the L2-norms of displacement and pressure errors in Fig.5.Note that the exact solution for the pressure ispe=0,therefore,for this example,p*cannot be normalized and represents the absolute error:p*=‖pe-ph‖=‖ph‖,that is the L2-norm of the numerical pressure.

Figure 5:Mixed formulation (QpQp-1).Convergence plot of the L2-norm of the displacement and pressure error for the quarter of an annulus with non-uniform body load and homogeneous boundary conditions problem for λ/μ=1 and λ/μ=104

Fig.5 shows that no difference is observed in the displacement errors when moving from the compressible to the nearly incompressible regime.As also reported in Table 1,moving fromλ/μ=1 toλ/μ=104,very similar error values are obtained using both coarse and fine meshes.

Table 1:Quarter of an annulus with non-uniform body load and homogeneous boundary conditions problem:L2-norm displacement and pressure error convergence rates and values for coarsest and finest meshes

Fig.5 also shows that,while the orders of convergence of the pressure error are not affected by the nearly incompressible regime,their values are higher.We may observe that the errors obtained by the Abaqusq1p0element for displacements are comparable with those of the isogeometric collocationQ3Q2case.However,the pressures of the Abaqusq1p0element are substantially worse both in convergence rate and absolute values.

3.2 Infinitely Long,Pressurized Thick-Walled Cylinder

As a second example,we consider the infinitely long and internally pressurized thick-walled cylinder,studied in [6].Exploiting the symmetry of the problem,we can consider only one quarter of the full domain,as represented in Fig.6.As the cylinder is infinitely long,the problem can be solved as a two-dimensional plane strain problem.As in the previous example,we select an external radiusRe=4 and an internal radiusRi=1.

Figure 6:Infinitely long,pressurized thick-walled cylinder

The boundary conditions,depicted in Fig.6,are

whereτis the unit tangent vector.The exact solution in terms of displacements (in radial and circumferential directions) is

whererindicates the radial position,andThe Young’s modulusEand the Poisson’s ratioνare defined in terms of the Lamé constants as:and

The aim of this test is to check whether mixed boundary conditions have an impact on the results.

3.2.1 Primal Formulation

The convergence plot of the L2-norm of the displacement erroru*rdefined as

whereueranduhrare the exact and approximate radial displacement components,respectively,is reported in Fig.7 forλ/μ=1 andλ/μ=104,along with the convergence plot of the L2-norm of the pressure errorp*,defined as

wherepeandphare the exact and approximate pressures,respectively.In the primal formulation case,pressures are derived from displacements.

Figure 7:Primal.Convergence plot of the L2-norm of the displacement and pressure error for the infinitely long,pressurized thick-walled cylinder problem for λ/μ=1 and λ/μ=104

As expected,for the primal formulation,in the nearly incompressible regime,i.e.,whenλ/μ=104,severe volumetric locking appears when using quadratic and cubic basis functions.Whenpis equal to 4,we observe that coarse meshes lead to significant errors and the obtained convergence rate is suboptimal both for displacements and pressures.For higher degrees (i.e.,p＞4),the rates of convergence are the same as those obtained forλ/μ=1,while the absolute errors are three orders of magnitude higher,as summarized by error values reported in Table 2.

Table 2:Infinitely long,pressurized thick-walled cylinder problem:L2-norm displacement and pressure error convergence rates and values for coarsest and finest meshes

(Continued)

Note:† Abaqus q1p0 element.

Figure 8:Mixed formulation (QpQp-1).Convergence plot of the L2-norm of the displacement and pressure error for the infinitely long,pressurized thick-walled cylinder problem problem

Figure 9:Line plots of displacement with different meshes vs.the radius for the case Q3Q2 and Q4Q3 (see square markers of the convergence plot in the top,λ/μ=104)

Figure 10:Line plots of pressure with different meshes vs.the radius for the case Q3Q2 and Q4Q3(see square markers of the convergence plot in the top,λ/μ=104)

3.2.2 Mixed u-p Formulation

The same trends observed in the previous benchmark are observed also in this case.In particular,optimal behavior in terms of both displacement and pressure errors is recovered in the nearly incompressible regime,as shown in Fig.8.Both displacement and pressure errors listed in Table 2 are very similar when moving fromλ/μ=1 toλ/μ=104.Displacement finite element errors (using again the Abaqusq1p0element) are approximately one order of magnitude smaller than those obtained with isogeometric collocation with theQ3Q2elements,but more than two orders of magnitude higher whenp＞3.Pressure finite element errors are always greater than those obtained with collocation.

To better understand the performance of theu-pmethod in the nearly-incompressible regime,since everything varies only radially,line plots for theQ3Q2andQ4Q3cases of displacement and pressure resultsvs.the radius are reported in Figs.9 and 10,respectively.While the radial displacement solution shows accurate results also for coarse meshes independent of the considered polynomial degree,the pressure results for theQ3Q2case are less accurate for coarse meshes,but quickly improve under mesh refinement.Much improved pressure behavior is observed for theQ4Q3case.For higher order cases there are no discernible errors on the scales of the plots for all meshes.

4 Numerical Results in 3D

We tested theu-pformulation on a 3D example and compared with Abaqus C3D8H elements,which are 8-node linear,hybrid,constant pressure elements.We again refer to this element as theq1p0element.The considered solid depicted in Fig.11 is characterized by the following dimensions:Ri=1,Re=4,H=1,and is fully clamped on its entire boundary.We assign a manufactured solution in terms of displacement components,which satisfy the prescribed boundary conditions,which isnotdivergence free.

The load is then calculated starting from the manufactured solution and imposing equilibrium.We note that in the nearly-incompressible case the pressure loading will be amplified by a factor of 104.We remark that the results obtained behaved similarly to a 2D version of this problem,indicating that the dimension of the problem is not a significant influencer.

Figure 11:Three dimensional quarter of annulus

Primal formulation.In Fig.12 the convergence plots of the L2-norm of the displacement3The L2-norm of the displacement error u* is now defined as

where

u

e

1

,

u

e

2

,and

u

e

3

are the exact displacement components,while

u

h

1

,

u

h

2

,and

u

h

3

are the approximated displacement ones. and pressure error are shown for the primal formulation.As observed in the 2D examples,the displacement error exhibits locking for

p

=2 and

p

=3 when

λ

/

μ

=10

4

,while pressure seems to not be significantly affected by the nearly-incompressible constraint,as also shown by error convergence rates provided in Table 3.In comparison with the two 2D problems previously solved,the present results are qualitatively similar with respect to displacements,clearly manifesting locking in the nearly-incompressible case;however,the pressure results are quite different,indicating only slightly larger errors for the nearly-incompressible case compared with the compressible case,

λ

/

μ

=1.

Figure 12:Primal formulation.Convergence plot of the L2-norm of the displacement and pressure error for the 3D quarter of an annulus with non-uniform body load and homogeneous Dirichlet boundary conditions problem for λ/μ=1 and λ/μ=104

Table 3:3D quarter of an annulus with non-uniform body load and homogeneous Dirichlet boundary conditions problem:L2-norm displacement and pressure error convergence rates and values for coarsest and finest meshes

Mixed u-p formulation.In Fig.13 the convergence plots of the L2-norm of the displacement and pressure error are shown for the mixed formulation.For the mixed collocation formulation,the displacement and pressure results forλ/μ=1 are about the same as for the primal formulation,as reported in Table 3.

Figure 13:Mixed formulation (QpQp-1).Convergence plot of the L2-norm of the displacement and pressure error for the 3D quarter of an annulus with non-uniform body load and homogeneous Dirichlet boundary conditions problem for λ/μ=1 and λ/μ=104

We note though that Abaqusq1p0element produces substantially greater errors for pressure.For the caseλ/μ=104the pressure results are all similar to those atλ/μ=1.However,the displacement results atλ/μ=104are significantly in error and generally worse than for the primal formulation.This was not seen in either of the 2D cases presented earlier.As mention earlier,this behavior is not attributable to the problem being 3D as an analogous version ofthisproblem in 2D behaved similarly.Surface plots obtained by cutting the domain with the plane highlighted in Fig.14 are shown to better understand theu-pmethod results in the nearly incompressible regime.The first component of the displacement solution is shown.The left column of Fig.14 highlights the inaccurate behavior of coarse meshes,while in the right column an improvement upon refinement is observed.Similar results are obtained for the other displacement components.

Figure 14:3D quarter of an annulus with non-uniform body load and homogeneous Dirichlet boundary conditions problem:Surface plots of the first displacement component u1 for different orders and meshes (λ/μ=104)

Normally,we are inclined to attribute inaccurate displacement results for the nearlyincompressible case to locking,but this does not seem to be what is occurring here,in fact,quite the opposite.Note the amplitude of the color bars in each frame and compare with that of the exact solution.The coarse and fine mesh results are converging toward exact from above.For theQ5Q4case the results are visibly the same as exact and the color bar values are likewise the same as exact.

The coarse mesh results are enormously in error,both qualitatively and with respect to amplitude.ForQ3Q2,the amplitude error is approximately one order of magnitude;the displacements are simply much too large,the antithesis of locking.

What is the explanation for the results of this problem? One might observe that the specification of the problem is somewhat contradictory.The material characterization is nearlyincompressible,indicated by the ratio of the Lamé parameters,λ/μ=104.However,the manufactured solution does not respect this characterization in that its divergence is of the same order as its deviatoric component.The upshot is that the applied forces have a pressure component four orders of magnitude greater than the forces proportional toμ.It is clear from the large displacement response that those forces are being balanced partially by deviatoric deformations.

In addition to the unphysical nature of this problem we may point out that the results of the Abaqusq1p0element are much the worst in error,both for displacements and pressure,and this element is widely used and considered appropriate for broad classes of real-world applications.However,we take the results of this problem to be cautionery,and indicative of the fact that both Galerkin and collocation results for mathematically well-posed problems may prove wanting in certain extreme cases.As a rule of thumb,it is clear that higher-order collocation results are significantly better than lower order,but even for higher order,sufficient mesh refinement of results is still necessary to obtain very accurate results.

We would now like to describe what we have seen from a mathematical point of view,termed“pressure robustness,” and thoroughly developed by Alexander Linke and his collaborators;see,e.g.,(18) and references therein.Pressure robustness is an additional attribute to inf-sup stability asked of a discrete scheme by Linke.It can be stated for our problem as follows:part of the force scaled byλis increased,holdingμfixed,the numerical displacement solution should beO(1),that is independent ofλ.Our results clearly indicate this is not the case,and so we conclude none of the methods considered is pressure robust.In all cases they areO(λ) but perhaps mollified by some power of the mesh parameterh.Over the small range of meshes considered for this problem,10×10×10 to 35×35×35,we might expect roughly three orders of magnitude deterioration of the solution withλand this is the case.In the array of standard primal and mixed finite elements used in structural mechanics,the combination of inf-sup stability and pressure robustness is elusive.We note,however,that stabilized methods seem to alleviate this issue when most others fail,but their use in structural mechanics has been limited despite widespread use in fluids.

5 Conclusions and Summary of Results

In this paper we have investigated the behavior of isogeometric collocation methods in linear,isotropic elasticity with emphasis on the behavior in the nearly-incompressible case,taken asλ/μ=104herein,whereλandμare the Lamé parameters.We also presented results for the compressible case,taken asλ/μ=1 herein,for comparison.Two formulations of the problem were considered;the primal formulation utilizing the Navier equations of elasticity,and the mixed,displacement-pressure formulation.For the primal formulation we investigated NURBS based models of polynomial degreep=2,3,...,7.For the mixed elements,we considered displacement fields of degreep=3,4,...,7 and pressures one degree lower,a combination that has been shown to be effective previously in isogeometric Galerkin methods [2,3].The primal formulation,when implemented in a standard,low-degree,Galerkin finite element method,is well known to suffer from displacement “locking.” The mixed formulation has been used to design both low-and high-degree finite elements that alleviate locking and have been successfully applied to engineering problems.The three numerical examples we considered herein asked similar questions of the isogeometric collocation methods.For context,we presented results for constant pressure Abaqus elements that are known not to lock and are widely used in engineering analysis.

A general statement about the results is that for the compressible case,there were no surprises.Both primal and mixed cases presented consistent convergence patterns,with higher-degree basis functions exhibiting greater accuracy.

For all problems studied,consistent patterns also emerged for the primal formulation.Locking was observed whenλ/μ=104forp=2 and 3.Higher-degree cases behaved better but the absolute values of displacement errors for the nearly-incompressible case were typically orders of magnitude higher than for the compressible case,despite convergence rates being the same.We conclude from these results that isogeometric collocation with the primal formulation is not effective for nearly-incompressible applications.

For the mixed formulation,applied to the two-dimensional problems,very good results were obtained in all cases.Even for the lowest degree case,Q3Q2,displacement errors were commensurate with the Abaqusq1p0element in magnitude and exactly the same in rate of convergence,namely,O(h2);for pressures,the mixed collocation results were more accurate in magnitude and in rate,and againO(h2).The higher-degree cases were much more accurate.Convergence rates for both displacement and pressure were the same and ofO(hp) andO(hp-1),for even and odd numbered elements,respectively,and there was very little difference in error amplitudes for the nearly incompressible and compressible cases.

The three-dimensional problem presented a more complex story for the mixed method.The loading in this case had bothλ-andμ-proportional terms,and theλterm dominated in the nearly-incompressible case.For theλ/μ=1 case,all results behaved as expected.The pressure results did not degrade substantially for theλ/μ=104case compared with theλ/μ=1 case.However,the displacement results were roughly three orders of magnitude worse in magnitude,although the rates of convergence were the same.It appears that the “constant” in the displacement error is orderλ,perhaps mollified by a power ofh.One might have initially assumed that this was due to locking,but scrutiny of the results indicated it was quite the opposite;the displacements were much too large.The explanation is as follows:Theλ-proportional body force term is precisely a pressure gradient,but the discrete pressure-gradient term in the mixed methods failed to fully balance it.Part of the burden was shifted to theμ-term,thus producing displacements of orderλ/μ.

The poor performance of mixed methods of the type studied here,when subjected to body forces that are gradients of a scalar,has been described as a lack of “pressure robustness” by Alexander Linke and collaborators [18].Unfortunately,lack of pressure robustness is the rule rather than the exception and it seems to be at odds with most “inf-sup” stable methods.It has been suggested that “stabilized methods,” which accommodate any combination of displacement and pressure basis functions,in particular equal-order,might be able to alleviate stability and pressure robustness issues simultaneously,and we hope to investigate this in future work;see,e.g.,Hughes et al.[26,27].

For all the problems studied,we presented corresponding results using Abaqus withq1p0elements.These elements are widely used in practice and provide a benchmark for what may be considered acceptable in engineering applications.The results for the lowest-degree mixed collocation method,Q3Q2,were commensurate in accuracy with theq1p0element for displacements,but considerably more accurate for pressures.The higher-degree mixed collocation methods were substantially more accurate for both,as might have been anticipated.We do not suggest this is an equitable comparison,but rather it serves as an indication that the mixed collocation methods achieve a level of accuracy that makes them viable for engineering computations.

Funding Statement:FF and LDL gratefully acknowledge the financial support of the German Research Foundation (DFG) within the DFG Priority Program SPP 1748 “Reliable Simulation Techniques in Solid Mechanics”.AR has been partially supported by the MIUR-PRIN project XFAST-SIMS (No.20173C478 N).

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.