A Virtual Element Method for transversely isotropic elasticity
AA Virtual Element Method for transversely isotropic elasticity
D. van Huyssteen B. D. Reddy Abstract
This work studies the approximation of plane problems concerning transverselyisotropic elasticity, using a low-order Virtual Element Method (VEM), with a focus on near-incompressibility and near-inextensibility. Additionally, both homogeneous problems, in whichthe plane of isotropy is fixed; and non-homogeneous problems, in which the fibre direction defin-ing the isotropy plane varies with position, are explored. In the latter case various options areconsidered for approximating the non-homogeneous fibre directions at element level. Througha range of numerical examples the VEM approximations are shown to be robust and locking-free for several element geometries and for fibre directions that correspond to mild and strongnon-homogeneity.
Keywords
Virtual Element Method · Transversely isotropic · Linear elasticity
The popular finite element method has the status of a classical approach for obtaining approxi-mate solutions to problems formulated as systems of partial differential equations or inequalities,or alternatively, in their variational form. Particularly in the domains of solid and fluid mechan-ics, the method has been used with great success for problems with high degrees of complexitysuch as non-linear problems and problems with intricate geometries (see for example the treat-ments in [1, 2]).A number of variants of standard conforming finite element methods have been developed overthe last four decades, with a range of motivations in mind. Mixed methods, for example, haveallowed all variables of interest to be approximated explicitly; and in addition have providedavenues through which stable and convergent finite element approximations can be developed insituations where the selection of the values of certain parameters might lead to non-convergence.Key examples are those of near-incompressibility, or problems in structural mechanics in whichthe geometry is characterized by a small length scale. These two features lead, in the contextof low-order standard finite element methods, to volumetric and shear locking respectively.Phenomena that may be circumvented by the use of mixed methods [3, 4]. Department of Mechanical Engineering and Centre for Research in Computational and Applied Mechanics,University of Cape Town, 7701 Rondebosch, South Africa. Email: [email protected] Department of Mathematics and Applied Mathematics and Centre for Research in Computational and AppliedMechanics, University of Cape Town, 7701 Rondebosch, South Africa. Email: [email protected] a r X i v : . [ m a t h . NA ] M a r et another variant of the standard conforming finite element method is the discontinuousGalerkin (DG) method, in which interelement continuity is abandoned in favour of greaterflexibility with regard to meshing (see for example [5]). In addition, DG methods, when designedappropriately, are stable and uniformly convergent in situations of near-incompressibility forlow-order approximations [6, 7, 8].A more recent development in the context of finite element methods is the Virtual ElementMethod (VEM). In contrast to the geometric restrictions on finite elements, which are mostgenerally triangular or quadrilateral in 2D, and tetrahedral or hexahedral in 3D, the VEMpermits elements to be arbitrary polygons in 2D or polyhedra 3D. Furthermore, there is no needfor elements to be convex, and degeneracies such as element sides having small interior anglesor arbitrarily small edges pose no problems. Some key works in a rapidly growing literatureinclude [9, 10, 11]. Applications of the VEM to nonlinear problems include works on nonlinearelasticity [12, 13], elastoplasticity [14, 15, 16], and contact [17].Applications of the VEM to elasticity have been largely confined to the isotropic problem, al-though there have been treatments of inextensible materials [18]. Problems involving anisotropicmaterials pose additional challenges in the context of VEM approaches, particularly for non-homogeneous materials in which the anisotropy varies with position. In [19] limiting extensibilityis investigated in an otherwise isotropic material using penalty, Lagrange multiplier, and per-turbed Lagrangian approaches. The work was considerably extended by Rasalofson et. al.[20]. This work presents a detailed treatment of the boundary value problem for transverselyisotropic linear elastic materials. Conditions for well-posedness are established, and finite ele-ment approximations are studied using both conforming and reduced integration approaches.An error analysis gives an indication of conditions under which low-order approximations areuniformly convergent in the incompressible and inextensible limits, and a set of numerical exper-iments provides further insight into the conditions under which locking-free behaviour occurs.Specifically, with the degree of anisotropy measured through the ratio of Young’s modulus in thefibre direction relative to that in the plane of isotropy, it is shown in this work that locking-freebehaviour occurs in conditions of moderate anisotropy for low-order conforming quadrilaterals,in contrast to the situation for isotropic materials. Furthermore, for high degrees of anisotropyleading to near-inextensibility, locking occurs, but is circumvented by the use of selective under-integration.The purpose of this work is to study low-order VEM approximations for plane problems con-cerning transversely isotropic elasticity. Of particular interest is the behaviour of VEM approx-imations for the limiting situations of near-incompressibility and near-inextensibility. Whereasin the case of conventional finite element approximations, as discussed above, some form ofmodification such as selective under-integration is necessary in order to circumvent locking, inthe case of VEM approximations locking-free behaviour is observed in the incompressible andinextensible limits. 2 further novel aspect of this work is the treatment of non-homogeneous transverse isotropy;that is, situations in which fibre directions vary with position. Here it becomes necessary toapproximate the non-homogeneous terms appropriately in order to preserve the simplicity ofthe VEM approach, in which integrals are evaluated only on element boundaries. The approx-imations adopted are shown to be robust, with the locking-free behaviour also evident for thenon-homogeneous problem.The structure of the rest of this work is as follows. Section 2 sets out the details of the constitutiverelations for transversely isotropic linear elastic materials, the set of governing equations, andthe associated weak formulation. The details of the Virtual Element Method are presented inSection 3, and the set of numerical results are presented and discussed in Section 4. This workconcludes with a summary of results and a discussion of open problems. Consider a linear elastic body which occupies a plane, polygonal bounded domain Ω ⊂ R withboundary ∂ Ω. The boundary comprises a non-trivial Dirichlet part Γ D and Neumann part Γ N such that Γ D ∩ Γ N = ∅ and Γ D ∪ Γ N = ∂ Ω. Transversely isotropic materials exhibit isotropic behaviour in a specified plane, this plane beingdefined by a normal vector, and referred to also as the fibre direction.The Cauchy stress tensor σ is related to the infinitesimal strain tensor ε through the elasticrelation σ = C ε . (2.1)Here C is a fourth-order tensor of elastic moduli. For a transversely isotropic material with thedirection of transverse isotropy defined by the unit vector a , (2.1) takes the form [20] σ = λ (tr ε ) I + 2 µ T ε + β ( M : ε ) M + α (( M : ε ) I + (tr ε ) M ) + 2( µ L − µ T )( εM + M ε ) . (2.2)Here M = a ⊗ a , λ and µ T are the conventional Lam´e parameters, µ L is the shear modulusin the fibre direction, and I denotes the second-order identity tensor. The material constants α and β do not have a direct interpretation, though it will be seen that β → ∞ in the limit ofinextensible behaviour in the fibre direction.The special case of an isotropic material is recovered by setting α = β = 0 and µ L = µ T .The five material constants in (2.2) may be related to the “engineering” constants, viz. Young’smoduli E L and E T in the fibre direction and plane of isotropy, respectively, and the corresponding3oisson’s ratios ν L and ν T , by inverting (2.2), specializing it to the case in which a = e , andcomparing with the compliance relation written in the form (see for example [21]) ε ε ε (cid:15) (cid:15) (cid:15) = E T − ν T E T − ν L E L − ν T E T E T − ν L E L − ν L E L − ν L E L E L µ L µ L
00 0 0 0 0 1 µ T σ σ σ σ σ σ . (2.3)In the remainder of this work we make the assumption, with little loss in generality, that ν T = ν L := ν and µ T = µ L := µ . Further, we set p = E L E T , (2.4)so that the parameter p measures the degree of transverse isotropy, with inextensible behaviourcorresponding to the limit p → ∞ .The parameters in (2.2) may then be expressed in terms of the engineering parameters as [20] λE T = ν ( ν + p ) D ,αE T = ν ( p − D , (2.5) βE T = p (1 − ν ) − p (1 + 2 ν ) + 3 ν D , in which the denominator D is given by D = (1 + ν )( p (1 − ν ) − ν ) . (2.6)We also have the relation µ T = E T ν ) . (2.7)Furthermore, noting that D → ν → when p →
1, which corre-4ponds to the isotropic limit, it is evident from (2.5) and (2.6) that (cid:40) λ is bounded as ν → , if p >
1, and as p → ∞ (inextensibility) λ → ∞ as ν → , for p = 1 (isotropy) (2.8) (cid:40) α is bounded as ν → , if p > α → p → (cid:40) β is bounded as ν → , if p > β → ∞ as p → ∞ (inextensibility) . (2.10)The elasticity tensor is assumed to be pointwise stable; that is, to satisfy the condition ε : C ε > ε . General conditions on the material constants for pointwise stability are somewhat complex [20],but the simple set of conditions λ + 23 µ > , µ > , p ≥ The body is subjected to a body force f , prescribed loading h on Γ N , and a prescribed displace-ment g on Γ D .The equation of equilibrium is − div σ = f on Ω . (2.12)Small displacements are assumed, and the strain displacement relation is ε ( u ) = ( ∇ u + [ ∇ u ] T ) or ε ij ( u ) = ( u i,j + u j,i ) . (2.13)Here u denotes the displacement, and ∇ u the displacement gradient with components u i,j . Hereand henceforth we choose a fixed Cartesian coordinate system x i with orthonormal basis e i .The boundary conditions are u = g on Γ D , (2.14a) σ · n = h on Γ N . (2.14b)Equations (2.12) – (2.14), together with the elastic relation (2.2), constitute the boundary-valueproblem for a transversely isotropic body. 5 .3 Weak formulation We denote by L (Ω) the space of square-integrable functions on Ω, and by H (Ω) the Sobolevspace of functions which, together with their generalized first derivatives, are square-integrable,and set V = [ H D (Ω)] d = { v | v i ∈ H (Ω) , v = on Γ D } .We also define the function u g ∈ [ H (Ω)] d such that u g | Γ D = g . The bilinear form a ( · , · ) and linear functional (cid:96) ( · ) are defined by a : [ H (Ω)] d × [ H (Ω)] d → R , a ( u , v ) = (cid:90) Ω σ ( u ) : ε ( v ) dx, (2.15a) l : [ H (Ω)] d → R , (cid:96) ( v ) = (cid:90) Ω f · v dx + (cid:90) Γ N h · v ds − a ( u g , v ) . (2.15b)The weak form of the problem is then as follows: given f ∈ [ L (Ω)] d and h ∈ [ L (Γ N )] d , find U ∈ [ H (Ω)] d such that U = u + u g , u ∈ V , and a ( u , v ) = (cid:96) ( v ) ∀ v ∈ V. (2.16)We write the bilinear form as a ( u , v ) = a iso ( u , v ) + a ti ( u , v ) , (2.17)where a iso ( u , v ) = λ (cid:90) Ω ( ∇ · u )( ∇ · v ) dx + 2 µ (cid:90) Ω ε ( u ) : ε ( v ) dx, (2.18a) a ti ( u , v ) = α (cid:90) Ω [( M : ε ( u ))( ∇ · v ) + ( ∇ · u )( M : ε ( v ))] dx + β (cid:90) Ω ( M : ε ( u ))( M : ε ( v )) dx . (2.18b)The bilinear form is clearly symmetric. The well-posedness of the weak problem requires thebilinear form to be continuous and coercive, and the linear functional continuous. With theassumptions (2.11), it is shown in [20] that the problem has a unique solution that dependscontinuously on the data. The domain Ω is partitioned into a mesh of elements comprising non-overlapping polygons E with ∪ E = ¯Ω. A typical polygonal element is shown in Figure 1. A generic edge in a polygonalelement is denoted by e and a vertex by i , so that e = 1 , . . . , N and i = 1 , . . . , N , where N isthe number of vertices of an element. 6 igure 1: An arbitrary virtual element We construct a conforming approximation in a space V h ⊂ V . The space V h comprises functionsthat are continuous on Ω, piecewise linear on the boundary ∂E of each element, and withdiv C ε ( v h ) vanishing on E [9, 10]: V h = { v h ∈ V | v h ∈ [ C (Ω)] , div C ε ( v h ) = on E, v h | e ∈ P ( e ) } . (3.1)Here and henceforth P k ( X ) denotes the space of polynomials of degree ≤ k on the set X ⊂ R d ( d = 1 , v h | E = ϕd (3.2)in which ϕ denotes a matrix of virtual basis functions, and d is the 2 N × e of elements, and it is convenient to writealso v h | ∂E = N d and ε ( v h ) = Bd , (3.3)in which N and B are respectively matrices of standard Lagrangian linear basis functions andtheir derivatives. Thus, the basis functions ϕ are not known, and are not required to be known;their traces on the boundary are however required, and are simple Lagrangian functions.We will require the projection Π : V h | E → P ( E ), defined on E by (cid:90) E Π v h dx = (cid:90) E ε ( v h ) dx . (3.4)Thus Π is the L -orthogonal projection onto constants of the strain associated with the displace-ment v h on an element E . From (3.3), and given that Π v h is constant we have, in componentform, (Π v h ) ij = 12 1 | E | (cid:90) E (( v h ) i,j + ( v h ) j,i ) dx = 12 1 | E | (cid:73) ∂E (( v h ) i n j + ( v h ) j n i ) ds = 12 1 | E | (cid:88) e ∈ ∂E (cid:90) e ( N iA d EA n j + N jA d EA n i ) ds . (3.5)7ere d EA denotes the degrees of freedom associated with element E , summation is implied overall repeated indices, and we have used integration by parts and the representation (3.3) . Theintegrals in (3.5) are readily evaluated as the edge basis functions are known. Thus the projectionΠ v h is available as a function of the degrees of freedom.To construct the virtual element formulation we start by writing a E ( u , v ) := a ( v , v ) | E = (cid:90) E ε ( u h ) : C ε ( v h ) dx , (3.6)so that a E ( · , · ) denotes the contribution of element E to the bilinear form a ( · , · ) defined in (2.16)and (3.7a). We have a E ( u h , v h ) = (cid:90) E Π u h : C Π v h dx + (cid:90) E ( ε ( u h ) − Π u h ) : C ( ε ( v h ) − Π v h ) dx (3.7a)+ (cid:90) E Π u h : C ( ε ( v h ) − Π v h ) dx + (cid:90) E ( ε ( u h ) − Π u h ) : C Π v h dx (3.7b)= (cid:90) E Π u h : C Π v h dx + (cid:90) E ( ε ( u h ) − Π u h ) : C ( ε ( v h ) − Π v h ) dx (3.7c)= (cid:90) E Π u h : C Π v h dx (cid:124) (cid:123)(cid:122) (cid:125) consistency term + (cid:90) E (cid:2) ε ( u h ) : C ε ( v h ) − Π u h : C Π v h (cid:3) dx (cid:124) (cid:123)(cid:122) (cid:125) stabilisation term (3.7d)The last line is obtained by noting the definition of the projection operator, so that the twoterms in (3.7b) are zero. Furthermore, the definition of the projection is invoked again in goingfrom (3.7c) to (3.7d). The terms in (3.7d) are referred to respectively as the consistency termand stabilisation term. The consistency term.
After substitution of (3.5) in the consistency term, evaluation ofthe integral leads to the expression (cid:90) E Π u h : C Π v h dx = (¯ d ) T K E con d (3.8)in which K E con is the consistency stiffness matrix for element E and d E and ¯ d E are respectivelythe vectors of nodal degrees of freedom of u h and v h on element E . The stabilisation term.
Use of the consistency term alone would lead to a rank-deficientstiffness matrix. The second term on the right hand side of (3.7d) serves the purpose of stabilizingthe formulation. The basic idea behind the VEM is that integrals are evaluated on the boundariesof elements only; the stabilisation term in its original form would require that integrals beevaluated on the elements. Nevertheless, it is not necessary to evaluate this term exactly, and itsuffices to replace it with an approximation. There are several methods that can be employed,8ee for example [15, 11]. However, we choose the stabilisation method presented in [22] as it hasproven very robust, a E stab ( u h , v h ) = τ d T [ I − D ( D T D ) − D T ] d . (3.9)Here d and d are, again, respectively the vectors of nodal degrees of freedom associated with v h and u h , and D is the matrix that relates the nodal degrees of freedom d of a linear vectorpolynomial to its degrees of freedom s relative to a scaled linear monomial basis. That is, foran element with N nodes, d = Ds . (3.10)Note that D has dimensions 2 N ×
6, and has the basis monomials M = { , ξ, η } = (cid:26) , x − x c d E , y − y c d E (cid:27) , (3.11)where d E is the diameter of element E , with x c and y c the x − and y − coordinates of the centroidof E respectively.This approximation may be motivated by seeking a stabilisation term of the form τ ( d T d − d T d ) (3.12)in which d are the nodal degrees of freedom of a linear polynomial that is closest to u h in somesense, and τ is a suitable scaler to be chosen. In the event that u h is a linear polynomial, thenthis term vanishes of course.From (3.10) we have d T d = ( s T D T )( Ds )= s T ( D T D )( D T D ) − ( D T D ) s = d T D ( D T D ) − D T d . (3.13)Then we obtain (3.9) by replacing d with the actual vector degrees of freedom.We need to choose a suitable value for the scalar τ , such that it is some value representative ofthe constitutive tensor. We consider the transversely isotropic material properties λ , α , β , µ L and µ T . As seen in Section 2.1 λ, β → ∞ as ν → .
5, to keep the VEM locking free we thereforereject these options. We choose τ = µ T as it is bounded and is representative of both isotropicand transversely isotropic materials. As we have set µ := µ L = µ T , we then have K E stab = µ (cid:104) I − D (cid:0) D T D (cid:1) − D T (cid:105) . (3.14)As we have used scaled coordinates, no area scaling of the stabilisation term is necessary. Thecomplete stiffness matrix is then given by K E = K E con + K E stab . (3.15)9 Numerical results
In this section we present numerical results for three model problems to illustrate the perfor-mance of the VEM. We consider homogeneous materials, for which the plane of isotropy is fixedacross the domain, and also non-homogenous materials, for which the plane of isotropy, as de-fined by the vector a , varies with position. Plane strain conditions are assumed. As in Section2.1 we set ν T = ν L = ν and µ T = µ L = µ . We consider values of p > ν = 0 . ν = 0 . a := (cid:92) ( Ox, a ) to be the angle between the x -axis and the fibre direction a . The resultsin the examples that follow are obtained for the following element types:Q The standard bilinear quadrilateralQ The standard biquadratic quadrilateralQuad The VEM formulation with four-noded elementsHex The VEM formulation with six-noded elementsVoronoi The VEM formulation with Voronoi elementsFigure 2 depicts patches of the meshes comprising six-noded and Voronoi elements for a meshdensity d of 7, where d = √ n elements . Meshes are constructed on a unit domain and then mappedto the problem domain, Figure 2 depicts meshes after this mapping. (a) Hex mesh (b) Voronoi meshFigure 2: Cook’s membrane problem, showing the hexagonal and Voronoi meshes for mesh density 7 We present results here for the case in which fibre directions are constant on the domain. Theemphasis is on near-incompressibility and near-inextensibility, either separately or combined. Inall examples Poisson’s ratio is set at ν = 0 . p > ook’s membrane problem. This problem consists of a trapezoidal panel fully fixed alongits left edge with a uniformly distributed load along its right edge, as shown in Figure 3. Theapplied load is P = 100N and E T = 250Pa. This test problem has no analytical solution. Thevertical displacement at point C is recorded. Figure 3: Cook’s membrane problem, showing fibres inclined at ˆ a = π Figure 4 shows a convergence plot of tip displacement vs mesh density for fibre angle ˆ a = π ,as illustrated in Figure 3, and with p = 5, for the VEM formulation with the three candidatemeshes, and for standard finite element formulations. The various VEM formulations are seento exhibit degrees of accuracy comparable to that of the Q approximation. . .
54 Mesh Density D i s p l a c e m e n t QuadHexVoronoi Q Q Figure 4: The Cook problem: convergence test for fibre angle ˆ a = π and p = 5 The results that follow have been generated using meshes with a density of d = 50.11igure 5 shows semilog plots of tip displacement vs p for 1 ≤ p ≤ and for fibre angles ˆ a = π and ˆ a = π . The well-known locking behaviour of Q is clear in the isotropic limit ( p → Q , and is locking-free. Quad, HexVoronoi, Q Q p D i s p l a c e m e n t QuadHexVoronoi Q Q Quad, HexVoronoi, Q Q p D i s p l a c e m e n t QuadHexVoronoi Q Q Figure 5: The Cook problem: tip displacement vs p for (a) fibre direction ˆ a = π ; (b) fibre directionˆ a = π Figure 6 shows a plot of tip displacement vs fibre orientation for a nearly inextensible material( p = 10 ). Again, we note the poor performance and locking behaviour of Q over most of therange, and on the other hand the robust behaviour of the VEM formulation. The Q elementdisplays sub-optimal accuracy for fibre angles greater than ˆ a = π and close to zero. This issomewhat surprising, in that the behaviour of this element in the near-inextensible limit wouldbe expected to mirror its good performance for near-incompressibility. On the other hand, whilethe element has been shown to be uniformly convergent for incompressible materials, there doesnot exist a corresponding analysis for near-inextensibility, to the best of the authors’ knowledge.Such an analysis could shed light on the behaviour seen in Figure 6.12
50 100 15002468 Q ˆ a (Degrees) D i s p l a c e m e n t QuadHexVoronoi Q Q Figure 6: The Cook problem: tip displacement vs fibre orientation, for p = 10 The beam problem.
This problem consists of a beam subject to a linearly varying load at itsright edge, and pinned at its left extrema, as depicted in Figure 7. The load has maximum andminimum values of F max = ± w = 10m, height h = 2m and Young’sModulus of E T = 1500Pa. The vertical displacement at point C is recorded.The displacement of point C is given by [20] u ( x, y ) = 2 F max h (cid:20) S xy + S (cid:18) y − h (cid:19)(cid:21) , (4.1) v ( x, y ) = F max h (cid:20) S (cid:18) y − h (cid:19) − S x (cid:21) ; (4.2)the coefficients S ij are lengthy functions of the material constants, and are given in the Appendixto [20]. Figure 7: The beam problem, showing fibres inclined at ˆ a = π Figure 8 shows a convergence plot of tip displacement vs mesh density for a fibre orientation ofˆ a = π , and with p = 5. It is seen that for the various VEM meshes the convergence behaviouris similar to that of the Q mesh for sufficiently fine meshes.13
10 20 30 40 500 . . . . D i s p l a c e m e n t QuadHexVoronoi Q Q Sol
Figure 8: The beam problem: convergence test for fibre angle ˆ a = π and p = 5 The results that follow have been generated using meshes with a density of d = 50.Figure 9 shows semilog plots of tip displacement vs p for 1 ≤ p ≤ for fibre angles of ˆ a = π and ˆ a = π . Again, as with the Cook problem, the VEM solutions are locking-free and displayhigh accuracy. As pointed out in [20], for mild anisotropy, that is, low values of p , the tendencyto lock for the Q mesh is mitigated as a result of the Lam´e parameter being bounded for p > ν very close to 0.5. This behaviour is evident inFigure 9, where the Q mesh is seen to be locking-free for p > p up to p ≈
10 for fibreangle ˆ a = π , and p ≈
100 for fibre angle ˆ a = π . . . . . QuadHexVoronoi Q Sol Q p D i s p l a c e m e n t QuadHexVoronoi Q Q Sol . . . . Quad, HexVoronoi, Q , Sol Q p D i s p l a c e m e n t QuadHexVoronoi Q Q Sol
Figure 9: Beam problem: tip displacement vs p for fibre angle (a) ˆ a = π ; (b) ˆ a = π Figure 10 shows a plot of tip displacement vs fibre orientation for the case of near-inextensibility( p = 10 ). Again, we note poor performance of Q and robust and accurate behaviour of14he VEM formulations. In contrast to the result for the Cook problem, here the Q elementdemonstrates equally accurate behaviour. . . . . QuadHex Voronoi Q , Sol Q ˆ a (Degrees) D i s p l a c e m e n t QuadHexVoronoi Q Q Sol
Figure 10: The beam problem: tip displacement vs fibre orientation, for p = 10 For a given distribution of fibre directions a ( x ) it follows that some approximation has to bemade for a ( x ) within each element so as to preserve the general approach to carrying out VEMcomputations. A simple option would be to approximate a by its centroidal value. However,such an approximation can be somewhat inaccurate for situations in which the fibre orientationvaries significantly across a length scale comparable to mesh size. A more reliable approach isto use the average fibre direction at the element nodes; this approach is observed to yield morestable and faster convergence. When dealing with rapidly varying fibre directions a more stableapproach was achieved by using a weighted average of the fibre direction at the centroid andthe average direction at the vertices. This approach applies an equal weighting to the centroidaldirection and the average of its nodal values for very coarse meshes; for finer meshes, that is,as the mesh density increases, the weighting of the centroidal value decreases rapidly. Thus forvery fine meshes and rapidly varying fibre directions, it is the nodal average that dominates.With the mesh density denoted by d , the centroidal weight w c is defined by w c = π + arctan( d cr − d )2 π . (4.3)This function is shown in Figure 11. Here d cr is a user-defined critical mesh density beyondwhich the value of the weight drops rapidly.The average fibre direction a is then given by a ave | E = w c a ( x c ) + (1 − w c ) 1 N N (cid:88) i =1 a ( x i ) , (4.4)15 . . . . . d cr Mesh Density w c Figure 11: The weight w c as a function of mesh density d where, as before, N denotes the number of nodes of element E , and x c and x i are respectivelythe coordinates of the element centroid and node i . We then approximate the elasticity tensoron an element E by C | E (cid:39) C ( a ave ) . (4.5)The critical density used is problem specific as it depends on the degree of variation in fibreorientation. However, for simplicity a critical density of d crit = 10 was used as it worked wellacross a range of problems.Except where otherwise stated, Poisson’s ratio is set at ν = 0 . y = c + f ( x )where f ( x ) is chosen to be, respectively, ( x − ( x − x −
36) and 2 sin x . The polynomialdistribution corresponds to mild variation with position, while the sinusoidal distribution is amore severe test of performance under conditions of rapidly varying direction. Figure 12 showsschematically the curves corresponding to these two cases for the Cook problem, one of theexamples considered in what follows. Figure 12: Cook’s membrane problem with curves showing variable fibre orientation for (a) quartic,and (b) sinusoidal distributions ook’s membrane problem. Figures 13 (a) and (b) show the tip displacement as a functionof mesh density for fibres corresponding to the quartic distribution, with p = 5, and in whichthe value of a is based respectively on its value at the element centroids and the average of itsvalues at the nodes. We note smooth and stable convergence of the VEM, with the quadrilateralVEM mesh performing somewhat more poorly for the case in which the average nodal value of a is used. D i s p l a c e m e n t QuadHexVoronoi Q Q (a) D i s p l a c e m e n t QuadHexVoronoi Q Q (b)Figure 13: Tip displacement vs mesh density for Cook’s membrane problem for p = 5, with fibredirections defined by quartic curves, and using (a) the centroidal value of a ; (b) the average of nodalvalues of a Figure 14 shows tip displacement as a function of mesh density for fibres corresponding to thesinusoidal distribution. In Figure 14 (a) we calculate a ave based on its value at the elementcentroid. Again we see good performance by the VEM elements, with an accuracy comparableto that of Q . In contrast to the results for polynomial variation in Figures 13 (a) and (b), thecoarse-mesh behaviour is somewhat erratic, with a smooth dependence on mesh density onlyafter d ≈
25. Figures 14 (c) and (d) present results for the cases in which, respectively, firstly anequal weighting of centroidal and nodal values of a is used, and secondly, using (4.4), a varyingweight is used. Similar behaviour is seen when compared with the results in Figure 13 (b),though the dependence on mesh density becomes smoother for coarser meshes at d ≈
10 20 30 40 50234567 Mesh Density D i s p l a c e m e n t Quad HexVoronoi Q Q (a) D i s p l a c e m e n t Quad HexVoronoi Q Q (b) D i s p l a c e m e n t Quad HexVoronoi Q Q (c) D i s p l a c e m e n t Quad HexVoronoi Q Q (d)Figure 14: Tip displacement vs mesh density for Cook’s membrane problem and with fibre directionsdefined by curves 2 sin x , and using (a) the centroidal value of a ; (b) the average of nodal values of a ; (c) equal weighting of nodal and centroidal values; and (d) a varying weighted average as in (4.3)and (4.4) Next, we consider behaviour in the near-incompressible limit, with ν = 0 . p , for the quartic and sinusoidal fibre distributionsrespectively. For the polynomial fibre distribution the centroidal values of fibre direction areused, while the weighted method is used for the sinusoidal distribution. There is little variationin the performance of the various VEM meshes, though for the Voronoi mesh and for near-inextensibility small scatter is observed. The sub-optimal behaviour of the Q mesh seen in theCook example in Figure 6 is not evident here. The Q mesh again displays locking behaviourexcept in a narrow range of mild anisotropy. 18 p D i s p l a c e m e n t Quad HexVoronoi Q Q (a) p D i s p l a c e m e n t Quad HexVoronoi Q Q (b)Figure 15: Tip displacement vs p for the Cook problem, for near-incompressibility and using (a) thequartic, and (b) the sinusoidal fibre distributions Beam in bending.
We consider next the problem of a beam in bending, shown in Figure16, with boundary conditions slightly different to those shown in Figure 7; the left edge is nowconstrained horizontally and pinned at the bottom left corner. The fibre distributions consideredare once again quartic and sinusoidal, as for the Cook problem, and are depicted in Figure 16.However, the quartic distribution is now defined by y = ( x − ( x − . x − .
5) + c , and thesinusoidal distribution is as defined previously. The vertical displacement at point C is recorded.A value of Poisson’s ratio of ν = 0 . Figure 16: Quartic and sinusoidal fibre distributions for the beam in bending problem
In Figures 17 (a) and (b) we present convergence plots of vertical displacement vs mesh density,for mild anisotropy; that is, p = 5, for the quartic distribution of fibres, and with a ave calculatedusing respectively the centroidal value, and the average nodal value. We note stable convergencein all cases. For the case in which the average nodal direction is used, it is seen that the standard Q element performs best for coarse meshes. 19
10 20 30 40 502345 Mesh Density D i s p l a c e m e n t QuadHexVoronoi Q Q D i s p l a c e m e n t QuadHexVoronoi Q Q Figure 17: Tip displacement vs mesh density for the beam in bending problem, with p = 5, andfibre directions defined by a quartic polynomial, using (a) the average of nodal values of a ; (b) thecentroidal value of a Figure 18 shows tip displacement as a function of mesh density for a sinusoidal distribution offibres. In Figure 18 (a) we calculate a ave based on its centroidal value. Again we see good perfor-mance by the VEM elements, though the Q mesh performs best. Rapid numerical convergenceis however observed after a density of d ≈
25. 20
10 20 30 40 5012345 Mesh Density D i s p l a c e m e n t Quad HexVoronoi Q Q (a) Fibre Direction At Centroid D i s p l a c e m e n t Quad HexVoronoi Q Q (b) Average Fibre Direction At Vertices D i s p l a c e m e n t Quad HexVoronoi Q Q (c) Constant Weighted Average - w = D i s p l a c e m e n t Quad HexVoronoi Q Q (d) Varying Weighted AverageFigure 18: Tip displacement vs mesh density for the beam in bending problem, with fibre directionsdefined by curves 2 sin x , and using (a) the centroidal value of a ; (b) the average of nodal values of a ; (c) equal weighting of nodal and centroidal values; and (d) a varying weighted average as in (4.3)and (4.4) As with the Cook problem, we next consider behaviour in the near-incompressible limit with ν = 0 . p , for the quartic and sinusoidalfibre distributions. For the polynomial fibre distribution the average value of fibre direction atthe vertices is used, while the weighted method is used for the sinusoidal distribution. Thereis little variation in the performance of the various VEM meshes, though for the Voronoi meshand for near-inextensibility small scatter is again observed. The Q mesh performs ratherpoorly, displaying some evidence of mild locking. This should be compared with the sub-optimal behaviour seen in Figure 6 for constant fibre directions. Except in a narrow range of21ild anisotropy, the Q mesh displays locking behaviour. p D i s p l a c e m e n t QuadHexVoronoi Q Q (a) p D i s p l a c e m e n t QuadHexVoronoi Q Q (b)Figure 19: Tip displacement vs p for the beam in bending problem, for near-incompressibility andusing (a) the quartic, and (b) the sinusoidal fibre distributions In this work we have formulated and implemented a Virtual Element Method for plane trans-versely isotropic elasticity, making provision for homogeneous as well as non-homogeneous bod-ies. In the latter case, various options for taking account of the non-constant elasticity tensorare investigated. The formulations have been studied numerically through two model problems,and for three different kinds of polygonal meshes. The results have been compared againstthose obtained using conventional conforming finite element approximations with bilinear andbiquadratic approximations.The VEM approximations are found to be locking-free for both near-incompressibility and near-inextensibility, without the need to make modifications to the formulation. In the case of finiteelement approximations, the well-known volumetric locking behaviour of bilinear approximationsis evident, except for a range of parameters corresponding to mild anisotropy. This behaviouris consistent with what has been shown in [20]; for mild anisotropy the Lam´e parameter relatedto the volumetric response is bounded. Locking does however occur in the inextensible limit.There have been few studies of transverse isotropy in the context of development of new finiteelement and related methods. The present study and the work cited above constitute two newcontributions. Further work is in progress on alternative formulations such as, for example, theuse of discontinuous Galerkin methods. The extension to problems involving nonlinear materialbehaviour and large deformations is also in progress. It would be of interest to investigate the22xtension of the work presented here to include higher order VEMs as well as problems in threedimensions.
Acknowledgement
This work was carried out with support from the National Research Foundation of South Africa,through the South African Research Chair in Computational Mechanics. The authors acknowl-edge with thanks this support.
References [1] Ted Belytschko, Wing Kam Liu, Brian Moran, and Khalil Elkhodary.
Nonlinear FiniteElements for Continua and Structures . Wiley, 2014.[2] P. Wriggers.
Nonlinear Finite Element Methods . Springer, Berlin, Heidelberg, New York,2008.[3] D. Boffi, F. Brezzi, and M. Fortin.
Mixed Finite Element Methods and Applications .Springer, New York, 2013.[4] T.J.R. Hughes.
The Finite Element Method. Linear Static and Dynamic Finite ElementAnalysis . Prentice-Hall, Englewood Cliffs, New Jersey, 1987.[5] D.N. Arnold, F. Brezzi, B. Cockburn, and L.D. Marini. Unified analysis of discontinuousGalerkin methods for elliptic problems.
SIAM Journal on Numerical Analysis , 39:1749–1779, 2002.[6] B. Grieshaber, A.T. McBride, and B.D. Reddy. Uniformly convergent interior penaltymethods using multilinear approximations for problems in elasticity.
SIAM Journal onNumerical Analysis , 53:2255–2278, 2015.[7] P. Hansbo and M.G. Larson. Discontinuous Galerkin methods for incompressible and nearlyincompressible elasticity by Nitsche’s method.
Computer Methods in Applied Mechanics andEngineering , 191:1895–1908, 2002.[8] T. P. Wihler. Locking-free DGFEM for elasticity problems in polygons.
IMA Journal ofNumerical Analysis , 24(1):45–75, jan 2004.[9] L. Beir˜ao da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L.D. Marini, and A. Russo. Ba-sic principles of virtual element methods.
Mathematical Models and Methods in AppliedSciences , 23(01):199–214, 2013. 2310] L. Beir˜ao da Veiga, F. Brezzi, L.D. Marini, and A. Russo. The hitchhiker’s guide to thevirtual element method.
Mathematical Models and Methods in Applied Sciences , 24:1541–1573, 2014.[11] A.L. Gain, C. Talischi, and G.H. Paulino. On the virtual element method for three-dimensional linear elasticity problems on arbitrary polyhedral meshes.
Computer Methodsin Applied Mechanics and Engineering , 282:132–160, 2014.[12] H. Chi, L. Beir˜ao da Veiga, and G. Paulino. Some basic formulations of the virtual ele-ment method (VEM) for finite deformations.
Computer Methods in Applied Mechanics andEngineering , 318:148–192, 2017.[13] P. Wriggers, B.D. Reddy, W. Rust, and B. Hudobivnik. Efficient virtual element formula-tions for compressible and incompressible finite deformations.
Computational Mechanics ,60:253–268, 2017.[14] E. Artioli, L. Beir˜ao da Veiga, C. Lovadina, and E. Sacco. Arbitrary order 2D virtual ele-ments for polygonal meshes: part II, inelastic problem.
Computational Mechanics , 60:643–657, 2017.[15] L. Beir˜ao da Veiga, C. Lovadina, and D. Mora. A virtual element method for elasticand inelastic problems on polytope meshes.
Computer Methods in Applied Mechanics andEngineering , 295:327–346, 2015.[16] P. Wriggers and B. Hudobivnik. A low order virtual element formulation for finite elasto-plastic deformations.
Computer Methods in Applied Mechanics and Engineering , 327:459–477, 2017.[17] P. Wriggers, W. Rust, and B.D. Reddy. A virtual element method for contact.
Computa-tional Mechanics , 58:1039–1050, 2016.[18] P. Wriggers, B. Hudobivnik, and J. Korelc. Efficient low order virtual elements foranisotropic materials at finite strains. In
Advances in Computational Plasticity , Computa-tional Methods in Applied Sciences, vol. 46. Springer, Berlin, 2018.[19] F. Auricchio, G. Scalet, and P. Wriggers. Fiber-reinforced materials: finite elements for thetreatment of the inextensibility constraint.
Computational Mechanics , 60:905–922, 2017.[20] Faraniaina Rasolofoson, BJ Grieshaber, and B Daya Reddy. Finite element approximationsfor near-incompressible and near-inextensible transversely isotropic bodies.
InternationalJournal for Numerical Methods in Engineering, , Advanced online publication, 2018.[21] G.E. Exadaktylos. On the constraints and relations of elastic constants of transverselyisotropic geomaterials.
International Journal of Rock Mechanics and Mining Sciences ,38:941–956, 2001. 2422] E. Artioli, L. Beir˜ao da Veiga, C. Lovadina, and E. Sacco. Arbitrary order 2D virtualelements for polygonal meshes: part I, elastic problem.