A three-dimensional nodal-based implementation of a family of discontinuous Galerkin methods for elasticity problems

doi:10.1016/j.compstruc.2008.11.009

Computers & Structures

Volume 87, Issues 3–4, February 2009, Pages 141-150

https://doi.org/10.1016/j.compstruc.2008.11.009 Get rights and content

Abstract

Four discontinuous Galerkin (DG) methods are proposed to enrich the resource of modeling elasticity problems as they are volume locking-free and allow hanging nodes in meshing. A detailed finite element formulation of these DG methods is presented. For implementation, we coded a three-dimensional nodal-based DG program in which the conventional nodal-based pure displacement finite element codes are fully exploited. The robustness and accuracy of each DG method are demonstrated and compared with mixed methods through solving a rubber beam problem. The coupled use of DG with continuous elements is proposed for some practical applications.

Introduction

In general, the popular pure displacement-based continuous Galerkin (CG) methods work well for linear elasticity problems with normal materials. However, their performance can be detrimental for cases where materials are highly incompressible [5]. In these cases, CG provides a smaller displacement solution, called volume locking. Methods for handling the volume locking include mixed finite element methods, reduced integrations, nonconforming methods, and stabilized finite elements [3], [6], [12], [14], [21], [22]. The locking may be completely removed using mixed finite element methods. However, the selection of the spaces for displacement and pressure is not independent and they must satisfy the BB condition [2]. Consequently, only 27-nodal elements with linear pressure variations (Q2/P1) and 27-nodal Taylor–Hood elements with continuous pressures (Q2/Q1) are the lowest order 3D elements that strictly satisfy the BB condition and exhibit optimal convergent rates [5]. The reduced integration method greatly reduces the accuracy. The lack of the coercivity in lower order nonconforming elements results in divergent solutions. Intensive trials for searching for an appropriate stabilized parameter required for stabilized finite element methods are impractical. We propose a family of DG methods for some practical applications to handle the locking for highly incompressible elasticity.

The main disadvantage for DG methods is the great increase in unknowns. Considering their many advantages over CG methods, however, DG methods can still be proposed for some practical applications. First, the feature of momentum conservation satisfied locally and weakly is very important to achieve a satisfactory stress result for elasticity problems with nearly incompressible materials. Second, DG allows hanging nodes that greatly facilitate the adaptive refinement on an original mesh in that each element can be refined or coarsened arbitrarily and independently, which could be a challenge for CG methods. Third, it has a potential to facilitate parallel computing. Fourth, it is a natural method to model some localization problems where the displacement field is discontinuous. Finally, the explosion in unknowns of DG methods can be alleviated by a coupled use of DG with CG methods in some cases, as explained in the section of numerical study.

The idea of DG finite element methods originated from Nitsche’s work in 1971 [15]. Rather than enforcing the Dirichlet boundary condition strongly in CG methods, he did it weakly. In 1976, the continuity of the stress in elliptic equations and the flux in parabolic equations across the interior faces between elements was weakly enforced by Douglas and Dupont [11]. The Interior Penalty Galerkin (IPG or SIPG) finite element method was first derived and formulated for second order elliptic equations by Wheeler [24] and Wheeler and Percell [19]. In their DG formulations, the stiffness contributed from interior faces is obtained by using the average of tractions and jumps on primary variables across the faces between elements. Also, the interior continuity condition of the primary variable is controlled by using a penalty method to add a penalty term involving the jump on the primary variable and its corresponding test function. By adjusting the penalty parameter, an accurate result may be obtained. If the penalty parameter is set to infinity, DG coincides with CG. The IPG method was exploited by Arnold [1] to solve parabolic equations and nonlinear elliptic problems. In addition, the penalty methods should be understood here. In CG framework, similar penalty techniques have been used to enforce the non-penetration condition for contact mechanics problems [16] and the rotation continuity condition for beam problems [4]. For example, in a popular and efficient pipe elbow element originally formulated by Bathe and Almelda [4], the jump of the first derivative of the continuous ovalization displacement on the interior jointed element face is penalized through incorporating a penalty energy function into the potential, in which the accuracy of the slope continuity across the joints is able to be improved by adjusting the penalty parameter.

Since finite element methods were developed, much research work has been focused on CG methods. Recently, DG methods have become more popular and have been exploited to deal with the challenges that CG methods are not able to handle easily [9]. A non-symmetric DG method proposed by Oden, Babuska and Baumann (OBB) [17] to solve diffusion problems results in a positive definite stiffness matrix. Therefore, the method is more stable and robust than the symmetric DG formulation, and it exhibits a property of local mass conservation at the element level, which is very important for solving convection and/or diffusion problems. In addition, Cockburn and Shu [8] developed the local DG (LDG), another family of DG which is different from the primal DG. This LDG was extended by Cockburn and Dawson [7] for convection–diffusion equations in multi-dimensions. In the LDG method, besides the primary variable, the flux also appears as an unknown.

Because of the symmetric advantage of the stiffness matrix, SIPG is often applied to problems that do not have strict stability requirements. SIPG was also used by Hansbo and Larson [13] to solve 2-D nearly incompressible elasticity and Stokes flow problems. Based on both IPG and OBB, a Non-Symmetric Interior Penalty Galerkin (NIPG) formulation was presented by Wheeler et al. [20] to solve elasticity problems. This formulation exhibits more stable and robust behavior than SIPG. Compared to OBB, NIPG allows one to adjust penalty parameters to achieve more accurate results. A new scheme, the Incomplete Interior Penalty Galerkin (IIPG) method proposed by Dawson et al. [10], has the best performance in terms of stability and accuracy of solutions for solving fluid problems. However, IIPG has not been studied for elasticity problems.

Recent effort in modeling elasticity problems by using DG methods [20], [13], [26] focused mostly on theoretical work on the convergence, stability, and locking free property of DG methods and on 2D implementations. The lack of effort in comparing the performance of different DG methods for 3D problems is obvious in the literature. In this paper, we compare the robustness and accuracy of four DG methods with the mixed method applied to 3D elasticity problems with highly incompressible materials. This work is important for selections on DG methods for practical engineering applications. Another important factor affecting DG application to practical problems is their implementations currently dominated by non-nodal based programs. Although these implementations are easy for full DG discretizations, it is difficult to use the coupling between continuous and discontinuous elements. In fact, the popular commercial finite element codes are nodal-based programs for CG methods. For our proposed four DG methods, we present a general 3D nodal-based implementation that is able to fully exploit a conventional 3D elasticity code and naturally handle the coupling of continuous and discontinuous elements. Adding interface stiffness into a CG code, which is actually computed using the same shape functions as CG programs, is only the additional work for this DG implementation. In Section 2, we formulate a family of DG methods for elasticity problems. In Section 3, we describe a 3D nodal-based implementation of DG methods and describe the procedures to construct interface stiffness. In Section 4, we compare the robustness and accuracy of our four DG methods with a mixed method using a nearly incompressible beam problem. The coupled use of DG with CG methods is also proposed in this section. Conclusions are provided in the final section.

Section snippets

DG formulation on elasticity problem

Let a body occupy a domain $Ω \subset R^{3}$ shown in Fig. 1. $Ω$ has a boundary surface $\partial Ω$ which is further divided into $Γ_{u}$ and $Γ_{t}$ with $\partial Ω = {\bar{Γ}}_{u} \cup {\bar{Γ}}_{t}, Γ_{u} \cap Γ_{t} = \emptyset,$ where the prescribed displacement $\bar{u} (x) \in H^{1} (Γ_{u})$ is given on the part $Γ_{u}$ of the boundary and the surface traction $\bar{t} (x) \in L^{2} (Γ_{t})$ is given on the remainder $Γ_{t}$ of the boundary. We denote the body force by $f (x) \in L^{2} (Ω)$ , the Cauchy stress tensor by $σ$ and the displacement field by variable u. The linear momentum equation for three-dimensional linear elasticity is $\nabla \cdot σ + f =$

A nodal-based DG implementation

The popular CG finite element codes are nodal-based. In general, a CG program uses geometry data containing the coordinates of nodes. The profile of an element is defined by its nodes. These nodes completely determine the order of the element. Moreover, as the solutions are just the values at nodes, they can be directly used for output and visualization without any further post-processing. The nodal-based CG computer programs for linear elasticity have been thoroughly implemented, tested, and

Numerical study

We consider a three-dimensional version of a cantilever beam problem with similar data used in [5]. This beam is shown in Fig. 5 but drawn only in y–z plane. The displacement boundary conditions are represented by the rollers. The constrains in the x-direction are on only the surfaces at $z = \pm 20 mm$ . The beam with 1 mm width is loaded by a uniform pressure on its top surface. We studied two cases, Poisson’s ratio $ν = 0.3$ and $ν = 0.499$ . In the following subsections, the more accurate results are

Conclusions

We have presented a formulation and implementation of a family of DG methods for solving elasticity problems. The key procedure for DG methods is to add the interface stiffness into a nodal-based CG program. The normal matrix for stiffness computation of a general curved interface has been introduced in this paper, which leads to a practical compacted Voigt implementation. Our numerical study demonstrates: (a) All proposed DG methods greatly improve the stress prediction when compared to CG

References (23)

I. Babuska et al.
The Babuska–Brezzi condition and the patch test: an example
Comput Methods Appl Mech Eng
(1997)
C.N. Dawson et al.
Compatible algorithms for coupled flow and transport
Comput Methods Appl Mech Engrg
(2004)
P. Hansbo et al.
Discontinuous Galerkin method for incompressible and nearly incompressible elasticity by Nitche’s methods
Comput Methods Appl Mech Eng
(2002)
D.S. Malkus et al.
Mixed finite element methods-reduced and selective integration techniques: a unification of concepts
Comput Methods Appl Mech Eng
(1978)
J.T. Oden et al.
A discontinuous hp finite element method for diffusion problems
J Comput Phys
(1998)
L.C. Wellford et al.
A theory of discontinuous finite element approximations for the analysis of shock waves in nonlinear elastic materials
J Comput Phys
(1975)
D.N. Arnold
An interior penalty finite element method with discontinuous elements
SIAM J Numer Anal
(1982)
G.A. Baker
Finite element methods for elliptic equations using nonconforming elements
Math Comput
(1977)
K.J. Bathe et al.
A simple and effective Pipe elbow element-interaction effects
J Appl Mech
(1982)
K.J. Bathe
Finite element procedures
(1996)

F. Brezzi et al.

Mixed and hybrid finite element methods

(1994)

Cited by (34)

A DG-based interface element method for modeling hydraulic fracturing in porous media
2020, Computer Methods in Applied Mechanics and Engineering
Citation Excerpt :
In [44], for pure solids undergoing finite deformation, DG formulations were performed and verified for both hypoelastoplasticity and hyperelastoplasticity problems. While all DG methods could be proposed to solve this hydraulic fracturing problem, we adopt the Incomplete Interior Penalty Galerkin (IIPG) [38] as this IIPG scheme exhibits a more robust performance in solving nearly incompressible elasticity problems [42] and a mathematically complete and consistent formulation can be nicely achieved for plasticity problems using this method in [44]. Moreover, the success of DG formulations is only observed in the IIPG scheme for poromechanics problems where the solid and fluid flow fields are coupled [43,46].
Fracture propagation coupled with fluid flow in porous media has very important applications in fracking in oil and gas reservoirs, corrosion in ceramics, and degradation of human bones. Modeling fracture failure in porous media governed by the coupled solid and fluid field equations is challenging and time-consuming. This paper extends discontinuous Galerkin (DG) finite element methods to model crack openings in porous media through exploiting an easy construction of interfaces for potential crack paths by locally breaking continuous elements. Consequently, a finite element mesh for fracture apertures is completely and gracefully constructed using DG interface objects as well as the definitions of jumps of displacements across the DG interfaces defined for bulk matrices. Furthermore, rigid cohesive laws often adopted for crack openings can be naturally implemented in DG formulations without introducing artificial stiffness, which may eventually improve the performance for implicit formulations in handling crack openings in porous media. In this work, we perform a consistent, fully implicit, and fully coupled hybrid DG/continuous finite element formulation for three field equations including the solid, bulk fluid, and fluid in fracture apertures resulting from crack openings. We verify our DG formulation and implementation using the Khristianovich-Geertsma-DeKlerk analytical solutions. Finally, we further demonstrate a good performance of the proposed DG method by modeling fracking in an oil reservoir containing two nature fractures.
Convergence in the incompressible limit of new discontinuous Galerkin methods with general quadrilateral and hexahedral elements
2020, Computer Methods in Applied Mechanics and Engineering
Citation Excerpt :
In contrast, DG primal formulations have been established as having optimal performance independent of material parameters for meshes of triangular elements [4,8–11], but not for meshes of quadrilateral elements. In a numerical investigation, Liu et al. [12] consider the matter of locking in the context of low-order hexahedral elements. This work, on a specific benchmark problem, shows the under-performance of the SG method and the superiority of three well-known primal IP methods, NIPG, SIPG and IIPG (Nonsymmetric, Symmetric and Incomplete Interior Penalty Galerkin methods, developed in [13–16]), as well as of the method of Oden, Babus̆ka and Baumann [17] (known as OBB), a penalty-free version of NIPG.
Standard low-order finite elements, which perform well for problems involving compressible elastic materials, are known to under-perform when nearly incompressible materials are involved, commonly exhibiting the locking phenomenon. Interior penalty (IP) discontinuous Galerkin methods have been shown to circumvent locking when simplicial elements are used. The same IP methods, however, result in locking on meshes of quadrilaterals. The authors have shown in earlier work that under-integration of specified terms in the IP formulation eliminates the locking problem for rectangular elements. Here it is demonstrated through an extensive numerical investigation that the effect of using under-integration carries over successfully to meshes of more general quadrilateral elements, as would likely be used in practical applications, and results in accurate displacement approximations. Uniform convergence with respect to the compressibility parameter is shown numerically. Additionally, a stress approximation obtained here by postprocessing shows good convergence in the incompressible limit.
A stabilized finite element method for enforcing stiff anisotropic cohesive laws using interface elements
2019, Computer Methods in Applied Mechanics and Engineering
We present a stabilized finite element method that generalizes Nitsche’s method for enforcing stiff anisotropic cohesive laws with different normal and tangential stiffness. For smaller values of cohesive stiffness, the stabilized method resembles the standard method, wherein the traction on the crack surface is enforced as a Neumann boundary condition. Conversely, for larger values of cohesive stiffness, the stabilized method resembles Nitsche’s method, wherein the cohesive law is enforced as a kinematic constraint. We present several numerical examples, in two-dimensions, to compare the performance of the stabilized and standard methods. Our results illustrate that the stabilized method enables accurate recovery of crack-face tractions for stiff isotropic and anisotropic cohesive laws; whereas, the standard method is less accurate due to spurious traction oscillations. Also, the stabilized method could mitigate spurious sensitivity of load–displacement results to displacement increment in mixed-mode fracture simulation, owing to its stability and accuracy.
A variational multiscale discontinuous Galerkin formulation for both implicit and explicit dynamic modeling of interfacial fracture
2019, Computer Methods in Applied Mechanics and Engineering
A stabilized Discontinuous Galerkin method is presented for modeling interfacial debonding under dynamic conditions. The method’s distinctive feature is its treatment of the inelastic gap at the interface as an internal variable that is evolved through a traction–separation constitutive model, enabling initially perfect adhesion to be captured. Numerical stability is derived using Variational Multiscale ideas, and comparisons are drawn with both the intrinsic cohesive zone method and hybrid Discontinuous Galerkin-cohesive zone method to expose subtle differences. The proposed formulation is applied to model the nonlinear response of a composite unit cell under impact loading. Both explicit and implicit time integration schemes provide similar accuracy and stability, and the critical time step for explicit analyses is only slightly reduced below the value for the continuous Galerkin mesh without an interface. Further studies of mesh refinement and adjustments to the initial elastic interface stiffness in the reference cohesive zone model reveal the improved accuracy of the Variational Multiscale Discontinuous Galerkin method on coarser grids.
On the algorithmic and implementational aspects of a Discontinuous Galerkin method at finite strains
2015, Computers and Mathematics with Applications
In this work, algorithmic modifications are proposed and analyzed for a recently developed stabilized finite strain Discontinuous Galerkin (DG) method. The distinguishing feature of the original method, referred to as VMDG, is a consistently derived expression for the numerical flux and stability tensor that account for evolving material and geometric nonlinearity in the vicinity of the interface. Herein, the proposed modifications involve simplifications to the residual force vector and tangent stiffness matrix of the VMDG method that lead to formulations similar to other existing DG methods but retain the enhanced definition for the stability parameters. The primary objective is to reduce the costs associated with implementing the method as well as executing simulations while retaining accuracy and flexibility, thereby making the formulation amenable to boarder material classes such as inelasticity. Each simplification has associated implications on the mathematical and algorithmic properties of the method, such as $L^{2}$ convergence rate, solution accuracy, continuity enforcement, and stability of the nonlinear equation solver. These implications are carefully quantified and assessed through a comprehensive numerical performance study. The range of two and three dimensional problems under consideration involves both isotropic and anisotropic materials. Both triangular and quadrilateral element types are employed along with $h$ and $p$ refinement. The ability of the proposed methods to produce stable and accurate results for such a broad class of problems is highlighted.
Volume-averaged nodal projection method for nearly-incompressible elasticity using meshfree and bubble basis functions
2015, Computer Methods in Applied Mechanics and Engineering
We present a displacement-based Galerkin meshfree method for the analysis of nearly-incompressible linear elastic solids, where low-order simplicial tessellations (i.e., 3-node triangular or 4-node tetrahedral meshes) are used as a background structure for numerical integration of the weak form integrals and to get the nodal information for the computation of the meshfree basis functions. In this approach, a volume-averaged nodal projection operator is constructed to project the dilatational strain into an approximation space of equal- or lower-order than the approximation space for the displacement field resulting in a locking-free method. The stability of the method is provided via bubble-like basis functions. Because the notion of an ‘element’ or ‘cell’ is not present in the computation of the meshfree basis functions such low-order tessellations can be used regardless of the order of the approximation spaces desired. First- and second-order meshfree basis functions are chosen as a particular case in the proposed method. Numerical examples are provided in two and three dimensions to demonstrate the robustness of the method, its ability to avoid volumetric locking in the nearly-incompressible regime, and its improved performance when compared to the MINI finite element scheme on the simplicial mesh.

View all citing articles on Scopus

View full text

A three-dimensional nodal-based implementation of a family of discontinuous Galerkin methods for elasticity problems

Abstract

Introduction

Section snippets

DG formulation on elasticity problem

A nodal-based DG implementation

Numerical study

Conclusions

Comput Methods Appl Mech Eng

Comput Methods Appl Mech Engrg

Comput Methods Appl Mech Eng

Comput Methods Appl Mech Eng

J Comput Phys

J Comput Phys

An interior penalty finite element method with discontinuous elements

SIAM J Numer Anal

Finite element methods for elliptic equations using nonconforming elements

Math Comput

A simple and effective Pipe elbow element-interaction effects

J Appl Mech

Finite element procedures

Mixed and hybrid finite element methods