- Email: [email protected]
A combined virtual element method and the scaled boundary finite element method for linear elastic fracture mechanics
Engineering Analysis with Boundary Elements 113 (2020) 9–16
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A combined virtual element method and the scaled boundary finite element method for linear elastic fracture mechanics Dibyendu Adak a, ALN Pramod a, Ean Tat Ooi b, Sundararajan Natarajan a,∗ a b
Integrated Modelling and Simulation Lab, Department of Mechanical Engineering, Indian Institute of Technology Madras, Chennai 600036, India School of Science, Engineering & Information Technology, Federation University Australia, Ballarat VIC3350, Australia
a r t i c l e
i n f o
Keywords: Stress intensity factor Virtual element method Scaled boundary finite element method Polygonal elements
a b s t r a c t In this paper, we propose a framework that combines the recently introduced virtual element method (VEM) and the scaled boundary finite element method (SBFEM) to evaluate the fracture parameters. The domain is discretized with arbitrary polygons and the element that contains the crack tip is treated within the framework of the SBFEM. This facilitates a semi-analytical treatment of the crack tip singularity allowing the fracture parameters are estimated directly from the definition. The VEM is employed for the rest of the domain. The salient feature of the VEM is that the terms in the stiffness matrix are computed without requiring higher order quadrature schemes. As both the methods satisfy partition of unity and the compatibility condition, the matrices are assembled as in the conventional FEM. The accuracy of the proposed formulation is demonstrated with two standard benchmark examples. The proposed VEM-SBFEM framework yields accurate results.
1. Introduction One of the computational challenges that the finite element method (FEM) faces is in modelling evolving discontinuities. Discontinuities such as cracks and notches require local mesh adaptation with special elements to treat the singular fields. Moreover, re-meshing [1,2] may be required as the discontinuities evolve in time. Another difficulty is that the conventional FEM is restricted to simplex elements, such as triangles and quadrilaterals. This poses a challenge for local mesh adaption. Thanks to the introduction of the partition of unity methods (PUM), such as the eXtended FEM (XFEM)/partition of unity FEM (PUFEM)/generalized FEM (GFEM) [3–6], element free Galerkin method (EFG) [7] and the polygonal FEM (PFEM) [8–10] which alleviates the meshing burden by detaching the geometry and the discretization. This allows the geometry to be represented independent of the underlying discretization. On the other hand, the PFEM allows the elements to take arbitrary number of sides and shapes. Both the approaches aim at relaxing the meshing burden and the topology constraint imposed by the conventional FEM. The PUM captures the discontinuity by supplementing the conventional FE basis with a priori known additional functions that capture the local behavior, whilst the PFEM allows the elements to take arbitrary shapes and can have arbitrary number of sides. Since the inception of the XFEM and the PFEM, these approaches have been applied to wider variety of problems. For a review on XFEM and PFEM, interested readers are referred to [3,11] and the references
∗
therein. The flexibility offered by the PFEM has led researchers to develop methods with polygonal discretizations, viz., virtual element method (VEM) [10,12–15], virtual node method [16], smoothed finite element method [17,18], scaled boundary finite element method (SBFEM) [19–21] and the discontinuous Galerkin method [22]. The VEM was recently developed in [10,12–15] and could be considered as a generalization of the FEM over arbitrary polytopes. Unlike the conventional FEM, the VEM does not require an explicit form for the basis functions to compute the discrete variational formulation. The VEM also alleviates the numerical integration difficulty encountered in the conventional polygonal FEM. As the method does not require the knowledge of the shape functions in explicit form, the implementation is computationally less intensive. However, since the VEM space contains only polynomial basis, the VEM cannot accurately capture the asymptotic crack tip fields in the vicinity of the crack tip. On another front, the scaled boundary finite element method (SBFEM) introduced by Song and Wolf has emerged as an attractive alternate to model problems with singularities [23,24]. It combines the best features of the FEM and the boundary element method (BEM). In [25] a recent review on the application of the SBFEM to linear elastic fracture mechanics can be found. Recently, Hyun and co-workers [26,27], proposed extended PFEM (XPFEM) by combining the XFEM with the PFEM. Within this formulation, the background discretization consists of arbitrary polygons and the XFEM allows for capturing the discontinuities without requiring a conforming mesh. However, the method requires a priori knowledge of the enrichment functions and special numerical integration schemes
Corresponding author. E-mail address: [email protected] (S. Natarajan).
https://doi.org/10.1016/j.enganabound.2019.12.008 Received 15 April 2019; Received in revised form 15 October 2019; Accepted 16 December 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
The corresponding weak form is: find 𝐮 ∈ 𝒰 such that
Ωh
Γt
∀𝐯 ∈ 𝒱 , 𝑎 (𝐮, 𝐯) =
Γc
Fig. 1. General elastic body with an internal discontinuity, Neumann and Dirichlet boundary conditions.
𝝈(𝐮) ∶ 𝜺(𝐯) d𝑉 , 𝓁(𝐯) =
∫Ω
𝐛 ⋅ 𝐯 d𝑉 +
∫Γ𝑡
𝐭̂ ⋅ 𝐯 d𝑆,
(3b)
where the space (Ω) includes the linear displacement fields. The domain is partitioned into elements Ωh . Using shape functions 𝜙a that span at least the linear space, we substitute vector-valued trial and test func∑ ∑ tions 𝐮ℎ = 𝑎 𝜙𝑎 𝐮𝑎 and 𝐯ℎ = 𝑏 𝜙𝑏 𝐯𝑏 , respectively, into Eq. (3) and apply a standard Galerkin procedure to obtain the discrete weak form: find 𝐮ℎ ∈ 𝒰 ℎ such that
have to be developed to integrate the enrichment functions. A different numerical formulation that also shows promise when applied to model fracture related problems is peridynamics. Peridynamics does not require an explicit description of the crack topology, an advantage that has attracted considerable amount of interest among researchers. However, peridynamics has a problem of spurious wave reflections and ghost forces between the particles [28,29]. The main objective of the paper is to combine the VEM with the SBFEM. The element that contains the crack tip is treated within the framework of the SBFEM. This facilitates a semi-analytical treatment of the crack tip singularity and the fracture parameters are estimated directly from the definition. The VEM is employed for the rest of the domain. The salient features of the proposed framework are: (a) does not require a priori knowledge of the enrichment functions, (b) the stress intensity factors are computed directly from the definition without requiring ‘any’ additional post processing, such as J-integral and (c) the SBFEM which requires a solution of an eigenvalue problem is restricted to a smaller region, thus increasing the computational efficiency. The rest of the manuscript is organized as follows: Section 2 presents the governing equation and the corresponding weak form. An overview of the VEM and the SBFEM is presented in Sections 3 and 4, respectively. Section 4 also presents the procedure to estimate the stress intensity factors directly from the definition. The accuracy of the estimated stress intensity factors for two standard problems in linear elastic fracture mechanics from the proposed framework are compared with available results in the literature in Section 5, followed by concluding remarks in the last section.
∀𝐯ℎ ∈ 𝒱 ℎ
𝑎(𝐮ℎ , 𝐯ℎ ) = 𝓁(𝐯ℎ ),
(5)
which leads to the following system of linear equations:
𝐊=
∑
𝐊𝐮= 𝐟 , ∑ 𝐊ℎ =
ℎ
𝐟=
∑
ℎ
𝐟 ℎ=
ℎ
∑
(6a) ∫Ωℎ (
𝐁T ℂ𝐁 d𝑉 ,
∫Ωℎ
ℎ
(6b) )
𝐍T 𝐛 d𝑉 +
∫Γℎ
𝐍T 𝐭̂ d𝑆 ,
(6c)
𝑡
where K is the assembled stiffness matrix, f the assembled nodal force vector, u the assembled vector of nodal displacements, N is the matrix of shape functions, ℂ is the constitutive matrix for an isotropic linear elastic material, and 𝐁 = 𝛁𝐍 is the strain-displacement matrix that is computed using the derivatives of the shape functions. In this study, the domain is assumed to be partitioned into non-overlapping elements with arbitrary number of sides. The main challenge is in computing the terms in the bilinear and linear form, as the shape functions over such arbitrary polygons are not simple polynomials [8,9]. 3. Basics of the virtual element method In this section, we present the virtual element formulation for the elastic deformation problem. Only important equations are presented here and for more detailed discussion and implementation, interested readers are referred to [30,31] and the references therein. The continuous variational formulation of the elastic deformation problem is given by Eq. (5). Since the virtual element space does not require explicit expression of the basis function, two projection operators 𝐸0 and 𝐸∇ are introduced to compute the bilinear and the linear form. Employing these projection operators, we project the virtual basis functions onto a computable polynomial subspace. The tensor valued L2 projection operator is defined as follows:
2. Governing equations and weak form In this section, the governing partial differential equations and the corresponding weak form for linear homogeneous elastic body containing a strong discontinuity (i.e., crack) is given. Consider a linear elastic body occupying Ω ⊂ IR2 bounded by a surface Γ, with an outward normal n and a discontinuity Γc as shown in Fig. 1. The boundary is assumed to admit decompositions: Γ = Γ𝑢 ∪ Γ𝑡 ∪ Γ𝑐 , where Γu is the Dirichlet boundary, Γt is the Neumann boundary and Γc , is the traction free surface representing the discontinuity. In the absence of interia and body forces, the boundary value problem for linear elasto-statics is: find 𝐮 ∶ Ω → ℝ2 , such that: in Ω on Γ𝑢 on Γ𝑡 , on Γ𝑐 ,
(3a)
where 𝜺 is the small strain tensor, and 𝒰 and 𝒱 are the displacement trial and test spaces: { } 𝒰 ∶= 𝐮(𝐱) ∈ [𝐶 0 (Ω)]𝑑 ∶ 𝐮 ∈ [(Ω)]𝑑 ⊆ [𝐻 1 (Ω)]𝑑 , 𝐮 = 𝐮̂ on Γ𝑢 , { } 𝒱 ∶= 𝐯(𝐱) ∈ [𝐶 0 (Ω)]𝑑 ∶ 𝐯 ∈ [(Ω)]𝑑 ⊆ [𝐻 1 (Ω)]𝑑 , 𝐯 = 𝟎 on Γ𝑢 ,
Γu
⎧ ⎪𝛁 ⋅ 𝝈 = 𝟎 ⎪𝐮 = 𝐮̂ ⎨𝝈 ⋅ 𝐧 = 𝐭̂ ⎪ ⎪𝝈 ⋅ 𝐧 = 𝟎 ⎩
∫Ω
𝑎(𝐮, 𝐯)= 𝓁(𝐯)
𝐸0 𝐐 ∶=
1 𝐐 |𝐸| ∫𝐸
∀ 𝐸 ∈ Ωℎ .
(7)
The global projection operator is defined as 0 (𝐐)|𝐸 = 𝐸0 (𝐐). In addition, we introduce a tensor valued elliptic type projection operator ∇ , defined on H1 (E)d as follows: for any 𝐯 ∈ ℎ we have 𝐸∇ 𝐯 ∈ 𝐏1 (𝐸) and it satisfies
(1)
⎧∇(𝐸∇ (𝐯)) ⎪ ⎨∑ ∇ 𝐸 (𝑙) ⎪ ⎩𝑙∈𝜕𝐸
where 𝛁 is the gradient operator, 𝝈 is the Cauchy stress tensor, n is the unit outward normal, u is the displacement field and 𝐮̂ and 𝐭̂ are the applied displacement and traction on the Dirichlet and Neumann boundary. For small displacements and strains, the strain displacement relation is given by: ) 1( 𝜺= 𝛁𝐮 + 𝛁𝐮T (2) 2
= 𝐸0 (∇𝐯|𝐸 ) =
∑
𝐯(𝑙)
(8)
𝑙∈𝜕𝐸
The local virtual element space is defined as ℎ,𝐸 ∶= {𝑣 ∈ [𝐻 1 (𝐸) ∩ 𝐶 0 (𝐸)]2 ∶ Δ𝐯 = 0 in 𝐸, 𝐯|𝐹 ∈ 𝐏1 (𝐹 ) ∀ 𝐹 ∈ 𝜕𝐸}. (9) 10
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
η
and the corresponding global virtual element space is defined as ℎ ∶= {𝐯 ∈ 𝐻10 (Ω)2 ∶ 𝐯|𝐸 ∈ ℎ,𝐸 }.
(10)
The expression for the stiffness matrix in the case of two dimensional linear elasticity can be written as: 𝐊𝐸 ℎ =
|𝐸|𝐖𝐶 ℂ𝐖T𝐶 ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟
Consist ency t erm=𝐊const
+ (𝐈 − 𝐏𝑝 )T 𝐒𝐸 (𝐈 − 𝐏𝑝 ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟
Line element
(x2,y2), η = 1
(x1,y1), η = −1
(11)
ξ=1
Stabilit y t erm=𝐊stab
S
where 𝐏𝑝 = 𝐏𝑅 + 𝐏𝐶
(12)
and 𝐏𝑅 = 𝐍𝑅 𝐖T𝑅 𝐏𝐶 = 𝐍𝐶 𝐖T𝐶 The NR and NC are expressed as: [ ] 1 0 (𝐱𝐼 − 𝐱 )2 (𝐍𝑅 )𝐼 = 0 1 −(𝐱𝐼 − 𝐱)1 [ ] (𝐱𝐼 − 𝐱 )1 0 (𝐱𝐼 − 𝐱 )2 (𝐍𝐶 )2 = 0 (𝐱𝐼 − 𝐱 )2 (𝐱𝐼 − 𝐱 )1
(13)
Radial Lines A
B
Scaling Centre
(14)
Fig. 2. Cracked polygon representation by the scaled boundary finite element method.
where xI is the coordinate of the node and 𝐱 is the polyhedron centroid. The block rows of WR and WC are expressed as: [ ] 1∕𝑛 0 (𝐪𝐼 )2 (𝐖𝑅 )𝐼 = 0 1∕𝑛 −(𝐪𝐼 )1 [ ] 2(𝐪𝐼 )1 0 (𝐪𝐼 )2 (𝐖𝐶 )𝐼 = (15) 0 2(𝐪𝐼 )2 (𝐪𝐼 )1
where N(𝜂) is the shape function matrix of the finite elements discretising the polygon boundary. The standard 1D Gauss–Lobatto shape functions or Lagrange shape functions can be used. In this study, we employ Lagrange shape functions. The displacements at a point in a polygon is approximated by: 𝐮(𝜉, 𝜂) = 𝐍(𝜂)𝐮(𝜉)
where n is the number of vertices and 𝐼 = 1, … , 𝑁 with N the number of nodes in the element and the subscript indicates the component of the associated vector and 𝐪𝐼 =
1 𝑁𝐼 𝐧 dΓ 2|𝐸| ∫
where u(𝜉) are the radial displacement functions. Substituting Eq. (18) in the definition of strain-displacement relations, the strains 𝜺(𝜉, 𝜂) in the scaled boundary coordinates are expressed as:
(16)
𝜕𝐸
𝜺(𝜉, 𝜂) = 𝐋𝐮(𝜉, 𝜂)
where 𝐒𝐸 = 𝛼𝐈 and 𝛼 = 𝛼 ∗ trace(|𝐸|𝐖𝐶 ℂ𝐖T𝐶 ) is a scaling coefficient and n is the unit outward normal vector. It can be seen that the computation of the stiffness matrix involves computing the matrices NR , NC , WR and WC . The calculation of the matrices WR and WC involves computing the surface integral of the basis functions.
(19)
where L is a linear operator matrix formulated in the scaled boundary coordinates as 𝐋 = 𝐛1 (𝜂)
4. Overview of the scaled boundary finite element method
𝜕 1 𝜕 + 𝐛 (𝜂) 𝜕𝜉 𝜉 2 𝜕𝜂
(20)
with
The VEM discussed in the last section can be applied to problems with strong discontinuity and singularity. However, to accurately capture the asymptotic fields at the crack tip, a very fine mesh in combination with singular elements at the crack tip is usually required. This poses additional difficulties when the crack evolves. Another possibility is to enrich the approximation space with functions that can capture the discontinuity and singularity [3,4]. In the literature, the latter method is referred to as the PUFEM. In this study, we propose to couple the VEM with the SBFEM to model problems with strong discontinuity and singularities. The SBFEM is a semi-analytical method and relies on defining a ‘scaling centre’ from which the entire boundary is visible. This is similar to the concept of ‘star convexity’. The boundary is divided into conventional linear finite elements, whilst the solution is sought analytically in the radial direction [32]. Moreover, by exploiting the special characteristics of the scaling centre, the stress intensity factors can be computed directly from their definition. When modelling a crack/notched surface the scaling centre is placed at the crack tip. The straight crack/notch edges are formed by scaling the nodes A and B on the boundary and the crack surfaces are not discretized (see Fig. 2). Displacement approximation The geometry of the element described by the coordinates on the boundary xb (𝜂) is expressed as: 𝐱𝑏 (𝜂) = 𝐍(𝜂)𝐱𝑏
(18)
𝐛1 (𝜂) =
⎡ 𝑦 (𝜂) 1 ⎢ 𝜂 ,𝜂 0 |𝐉(𝜂)| ⎢ ⎣−𝑥𝜂 (𝜂),𝜂
𝐛2 (𝜂) =
⎡−𝑦 (𝜂) 1 ⎢ 𝜂 0 |𝐉(𝜂)| ⎢ ⎣ 𝑥𝜂 (𝜂)
0 ⎤ −𝑥𝜂 (𝜂),𝜂 ⎥ ⎥ 𝑦𝜂 (𝜂),𝜂 ⎦ 0 ⎤ 𝑥𝜂 (𝜂) ⎥ ⎥ −𝑦𝜂 (𝜂)⎦
(21)
Using Eq. (19) and Hooke’s law 𝝈 = ℂ𝜺, the stresses 𝝈(𝜉, 𝜂) is expressed as ( ) 𝝈(𝜉, 𝜂) = ℂ 𝐁1 (𝜂)𝐮(𝜉),𝜉 + 𝜉 −1 𝐁2 (𝜂)𝐮(𝜉)
(22)
where ℂ is the material constitutive matrix. B1 (𝜂) and B2 (𝜂) are the strain-displacement matrices in the SBFEM and is defined as 𝐁1 (𝜂) = 𝐛1 (𝜂)𝐍(𝜂) 𝐁2 (𝜂) = 𝐛2 (𝜂)𝐍(𝜂),𝜂
(23)
By following the procedure outlined in [33,34], the following ODE is obtained: 𝐄0 𝜉 2 𝐮(𝜉),𝜉𝜉 + (𝐄0 + 𝐄T1 − 𝐄1 )𝜉𝐮(𝜉),𝜉 − 𝐄2 𝐮(𝜉) = 0 where E0 , E1 and E2 are the coefficient matrices given by:
(17) 11
(24)
D. Adak, A. Pramod and E.T. Ooi et al.
𝐄0 = 𝐄1 = 𝐄2 =
∫𝜂 ∫𝜂 ∫𝜂
Engineering Analysis with Boundary Elements 113 (2020) 9–16
𝐁1 (𝜂)T ℂ𝐁1 (𝜂)|𝐉(𝜂)|𝑑𝜂, 𝐁2 (𝜂)T ℂ𝐁1 (𝜂)|𝐉(𝜂)|𝑑𝜂, 𝐁2 (𝜂)T ℂ𝐁2 (𝜂)|𝐉(𝜂)|𝑑𝜂.
(25)
y
x
and the determinant of the Jacobian matrix is: |𝐉(𝜂)| = 𝑥𝑏 (𝜂)𝑦𝑏 (𝜂),𝜂 − 𝑦𝑏 (𝜂)𝑥𝑏 (𝜂),𝜂
O
(26)
The coefficient matrices are evaluated element-by-element on the polygon boundary and assembled over a polygon. This process is similar to the standard FE procedure of assemblage. Eq. (24) is a homogeneous second-order ordinary differential equation. Its solution is obtained by introducing the variable 𝝌(𝜉) { } 𝐮(𝜉) 𝝌(𝜉) = (27) 𝐪(𝜉)
ro
Fig. 3. A cracked domain modelled by the SBFEM and the definition of local coordinate system, where the ‘black’ dots represent the nodes.
where q(𝜉) is the internal load vector 𝐪(𝜉) = 𝐄0 𝜉𝐮(𝜉),𝜉 + 𝐄T1 𝐮(𝜉)
Lo
(28)
The boundary nodal forces are related to the displacement functions by:
From Eq. (39), the stiffness matrix K can be identified to be given by the expression
𝐟 = 𝐪(𝜉 = 1) = (𝐄0 𝜉𝐮(𝜉),𝜉 + 𝐄T1 𝐮(𝜉))|𝜉=1
𝐊 = 𝚽q 𝚽−1 u
(29)
Remark 4.1. The stiffness computed by employing the SBFEM is positive definite and symmetric. Hence, the stiffness matrix can be assembled in the conventional FEM approach. A simple Matlab ○R function is given in [35] to compute the stiffness matrix using the SBFEM.
This allows Eq. (24) to be transformed into a first order ordinary differential equation with twice the number of unknowns in an element as: 𝜉𝝌(𝜉),𝜉 = −𝐙𝝌(𝜉) where Z is the Hamiltonian matrix [ ] 𝐄−1 𝐄T1 −𝐄−1 0 0 𝐙= 𝐄1 𝐄−1 𝐄T1 − 𝐄2 −𝐄1 𝐄−1 0 0
(30)
Calculation of the stress intensity factors A unique feature of the SBFEM is that stress singularities, if present, are analytically represented in the radial displacement functions u(𝜉). When a crack is modelled by a polygon with its scaling centre chosen at the crack tip in Fig. 2, some of (s) the eigenvalues 𝚲(s) n ⊂ 𝚲n satisfy −1 < 𝚲n < 0. These eigenvalues lead to singular stresses at the crack tip i.e. when 𝜉 → 0. Using 𝚲(s) n , the singular stress field 𝝈 (s) (𝜉, 𝜂) can be defined as:
(31)
An eigenvalue decomposition of Z is performed. The blocks of eigenvalues and transformation matrices necessary are: [ ] [ ] 𝚽 𝚽u 𝐙 u = 𝚲 (32) 𝚽q 𝚽q n ( ) In Eq. (32), 𝚲n = diag 𝜆1 , 𝜆2 , … , 𝜆𝑛 contains the eigenvalues with negative real part. 𝚽u and 𝚽q are the corresponding transformation matrices of 𝚲n . They represent the modal displacements and forces, respectively. The general solution of Eq. (30) is given by: 𝐮(𝜉) = 𝚽u 𝜉 −𝚲n 𝐜
(33)
𝐪(𝜉) = 𝚽q 𝜉 −𝚲n 𝐜
(34)
(s)
− 𝚲n 𝝈 (s) (𝜉, 𝜂) =𝚿(s) 𝜎 (𝜂(𝜃))𝜉
𝚽(s) u
Taking the derivative of u(𝜉) with respect to 𝜉 and substituting into Eq. (22), the stress field 𝝈(𝜉, 𝜂) in the scaled boundary coordinates can be expressed as:
where the stress mode 𝚿𝜎 (𝜂) is defined as: ( ) 𝚿𝜎 (𝜂) = ℂ −𝐁1 (𝜂)𝚽u 𝚲n + 𝐁2 (𝜂)𝚽u
(38)
{
𝐾𝐼 𝐾𝐼𝐼
=
(42)
c(s)
}
⎧ (s) (s) ⎫ √ ⎪𝚿𝜎𝑦𝑦 (𝜂(𝜃 = 0))𝐜 ⎪ = 2𝜋𝐿𝑜 ⎨ ⎬ (s) ⎪𝚿𝜏𝑥𝑦 (𝜂(𝜃 = 0))𝐜(s) ⎪ ⎩ ⎭
(44)
5. Numerical examples In this section, the accuracy and the convergence properties of the proposed framework is demonstrated by solving two standard benchmark problems within the framework of linear elastic fracture mechanics. The discretization is based on centroid Voronoi tessellation. The two
The stiffness matrix of an element is obtained by first substituting Eq. (35) into Eq. (34) at 𝜉 = 1. This results in: 𝐟 = 𝚽q 𝚽−1 u 𝐮b
mode
Substituting the stress components in Eq. (41) at angle 𝜃 = 0 into Eq. (43) and using the relation 𝜉 = 𝑟∕𝐿𝑜 at 𝜃 = 0, the stress intensity factors are
(36)
(37)
(41) 𝚿(s) 𝜎 (𝜂(𝜃))
In Eq. (42) ⊂ 𝚽u and ⊂ c, contain the displacement modes and the integration constants corresponding to 𝚲(s) n . It can be discerned from Eq. (41) that 𝚲(s) n leads to singular stresses at the crack tip. This enables the stress intensity factors to be computed directly from their definitions. The stress intensity factors for a crack that is aligned with the Cartesian coordinate axes shown in Fig. 3 are defined as ⎧√ ⎫ { } lim ⎪ 2𝜋𝑟𝜎𝑦𝑦 |𝜃=0 ⎪ 𝐾𝐼 = 𝑟 → 0 ⎨√ (43) ⎬ 𝐾𝐼𝐼 ⎪ 2𝜋𝑟𝜏𝑥𝑦 |𝜃=0 ⎪ ⎩ ⎭
(35)
𝝈(𝜉, 𝜂) = 𝚿𝜎 (𝜂)𝜉 −𝚲n −𝐈 𝐜
𝐜
(s) (s) (s) 𝚿(s) 𝜎 (𝜂(𝜃)) =ℂ(−𝐁1 (𝜂(𝜃))𝚽u 𝚲n + 𝐁2 (𝜂(𝜃))𝚽u )
The complete displacement field of a point defined by the sector covered by a line element on the element is obtained by substituting Eq. (33) into Eq. (18) resulting in: 𝐮(𝜉, 𝜂) = 𝐍(𝜂)𝚽u 𝜉 −𝚲n 𝐜
−𝐈 (s)
where the singular stress [ ]T (s) (s) (s) 𝚿𝜎𝑥𝑥 (𝜂(𝜃)) 𝚿𝜎𝑦𝑦 (𝜂(𝜃)) 𝚿𝜏𝑥𝑦 (𝜂(𝜃)) is
where c are the integration constants that are obtained from the nodal displacements 𝐮b = 𝐮(𝜉 = 1) as: 𝐜 = 𝚽−1 u 𝐮b
(40)
(39) 12
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
Fig. 4. Plate with a double edge crack under tension: (a) geometry and boundary conditions and (b-c) domain discretized with polygonal elements for the pure virtual element method and the combined VEM and the SBFEM.
σ
a
a
L=2
H=1
σ
(a)
(b)
(c)
dimensional region is discretized with polygons using Polymesher [36], a MATLAB based meshing tool. The results from the present framework, especially mode-I and mode-II stress intensity factors are compared with empirical solutions. For all the examples, a state of plane strain condition is assumed.
of the crack length, a, to the width of the plate, H, is 𝑎∕𝐻 = 0.25. The material properties of the plate are: Young’s modulus, 𝐸 = 200 GPa and Poisson’s ratio, 𝜈 = 0.3. In this example, a state of plane strain condition is assumed. The empirical mode I SIF for the present configuration is given by:
5.1. Plate with double edge crack in tension
√ 𝐾Iref = 𝐶𝜎 𝜋𝑎
The plate with double edge crack subjected to a uniform tension at both ends as shown in Fig. 4 is considered. In the computations, the ratio
where C is a correction factor. For 𝑎∕𝑏 > 0.4, 𝑏 = 𝐻∕2, the correction factor is given by [37]: 13
(45)
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
Fig. 5. Plate with an oblique crack: geometry and boundary conditions.
σ2
x2
β σ1
B
2w
x1
σ1
2a A 2w
σ2
Fig. 6. Variation of the SIFs for a plate with an inclined crack with crack angle 𝛽.
14
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
Table 1 Plate with an edge crack in remote tension: convergence of mode I SIF. h
0.25 0.125 0.0625 0.03125 0.015625
Number
Number
CPU time (s)
√ 𝐾𝐼 ∕(𝜎 𝜋𝑎)
of Polygons
of Nodes
SBFEM
Present
SBFEM
Present
52 172 633 2258 8714
188 456 1422 4780 17,896
0.38 0.65 1.82 8.33 75.36
0.23 0.26 0.61 4.62 59.89
1.1702 1.1760 1.1766 1.1779 1.1774
1.1941 1.1892 1.1855 1.1828 1.1812
Table 2 Convergence of mode I SIF for a plate with an edge crack in remote tension: a comparison between VEM and the present approach. h
0.25 0.125 0.0625 0.03125 0.015625
𝐶 = 1.12 + 0.203
( ) ( )2 ( )3 𝑎 𝑎 𝑎 − 1.197 + 1.930 𝑏 𝑏 𝑏
Number of Polygons
CPU time(s)
√ 𝐾𝐼 ∕(𝜎 𝜋𝑎)
VEM
Present
VEM
Present
VEM
Present
58 178 639 2264 8720
52 172 633 2258 8714
0.20 0.32 1.16 5.62 55.66
0.23 0.26 0.61 4.62 59.89
1.2028 1.1922 1.1914 1.1909 1.1916
1.1941 1.1892 1.1855 1.1828 1.1812
√ 𝐾II = (𝜎2 − 𝜎1 ) sin 𝛽 cos 𝛽 𝜋𝑎
(46)
(47)
The material properties of the plate are: Young’s modulus, 𝐸 = 200 GPa and Poisson’s ratio, 𝜈 = 0.3. In this example, the plate is discretized with polygon meshes (containing 300 polygons). The variations of the mode I and mode II SIF with the crack orientation 𝛽 are presented in Fig. 6. It can be observed that the results from the present method agree very well with the reference solution.
The above factor corrects for an infinite plate with an accuracy√of 2%. For the chosen parameters, the reference normalized SIF is 𝐾𝐼 ∕(𝜎 𝜋𝑎) = 1.1635. The plate is discretized with a polygonal mesh. For the polygons containing the crack tip, we employ the SBFEM technique to capture the singularity. In these polygons, each edge is further discretized with 5 linear elements so that the angular variation of the SIF can be computed accurately [38] (c.f. Fig. 4c). For the elements that do not contain the crack tip, we employ the VEM with stabilization to compute the stiffness matrix. The convergence of the mode I SIF with mesh refinement is give in Table 1 for the present approach and for the conventional SBFEM. For comparison, same discretization density is used. In case of SBFEM, all the polygons are treated within the SBFEM framework, whilst in the present formulation, only the element containing the crack tip is treated using the SBFEM and the rest of the elements are treated withing VEM framework. The computational time (in sec) is also given. It is seen that the present method yields comparable and accurate results and it requires slightly lower computational effort. Table 2 presents a comparison of the CPU time (in sec) and the numerically estimated mode I SIF based on the present work and the VEM. In the case of the VEM, a conforming mesh is adopted (c.f. Fig. 4b). The SBFEM polygon used at the crack tip is subdivided into conforming polygons modelled by the VEM. An interaction integral is adopted to estimate the numerical mode I SIF. The size of the 𝐽 −domain is taken as half the crack length. It can be inferred that, although the VEM takes approximately same time as compared to the present work1 , the numerically computed SIF is not as accurate. This is due to the fact that the VEM basis are polynomial functions, that does not capture the singular stress fields accurately.
6. Conclusions A new framework to estimate the fracture parameters accurately is presented by combining the VEM and the SBFEM. The region in the vicinity of the crack tip is treated semi-analytically by the SBFEM and the rest of the domain is treated by the VEM framework. The salient features of the proposed framework are: (a) does not require structured elements with special treatment to capture singularity, such as quarter point elements and (b) does not require a priori knowledge of the asymptotic crack tip fields like the XFEM. From the systematic numerical study, it can be inferred that the proposed technique yields accurate results. Crack propagation and cracks in heterogeneous materials can be handled without much complexity, which are scope for future work. References [1] Areias P, Rabczuk T. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elem Anal Des 2017;132:27–41. [2] Areias P, Reinoso J, Camanho PP, de Sá JC, Rabczuk T. Effective 2D and 3D crack propagation with local mesh refinement and the screened poisson equation. Eng Fract Mech 2018;189:339–60. [3] Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods for material modelling. Model Simul Mater Sci. Eng 2009;17:043001–1–043001–24. [4] Melenk J, Babuška I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech. Eng 1996;139:289–314. [5] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 1999;45(5):601–20. [6] Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modelling. Int J Numer Methods Eng 2000;48(11):1549–70. [7] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Methods Eng 2004;61:2316–43. [8] Sukumar N, Malsch E. Recent advances in the construction of polygonal finite element interpolants. Arch Comput Methods Eng 2006;13(1):129–63. [9] Sukumar N, Tabarraei A. Conforming polygonal finite elements. Int J Numer Methods Eng 2004;61:2045–66. [10] Beirão Da Veiga L, Brezzi F, Cangiani A, Manzini G, Marini L, Russo A. The basic principles of virtual element methods. Math Models Methods Appl. Sci 2013;23:199–214. [11] Perumal L. A brief review on polygonal/polyhedral finite element methods. Math Probl Eng 2018;2018:1–22 Article ID 5792372.
5.2. Angled crack in an isotropic material In this example, a plate with an angled crack subjected to far field bi-axial stress field, 𝝈 (see Fig. 5) with 𝑎∕𝑤 = 0.1, 𝜎1 = 1 and 𝜎2 = 2 is considered. In this example, the mode I and the mode II SIFs, KI and KII , respectively, are obtained as a function of the crack angle 𝛽. For the loads shown, the analytical SIF for an infinite plate are given by Aliabadi et al. [39]: √ 𝐾I = (𝜎2 sin2 𝛽 + 𝜎1 cos2 𝛽) 𝜋𝑎 1 this can be attributed to the fact that in the present work, all elements are treated within VEM framework except the element containing the crack tip.
15
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
[12] L Beirão da Veiga, Manzini G. The mimetic finite difference method and the virtual element method for elliptic problems with arbitrary regularity. Tech. Rep. LA-UR-12-22977. Los Alamos National Laboratory; 2012. [13] Beirão Da Veiga L, Brezzi F, Marini L, Russo A. The Hitchhiker’s guide to the virtual element method. Math Models Methods Appl Sci 2014;24:1541. [14] Manzini G, Russo A, Sukumar N. New perspectives on polygonal and polyhedral finite element methods. Math Models Methods Appl Sci 2014;24:1665. [15] Nguyen-Thanh VM, Zhuang X, Nguyen-Xuan H, Rabczuk T, Wriggers P. A virtual element method for 2d linear elastic fracture analysis. Comput Methods Appl Mech Eng 2018;340:366–95. [16] hai Tang X, Wu S-C, Zheng C, hai Zhang J. A novel virtual node method for polygonal elements. Appl Math Mech 2009;30:1233–46. [17] Lee C, Kim H, Im S. Polyhedral elements by means of node/edge-based smoothed finite element method. Int J Numer Methods Eng 2016. doi:10.1002/nme.5449. [18] Francis A, AOrtiz-Bernardin, Bordas SPA, Natarajan S. Linear smoothed polygonal and polyhedral finite elements. Int J Numer Methods Eng 2017;109(9):1263– 1288. [19] Ooi E, Song C, Natarajan S. Construction of high-order complete scaled boundary shape functions over arbitrary polygons with bubble functions. Int J Numer Methods Eng 2016;108(9):1086–120. [20] Natarajan S, Ooi ET, Saputra A, Song C. A scaled boundary finite element formulation over arbitrary faceted star convex polyhedra. Engi Anal Bound Elem 2017;80:218–29. [21] Talebi H, Saputra A, Song C. Stress analysis of 3d complex geometries using the scaled boundary polyhedral finite elements. Comput Mech 2016;58(4):697–715. doi:10.1007/s00466-016-1312-0. [22] Cangiani A, Georgoulis EH, Houston P. Hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes. Math Models Methods Appl Sci 2014;24(10):2009–41. [23] Song C. Evaluation of power-logarithmic singularities, t-stresses and higher order terms of in-plane singular stress fields at cracks and multi-material corners. Eng Fract Mech 2005;72:1498–530. [24] Song C. Analysis of singular stress fields at multi-material corners under thermal loading. Int J Numer Methods Eng 2006;65:620–52. [25] Song C, Ooi ET, Natarajan S. A review of the scaled boundary finite element method for two-dimensional linear elastic fracture mechanics. Eng Fract Mech 2018;187:45–73.
[26] Huynh HD, Tran P, Zhuang X, Nguyen-Xuan H. An extended polygonal finite element method for large deformation fracture analysis. Eng Fract Mech 2019;209:344–68. [27] Huynh HD, Nguyen MN, Cusatis G, Tanaka S, Bui TQ. A polygonal XFEM with new numerical integration for linear elastic fracture mechanics. Eng Fract Mech 2019. doi:10.1016/j.engfracmech.2019.04.002. [28] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. Int J Numer Methods Eng 2016;108:1451–76. [29] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: a stable solution to varying horizons. Comput Methods Appl Mech Eng 2017;318:762–82. [30] Beirão Da Veiga L, Brezzi F, Marini L. Virtual elements for linear elasticity problems. SIAM J Numer Anal 2013;51:794–812. [31] Gain AL, Talischi C, Paulino GH. On the virtual element method for threedimensional linear elasticity problems on arbitrary polyhedra meshes. Comput Methods Appl Mech Eng 2014. doi:10.1016/j.cma.2014.05.005. [32] Wolf J, Song C. The scaled boundary finite-element method - a premier derivation. Comput Struct 2000;78:191–210. [33] Wolf J, Song C. The scaled boundary finite-element method - a fundamental solution-less boundary element method. Comput Methods Appl Mech Eng 2001;190:5551–68. [34] Deeks A, Wolf J. A virtual work derivation of the scaled boundary finite-element method for elastostatics. Comput Mech 2002;28:489–504. [35] Natarajan S, Ooi ET, Chiong I, Song C. Convergence and accuracy of displacement based finite element formulation over arbitrary polygons: laplace interpolants, strain smoothing and scaled boundary polygon formulation. Finite Elem Anal Des 2014;85:101–22. [36] Talischi C, Paulino GH, Pereira A, Menezes IF. Polytop: a matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes. Struct Multidisc Optim 2012;45:329–57. [37] Tada H, Paris P, Irwin G. The stress analysis of cracks. The Americal Society of Mechanical Engineers, New York; 2000. [38] Ooi ET, Song C, Tin-Loi F, Yang Z. Polygon scaled boundary finite elements for crack propagation modelling. Int J Numer Methods Eng 2012;91:319–42. [39] Aliabadi M, Rooke D, Cartwright D. Mixed-mode Bueckner weight functions using boundary element analysis. Int J Fract 1987;34:131–47.
16
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A combined virtual element method and the scaled boundary finite element method for linear elastic fracture mechanics Dibyendu Adak a, ALN Pramod a, Ean Tat Ooi b, Sundararajan Natarajan a,∗ a b
Integrated Modelling and Simulation Lab, Department of Mechanical Engineering, Indian Institute of Technology Madras, Chennai 600036, India School of Science, Engineering & Information Technology, Federation University Australia, Ballarat VIC3350, Australia
a r t i c l e
i n f o
Keywords: Stress intensity factor Virtual element method Scaled boundary finite element method Polygonal elements
a b s t r a c t In this paper, we propose a framework that combines the recently introduced virtual element method (VEM) and the scaled boundary finite element method (SBFEM) to evaluate the fracture parameters. The domain is discretized with arbitrary polygons and the element that contains the crack tip is treated within the framework of the SBFEM. This facilitates a semi-analytical treatment of the crack tip singularity allowing the fracture parameters are estimated directly from the definition. The VEM is employed for the rest of the domain. The salient feature of the VEM is that the terms in the stiffness matrix are computed without requiring higher order quadrature schemes. As both the methods satisfy partition of unity and the compatibility condition, the matrices are assembled as in the conventional FEM. The accuracy of the proposed formulation is demonstrated with two standard benchmark examples. The proposed VEM-SBFEM framework yields accurate results.
1. Introduction One of the computational challenges that the finite element method (FEM) faces is in modelling evolving discontinuities. Discontinuities such as cracks and notches require local mesh adaptation with special elements to treat the singular fields. Moreover, re-meshing [1,2] may be required as the discontinuities evolve in time. Another difficulty is that the conventional FEM is restricted to simplex elements, such as triangles and quadrilaterals. This poses a challenge for local mesh adaption. Thanks to the introduction of the partition of unity methods (PUM), such as the eXtended FEM (XFEM)/partition of unity FEM (PUFEM)/generalized FEM (GFEM) [3–6], element free Galerkin method (EFG) [7] and the polygonal FEM (PFEM) [8–10] which alleviates the meshing burden by detaching the geometry and the discretization. This allows the geometry to be represented independent of the underlying discretization. On the other hand, the PFEM allows the elements to take arbitrary number of sides and shapes. Both the approaches aim at relaxing the meshing burden and the topology constraint imposed by the conventional FEM. The PUM captures the discontinuity by supplementing the conventional FE basis with a priori known additional functions that capture the local behavior, whilst the PFEM allows the elements to take arbitrary shapes and can have arbitrary number of sides. Since the inception of the XFEM and the PFEM, these approaches have been applied to wider variety of problems. For a review on XFEM and PFEM, interested readers are referred to [3,11] and the references
∗
therein. The flexibility offered by the PFEM has led researchers to develop methods with polygonal discretizations, viz., virtual element method (VEM) [10,12–15], virtual node method [16], smoothed finite element method [17,18], scaled boundary finite element method (SBFEM) [19–21] and the discontinuous Galerkin method [22]. The VEM was recently developed in [10,12–15] and could be considered as a generalization of the FEM over arbitrary polytopes. Unlike the conventional FEM, the VEM does not require an explicit form for the basis functions to compute the discrete variational formulation. The VEM also alleviates the numerical integration difficulty encountered in the conventional polygonal FEM. As the method does not require the knowledge of the shape functions in explicit form, the implementation is computationally less intensive. However, since the VEM space contains only polynomial basis, the VEM cannot accurately capture the asymptotic crack tip fields in the vicinity of the crack tip. On another front, the scaled boundary finite element method (SBFEM) introduced by Song and Wolf has emerged as an attractive alternate to model problems with singularities [23,24]. It combines the best features of the FEM and the boundary element method (BEM). In [25] a recent review on the application of the SBFEM to linear elastic fracture mechanics can be found. Recently, Hyun and co-workers [26,27], proposed extended PFEM (XPFEM) by combining the XFEM with the PFEM. Within this formulation, the background discretization consists of arbitrary polygons and the XFEM allows for capturing the discontinuities without requiring a conforming mesh. However, the method requires a priori knowledge of the enrichment functions and special numerical integration schemes
Corresponding author. E-mail address: [email protected] (S. Natarajan).
https://doi.org/10.1016/j.enganabound.2019.12.008 Received 15 April 2019; Received in revised form 15 October 2019; Accepted 16 December 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
The corresponding weak form is: find 𝐮 ∈ 𝒰 such that
Ωh
Γt
∀𝐯 ∈ 𝒱 , 𝑎 (𝐮, 𝐯) =
Γc
Fig. 1. General elastic body with an internal discontinuity, Neumann and Dirichlet boundary conditions.
𝝈(𝐮) ∶ 𝜺(𝐯) d𝑉 , 𝓁(𝐯) =
∫Ω
𝐛 ⋅ 𝐯 d𝑉 +
∫Γ𝑡
𝐭̂ ⋅ 𝐯 d𝑆,
(3b)
where the space (Ω) includes the linear displacement fields. The domain is partitioned into elements Ωh . Using shape functions 𝜙a that span at least the linear space, we substitute vector-valued trial and test func∑ ∑ tions 𝐮ℎ = 𝑎 𝜙𝑎 𝐮𝑎 and 𝐯ℎ = 𝑏 𝜙𝑏 𝐯𝑏 , respectively, into Eq. (3) and apply a standard Galerkin procedure to obtain the discrete weak form: find 𝐮ℎ ∈ 𝒰 ℎ such that
have to be developed to integrate the enrichment functions. A different numerical formulation that also shows promise when applied to model fracture related problems is peridynamics. Peridynamics does not require an explicit description of the crack topology, an advantage that has attracted considerable amount of interest among researchers. However, peridynamics has a problem of spurious wave reflections and ghost forces between the particles [28,29]. The main objective of the paper is to combine the VEM with the SBFEM. The element that contains the crack tip is treated within the framework of the SBFEM. This facilitates a semi-analytical treatment of the crack tip singularity and the fracture parameters are estimated directly from the definition. The VEM is employed for the rest of the domain. The salient features of the proposed framework are: (a) does not require a priori knowledge of the enrichment functions, (b) the stress intensity factors are computed directly from the definition without requiring ‘any’ additional post processing, such as J-integral and (c) the SBFEM which requires a solution of an eigenvalue problem is restricted to a smaller region, thus increasing the computational efficiency. The rest of the manuscript is organized as follows: Section 2 presents the governing equation and the corresponding weak form. An overview of the VEM and the SBFEM is presented in Sections 3 and 4, respectively. Section 4 also presents the procedure to estimate the stress intensity factors directly from the definition. The accuracy of the estimated stress intensity factors for two standard problems in linear elastic fracture mechanics from the proposed framework are compared with available results in the literature in Section 5, followed by concluding remarks in the last section.
∀𝐯ℎ ∈ 𝒱 ℎ
𝑎(𝐮ℎ , 𝐯ℎ ) = 𝓁(𝐯ℎ ),
(5)
which leads to the following system of linear equations:
𝐊=
∑
𝐊𝐮= 𝐟 , ∑ 𝐊ℎ =
ℎ
𝐟=
∑
ℎ
𝐟 ℎ=
ℎ
∑
(6a) ∫Ωℎ (
𝐁T ℂ𝐁 d𝑉 ,
∫Ωℎ
ℎ
(6b) )
𝐍T 𝐛 d𝑉 +
∫Γℎ
𝐍T 𝐭̂ d𝑆 ,
(6c)
𝑡
where K is the assembled stiffness matrix, f the assembled nodal force vector, u the assembled vector of nodal displacements, N is the matrix of shape functions, ℂ is the constitutive matrix for an isotropic linear elastic material, and 𝐁 = 𝛁𝐍 is the strain-displacement matrix that is computed using the derivatives of the shape functions. In this study, the domain is assumed to be partitioned into non-overlapping elements with arbitrary number of sides. The main challenge is in computing the terms in the bilinear and linear form, as the shape functions over such arbitrary polygons are not simple polynomials [8,9]. 3. Basics of the virtual element method In this section, we present the virtual element formulation for the elastic deformation problem. Only important equations are presented here and for more detailed discussion and implementation, interested readers are referred to [30,31] and the references therein. The continuous variational formulation of the elastic deformation problem is given by Eq. (5). Since the virtual element space does not require explicit expression of the basis function, two projection operators 𝐸0 and 𝐸∇ are introduced to compute the bilinear and the linear form. Employing these projection operators, we project the virtual basis functions onto a computable polynomial subspace. The tensor valued L2 projection operator is defined as follows:
2. Governing equations and weak form In this section, the governing partial differential equations and the corresponding weak form for linear homogeneous elastic body containing a strong discontinuity (i.e., crack) is given. Consider a linear elastic body occupying Ω ⊂ IR2 bounded by a surface Γ, with an outward normal n and a discontinuity Γc as shown in Fig. 1. The boundary is assumed to admit decompositions: Γ = Γ𝑢 ∪ Γ𝑡 ∪ Γ𝑐 , where Γu is the Dirichlet boundary, Γt is the Neumann boundary and Γc , is the traction free surface representing the discontinuity. In the absence of interia and body forces, the boundary value problem for linear elasto-statics is: find 𝐮 ∶ Ω → ℝ2 , such that: in Ω on Γ𝑢 on Γ𝑡 , on Γ𝑐 ,
(3a)
where 𝜺 is the small strain tensor, and 𝒰 and 𝒱 are the displacement trial and test spaces: { } 𝒰 ∶= 𝐮(𝐱) ∈ [𝐶 0 (Ω)]𝑑 ∶ 𝐮 ∈ [(Ω)]𝑑 ⊆ [𝐻 1 (Ω)]𝑑 , 𝐮 = 𝐮̂ on Γ𝑢 , { } 𝒱 ∶= 𝐯(𝐱) ∈ [𝐶 0 (Ω)]𝑑 ∶ 𝐯 ∈ [(Ω)]𝑑 ⊆ [𝐻 1 (Ω)]𝑑 , 𝐯 = 𝟎 on Γ𝑢 ,
Γu
⎧ ⎪𝛁 ⋅ 𝝈 = 𝟎 ⎪𝐮 = 𝐮̂ ⎨𝝈 ⋅ 𝐧 = 𝐭̂ ⎪ ⎪𝝈 ⋅ 𝐧 = 𝟎 ⎩
∫Ω
𝑎(𝐮, 𝐯)= 𝓁(𝐯)
𝐸0 𝐐 ∶=
1 𝐐 |𝐸| ∫𝐸
∀ 𝐸 ∈ Ωℎ .
(7)
The global projection operator is defined as 0 (𝐐)|𝐸 = 𝐸0 (𝐐). In addition, we introduce a tensor valued elliptic type projection operator ∇ , defined on H1 (E)d as follows: for any 𝐯 ∈ ℎ we have 𝐸∇ 𝐯 ∈ 𝐏1 (𝐸) and it satisfies
(1)
⎧∇(𝐸∇ (𝐯)) ⎪ ⎨∑ ∇ 𝐸 (𝑙) ⎪ ⎩𝑙∈𝜕𝐸
where 𝛁 is the gradient operator, 𝝈 is the Cauchy stress tensor, n is the unit outward normal, u is the displacement field and 𝐮̂ and 𝐭̂ are the applied displacement and traction on the Dirichlet and Neumann boundary. For small displacements and strains, the strain displacement relation is given by: ) 1( 𝜺= 𝛁𝐮 + 𝛁𝐮T (2) 2
= 𝐸0 (∇𝐯|𝐸 ) =
∑
𝐯(𝑙)
(8)
𝑙∈𝜕𝐸
The local virtual element space is defined as ℎ,𝐸 ∶= {𝑣 ∈ [𝐻 1 (𝐸) ∩ 𝐶 0 (𝐸)]2 ∶ Δ𝐯 = 0 in 𝐸, 𝐯|𝐹 ∈ 𝐏1 (𝐹 ) ∀ 𝐹 ∈ 𝜕𝐸}. (9) 10
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
η
and the corresponding global virtual element space is defined as ℎ ∶= {𝐯 ∈ 𝐻10 (Ω)2 ∶ 𝐯|𝐸 ∈ ℎ,𝐸 }.
(10)
The expression for the stiffness matrix in the case of two dimensional linear elasticity can be written as: 𝐊𝐸 ℎ =
|𝐸|𝐖𝐶 ℂ𝐖T𝐶 ⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟
Consist ency t erm=𝐊const
+ (𝐈 − 𝐏𝑝 )T 𝐒𝐸 (𝐈 − 𝐏𝑝 ) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟
Line element
(x2,y2), η = 1
(x1,y1), η = −1
(11)
ξ=1
Stabilit y t erm=𝐊stab
S
where 𝐏𝑝 = 𝐏𝑅 + 𝐏𝐶
(12)
and 𝐏𝑅 = 𝐍𝑅 𝐖T𝑅 𝐏𝐶 = 𝐍𝐶 𝐖T𝐶 The NR and NC are expressed as: [ ] 1 0 (𝐱𝐼 − 𝐱 )2 (𝐍𝑅 )𝐼 = 0 1 −(𝐱𝐼 − 𝐱)1 [ ] (𝐱𝐼 − 𝐱 )1 0 (𝐱𝐼 − 𝐱 )2 (𝐍𝐶 )2 = 0 (𝐱𝐼 − 𝐱 )2 (𝐱𝐼 − 𝐱 )1
(13)
Radial Lines A
B
Scaling Centre
(14)
Fig. 2. Cracked polygon representation by the scaled boundary finite element method.
where xI is the coordinate of the node and 𝐱 is the polyhedron centroid. The block rows of WR and WC are expressed as: [ ] 1∕𝑛 0 (𝐪𝐼 )2 (𝐖𝑅 )𝐼 = 0 1∕𝑛 −(𝐪𝐼 )1 [ ] 2(𝐪𝐼 )1 0 (𝐪𝐼 )2 (𝐖𝐶 )𝐼 = (15) 0 2(𝐪𝐼 )2 (𝐪𝐼 )1
where N(𝜂) is the shape function matrix of the finite elements discretising the polygon boundary. The standard 1D Gauss–Lobatto shape functions or Lagrange shape functions can be used. In this study, we employ Lagrange shape functions. The displacements at a point in a polygon is approximated by: 𝐮(𝜉, 𝜂) = 𝐍(𝜂)𝐮(𝜉)
where n is the number of vertices and 𝐼 = 1, … , 𝑁 with N the number of nodes in the element and the subscript indicates the component of the associated vector and 𝐪𝐼 =
1 𝑁𝐼 𝐧 dΓ 2|𝐸| ∫
where u(𝜉) are the radial displacement functions. Substituting Eq. (18) in the definition of strain-displacement relations, the strains 𝜺(𝜉, 𝜂) in the scaled boundary coordinates are expressed as:
(16)
𝜕𝐸
𝜺(𝜉, 𝜂) = 𝐋𝐮(𝜉, 𝜂)
where 𝐒𝐸 = 𝛼𝐈 and 𝛼 = 𝛼 ∗ trace(|𝐸|𝐖𝐶 ℂ𝐖T𝐶 ) is a scaling coefficient and n is the unit outward normal vector. It can be seen that the computation of the stiffness matrix involves computing the matrices NR , NC , WR and WC . The calculation of the matrices WR and WC involves computing the surface integral of the basis functions.
(19)
where L is a linear operator matrix formulated in the scaled boundary coordinates as 𝐋 = 𝐛1 (𝜂)
4. Overview of the scaled boundary finite element method
𝜕 1 𝜕 + 𝐛 (𝜂) 𝜕𝜉 𝜉 2 𝜕𝜂
(20)
with
The VEM discussed in the last section can be applied to problems with strong discontinuity and singularity. However, to accurately capture the asymptotic fields at the crack tip, a very fine mesh in combination with singular elements at the crack tip is usually required. This poses additional difficulties when the crack evolves. Another possibility is to enrich the approximation space with functions that can capture the discontinuity and singularity [3,4]. In the literature, the latter method is referred to as the PUFEM. In this study, we propose to couple the VEM with the SBFEM to model problems with strong discontinuity and singularities. The SBFEM is a semi-analytical method and relies on defining a ‘scaling centre’ from which the entire boundary is visible. This is similar to the concept of ‘star convexity’. The boundary is divided into conventional linear finite elements, whilst the solution is sought analytically in the radial direction [32]. Moreover, by exploiting the special characteristics of the scaling centre, the stress intensity factors can be computed directly from their definition. When modelling a crack/notched surface the scaling centre is placed at the crack tip. The straight crack/notch edges are formed by scaling the nodes A and B on the boundary and the crack surfaces are not discretized (see Fig. 2). Displacement approximation The geometry of the element described by the coordinates on the boundary xb (𝜂) is expressed as: 𝐱𝑏 (𝜂) = 𝐍(𝜂)𝐱𝑏
(18)
𝐛1 (𝜂) =
⎡ 𝑦 (𝜂) 1 ⎢ 𝜂 ,𝜂 0 |𝐉(𝜂)| ⎢ ⎣−𝑥𝜂 (𝜂),𝜂
𝐛2 (𝜂) =
⎡−𝑦 (𝜂) 1 ⎢ 𝜂 0 |𝐉(𝜂)| ⎢ ⎣ 𝑥𝜂 (𝜂)
0 ⎤ −𝑥𝜂 (𝜂),𝜂 ⎥ ⎥ 𝑦𝜂 (𝜂),𝜂 ⎦ 0 ⎤ 𝑥𝜂 (𝜂) ⎥ ⎥ −𝑦𝜂 (𝜂)⎦
(21)
Using Eq. (19) and Hooke’s law 𝝈 = ℂ𝜺, the stresses 𝝈(𝜉, 𝜂) is expressed as ( ) 𝝈(𝜉, 𝜂) = ℂ 𝐁1 (𝜂)𝐮(𝜉),𝜉 + 𝜉 −1 𝐁2 (𝜂)𝐮(𝜉)
(22)
where ℂ is the material constitutive matrix. B1 (𝜂) and B2 (𝜂) are the strain-displacement matrices in the SBFEM and is defined as 𝐁1 (𝜂) = 𝐛1 (𝜂)𝐍(𝜂) 𝐁2 (𝜂) = 𝐛2 (𝜂)𝐍(𝜂),𝜂
(23)
By following the procedure outlined in [33,34], the following ODE is obtained: 𝐄0 𝜉 2 𝐮(𝜉),𝜉𝜉 + (𝐄0 + 𝐄T1 − 𝐄1 )𝜉𝐮(𝜉),𝜉 − 𝐄2 𝐮(𝜉) = 0 where E0 , E1 and E2 are the coefficient matrices given by:
(17) 11
(24)
D. Adak, A. Pramod and E.T. Ooi et al.
𝐄0 = 𝐄1 = 𝐄2 =
∫𝜂 ∫𝜂 ∫𝜂
Engineering Analysis with Boundary Elements 113 (2020) 9–16
𝐁1 (𝜂)T ℂ𝐁1 (𝜂)|𝐉(𝜂)|𝑑𝜂, 𝐁2 (𝜂)T ℂ𝐁1 (𝜂)|𝐉(𝜂)|𝑑𝜂, 𝐁2 (𝜂)T ℂ𝐁2 (𝜂)|𝐉(𝜂)|𝑑𝜂.
(25)
y
x
and the determinant of the Jacobian matrix is: |𝐉(𝜂)| = 𝑥𝑏 (𝜂)𝑦𝑏 (𝜂),𝜂 − 𝑦𝑏 (𝜂)𝑥𝑏 (𝜂),𝜂
O
(26)
The coefficient matrices are evaluated element-by-element on the polygon boundary and assembled over a polygon. This process is similar to the standard FE procedure of assemblage. Eq. (24) is a homogeneous second-order ordinary differential equation. Its solution is obtained by introducing the variable 𝝌(𝜉) { } 𝐮(𝜉) 𝝌(𝜉) = (27) 𝐪(𝜉)
ro
Fig. 3. A cracked domain modelled by the SBFEM and the definition of local coordinate system, where the ‘black’ dots represent the nodes.
where q(𝜉) is the internal load vector 𝐪(𝜉) = 𝐄0 𝜉𝐮(𝜉),𝜉 + 𝐄T1 𝐮(𝜉)
Lo
(28)
The boundary nodal forces are related to the displacement functions by:
From Eq. (39), the stiffness matrix K can be identified to be given by the expression
𝐟 = 𝐪(𝜉 = 1) = (𝐄0 𝜉𝐮(𝜉),𝜉 + 𝐄T1 𝐮(𝜉))|𝜉=1
𝐊 = 𝚽q 𝚽−1 u
(29)
Remark 4.1. The stiffness computed by employing the SBFEM is positive definite and symmetric. Hence, the stiffness matrix can be assembled in the conventional FEM approach. A simple Matlab ○R function is given in [35] to compute the stiffness matrix using the SBFEM.
This allows Eq. (24) to be transformed into a first order ordinary differential equation with twice the number of unknowns in an element as: 𝜉𝝌(𝜉),𝜉 = −𝐙𝝌(𝜉) where Z is the Hamiltonian matrix [ ] 𝐄−1 𝐄T1 −𝐄−1 0 0 𝐙= 𝐄1 𝐄−1 𝐄T1 − 𝐄2 −𝐄1 𝐄−1 0 0
(30)
Calculation of the stress intensity factors A unique feature of the SBFEM is that stress singularities, if present, are analytically represented in the radial displacement functions u(𝜉). When a crack is modelled by a polygon with its scaling centre chosen at the crack tip in Fig. 2, some of (s) the eigenvalues 𝚲(s) n ⊂ 𝚲n satisfy −1 < 𝚲n < 0. These eigenvalues lead to singular stresses at the crack tip i.e. when 𝜉 → 0. Using 𝚲(s) n , the singular stress field 𝝈 (s) (𝜉, 𝜂) can be defined as:
(31)
An eigenvalue decomposition of Z is performed. The blocks of eigenvalues and transformation matrices necessary are: [ ] [ ] 𝚽 𝚽u 𝐙 u = 𝚲 (32) 𝚽q 𝚽q n ( ) In Eq. (32), 𝚲n = diag 𝜆1 , 𝜆2 , … , 𝜆𝑛 contains the eigenvalues with negative real part. 𝚽u and 𝚽q are the corresponding transformation matrices of 𝚲n . They represent the modal displacements and forces, respectively. The general solution of Eq. (30) is given by: 𝐮(𝜉) = 𝚽u 𝜉 −𝚲n 𝐜
(33)
𝐪(𝜉) = 𝚽q 𝜉 −𝚲n 𝐜
(34)
(s)
− 𝚲n 𝝈 (s) (𝜉, 𝜂) =𝚿(s) 𝜎 (𝜂(𝜃))𝜉
𝚽(s) u
Taking the derivative of u(𝜉) with respect to 𝜉 and substituting into Eq. (22), the stress field 𝝈(𝜉, 𝜂) in the scaled boundary coordinates can be expressed as:
where the stress mode 𝚿𝜎 (𝜂) is defined as: ( ) 𝚿𝜎 (𝜂) = ℂ −𝐁1 (𝜂)𝚽u 𝚲n + 𝐁2 (𝜂)𝚽u
(38)
{
𝐾𝐼 𝐾𝐼𝐼
=
(42)
c(s)
}
⎧ (s) (s) ⎫ √ ⎪𝚿𝜎𝑦𝑦 (𝜂(𝜃 = 0))𝐜 ⎪ = 2𝜋𝐿𝑜 ⎨ ⎬ (s) ⎪𝚿𝜏𝑥𝑦 (𝜂(𝜃 = 0))𝐜(s) ⎪ ⎩ ⎭
(44)
5. Numerical examples In this section, the accuracy and the convergence properties of the proposed framework is demonstrated by solving two standard benchmark problems within the framework of linear elastic fracture mechanics. The discretization is based on centroid Voronoi tessellation. The two
The stiffness matrix of an element is obtained by first substituting Eq. (35) into Eq. (34) at 𝜉 = 1. This results in: 𝐟 = 𝚽q 𝚽−1 u 𝐮b
mode
Substituting the stress components in Eq. (41) at angle 𝜃 = 0 into Eq. (43) and using the relation 𝜉 = 𝑟∕𝐿𝑜 at 𝜃 = 0, the stress intensity factors are
(36)
(37)
(41) 𝚿(s) 𝜎 (𝜂(𝜃))
In Eq. (42) ⊂ 𝚽u and ⊂ c, contain the displacement modes and the integration constants corresponding to 𝚲(s) n . It can be discerned from Eq. (41) that 𝚲(s) n leads to singular stresses at the crack tip. This enables the stress intensity factors to be computed directly from their definitions. The stress intensity factors for a crack that is aligned with the Cartesian coordinate axes shown in Fig. 3 are defined as ⎧√ ⎫ { } lim ⎪ 2𝜋𝑟𝜎𝑦𝑦 |𝜃=0 ⎪ 𝐾𝐼 = 𝑟 → 0 ⎨√ (43) ⎬ 𝐾𝐼𝐼 ⎪ 2𝜋𝑟𝜏𝑥𝑦 |𝜃=0 ⎪ ⎩ ⎭
(35)
𝝈(𝜉, 𝜂) = 𝚿𝜎 (𝜂)𝜉 −𝚲n −𝐈 𝐜
𝐜
(s) (s) (s) 𝚿(s) 𝜎 (𝜂(𝜃)) =ℂ(−𝐁1 (𝜂(𝜃))𝚽u 𝚲n + 𝐁2 (𝜂(𝜃))𝚽u )
The complete displacement field of a point defined by the sector covered by a line element on the element is obtained by substituting Eq. (33) into Eq. (18) resulting in: 𝐮(𝜉, 𝜂) = 𝐍(𝜂)𝚽u 𝜉 −𝚲n 𝐜
−𝐈 (s)
where the singular stress [ ]T (s) (s) (s) 𝚿𝜎𝑥𝑥 (𝜂(𝜃)) 𝚿𝜎𝑦𝑦 (𝜂(𝜃)) 𝚿𝜏𝑥𝑦 (𝜂(𝜃)) is
where c are the integration constants that are obtained from the nodal displacements 𝐮b = 𝐮(𝜉 = 1) as: 𝐜 = 𝚽−1 u 𝐮b
(40)
(39) 12
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
Fig. 4. Plate with a double edge crack under tension: (a) geometry and boundary conditions and (b-c) domain discretized with polygonal elements for the pure virtual element method and the combined VEM and the SBFEM.
σ
a
a
L=2
H=1
σ
(a)
(b)
(c)
dimensional region is discretized with polygons using Polymesher [36], a MATLAB based meshing tool. The results from the present framework, especially mode-I and mode-II stress intensity factors are compared with empirical solutions. For all the examples, a state of plane strain condition is assumed.
of the crack length, a, to the width of the plate, H, is 𝑎∕𝐻 = 0.25. The material properties of the plate are: Young’s modulus, 𝐸 = 200 GPa and Poisson’s ratio, 𝜈 = 0.3. In this example, a state of plane strain condition is assumed. The empirical mode I SIF for the present configuration is given by:
5.1. Plate with double edge crack in tension
√ 𝐾Iref = 𝐶𝜎 𝜋𝑎
The plate with double edge crack subjected to a uniform tension at both ends as shown in Fig. 4 is considered. In the computations, the ratio
where C is a correction factor. For 𝑎∕𝑏 > 0.4, 𝑏 = 𝐻∕2, the correction factor is given by [37]: 13
(45)
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
Fig. 5. Plate with an oblique crack: geometry and boundary conditions.
σ2
x2
β σ1
B
2w
x1
σ1
2a A 2w
σ2
Fig. 6. Variation of the SIFs for a plate with an inclined crack with crack angle 𝛽.
14
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
Table 1 Plate with an edge crack in remote tension: convergence of mode I SIF. h
0.25 0.125 0.0625 0.03125 0.015625
Number
Number
CPU time (s)
√ 𝐾𝐼 ∕(𝜎 𝜋𝑎)
of Polygons
of Nodes
SBFEM
Present
SBFEM
Present
52 172 633 2258 8714
188 456 1422 4780 17,896
0.38 0.65 1.82 8.33 75.36
0.23 0.26 0.61 4.62 59.89
1.1702 1.1760 1.1766 1.1779 1.1774
1.1941 1.1892 1.1855 1.1828 1.1812
Table 2 Convergence of mode I SIF for a plate with an edge crack in remote tension: a comparison between VEM and the present approach. h
0.25 0.125 0.0625 0.03125 0.015625
𝐶 = 1.12 + 0.203
( ) ( )2 ( )3 𝑎 𝑎 𝑎 − 1.197 + 1.930 𝑏 𝑏 𝑏
Number of Polygons
CPU time(s)
√ 𝐾𝐼 ∕(𝜎 𝜋𝑎)
VEM
Present
VEM
Present
VEM
Present
58 178 639 2264 8720
52 172 633 2258 8714
0.20 0.32 1.16 5.62 55.66
0.23 0.26 0.61 4.62 59.89
1.2028 1.1922 1.1914 1.1909 1.1916
1.1941 1.1892 1.1855 1.1828 1.1812
√ 𝐾II = (𝜎2 − 𝜎1 ) sin 𝛽 cos 𝛽 𝜋𝑎
(46)
(47)
The material properties of the plate are: Young’s modulus, 𝐸 = 200 GPa and Poisson’s ratio, 𝜈 = 0.3. In this example, the plate is discretized with polygon meshes (containing 300 polygons). The variations of the mode I and mode II SIF with the crack orientation 𝛽 are presented in Fig. 6. It can be observed that the results from the present method agree very well with the reference solution.
The above factor corrects for an infinite plate with an accuracy√of 2%. For the chosen parameters, the reference normalized SIF is 𝐾𝐼 ∕(𝜎 𝜋𝑎) = 1.1635. The plate is discretized with a polygonal mesh. For the polygons containing the crack tip, we employ the SBFEM technique to capture the singularity. In these polygons, each edge is further discretized with 5 linear elements so that the angular variation of the SIF can be computed accurately [38] (c.f. Fig. 4c). For the elements that do not contain the crack tip, we employ the VEM with stabilization to compute the stiffness matrix. The convergence of the mode I SIF with mesh refinement is give in Table 1 for the present approach and for the conventional SBFEM. For comparison, same discretization density is used. In case of SBFEM, all the polygons are treated within the SBFEM framework, whilst in the present formulation, only the element containing the crack tip is treated using the SBFEM and the rest of the elements are treated withing VEM framework. The computational time (in sec) is also given. It is seen that the present method yields comparable and accurate results and it requires slightly lower computational effort. Table 2 presents a comparison of the CPU time (in sec) and the numerically estimated mode I SIF based on the present work and the VEM. In the case of the VEM, a conforming mesh is adopted (c.f. Fig. 4b). The SBFEM polygon used at the crack tip is subdivided into conforming polygons modelled by the VEM. An interaction integral is adopted to estimate the numerical mode I SIF. The size of the 𝐽 −domain is taken as half the crack length. It can be inferred that, although the VEM takes approximately same time as compared to the present work1 , the numerically computed SIF is not as accurate. This is due to the fact that the VEM basis are polynomial functions, that does not capture the singular stress fields accurately.
6. Conclusions A new framework to estimate the fracture parameters accurately is presented by combining the VEM and the SBFEM. The region in the vicinity of the crack tip is treated semi-analytically by the SBFEM and the rest of the domain is treated by the VEM framework. The salient features of the proposed framework are: (a) does not require structured elements with special treatment to capture singularity, such as quarter point elements and (b) does not require a priori knowledge of the asymptotic crack tip fields like the XFEM. From the systematic numerical study, it can be inferred that the proposed technique yields accurate results. Crack propagation and cracks in heterogeneous materials can be handled without much complexity, which are scope for future work. References [1] Areias P, Rabczuk T. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elem Anal Des 2017;132:27–41. [2] Areias P, Reinoso J, Camanho PP, de Sá JC, Rabczuk T. Effective 2D and 3D crack propagation with local mesh refinement and the screened poisson equation. Eng Fract Mech 2018;189:339–60. [3] Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods for material modelling. Model Simul Mater Sci. Eng 2009;17:043001–1–043001–24. [4] Melenk J, Babuška I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech. Eng 1996;139:289–314. [5] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 1999;45(5):601–20. [6] Sukumar N, Moës N, Moran B, Belytschko T. Extended finite element method for three-dimensional crack modelling. Int J Numer Methods Eng 2000;48(11):1549–70. [7] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Methods Eng 2004;61:2316–43. [8] Sukumar N, Malsch E. Recent advances in the construction of polygonal finite element interpolants. Arch Comput Methods Eng 2006;13(1):129–63. [9] Sukumar N, Tabarraei A. Conforming polygonal finite elements. Int J Numer Methods Eng 2004;61:2045–66. [10] Beirão Da Veiga L, Brezzi F, Cangiani A, Manzini G, Marini L, Russo A. The basic principles of virtual element methods. Math Models Methods Appl. Sci 2013;23:199–214. [11] Perumal L. A brief review on polygonal/polyhedral finite element methods. Math Probl Eng 2018;2018:1–22 Article ID 5792372.
5.2. Angled crack in an isotropic material In this example, a plate with an angled crack subjected to far field bi-axial stress field, 𝝈 (see Fig. 5) with 𝑎∕𝑤 = 0.1, 𝜎1 = 1 and 𝜎2 = 2 is considered. In this example, the mode I and the mode II SIFs, KI and KII , respectively, are obtained as a function of the crack angle 𝛽. For the loads shown, the analytical SIF for an infinite plate are given by Aliabadi et al. [39]: √ 𝐾I = (𝜎2 sin2 𝛽 + 𝜎1 cos2 𝛽) 𝜋𝑎 1 this can be attributed to the fact that in the present work, all elements are treated within VEM framework except the element containing the crack tip.
15
D. Adak, A. Pramod and E.T. Ooi et al.
Engineering Analysis with Boundary Elements 113 (2020) 9–16
[12] L Beirão da Veiga, Manzini G. The mimetic finite difference method and the virtual element method for elliptic problems with arbitrary regularity. Tech. Rep. LA-UR-12-22977. Los Alamos National Laboratory; 2012. [13] Beirão Da Veiga L, Brezzi F, Marini L, Russo A. The Hitchhiker’s guide to the virtual element method. Math Models Methods Appl Sci 2014;24:1541. [14] Manzini G, Russo A, Sukumar N. New perspectives on polygonal and polyhedral finite element methods. Math Models Methods Appl Sci 2014;24:1665. [15] Nguyen-Thanh VM, Zhuang X, Nguyen-Xuan H, Rabczuk T, Wriggers P. A virtual element method for 2d linear elastic fracture analysis. Comput Methods Appl Mech Eng 2018;340:366–95. [16] hai Tang X, Wu S-C, Zheng C, hai Zhang J. A novel virtual node method for polygonal elements. Appl Math Mech 2009;30:1233–46. [17] Lee C, Kim H, Im S. Polyhedral elements by means of node/edge-based smoothed finite element method. Int J Numer Methods Eng 2016. doi:10.1002/nme.5449. [18] Francis A, AOrtiz-Bernardin, Bordas SPA, Natarajan S. Linear smoothed polygonal and polyhedral finite elements. Int J Numer Methods Eng 2017;109(9):1263– 1288. [19] Ooi E, Song C, Natarajan S. Construction of high-order complete scaled boundary shape functions over arbitrary polygons with bubble functions. Int J Numer Methods Eng 2016;108(9):1086–120. [20] Natarajan S, Ooi ET, Saputra A, Song C. A scaled boundary finite element formulation over arbitrary faceted star convex polyhedra. Engi Anal Bound Elem 2017;80:218–29. [21] Talebi H, Saputra A, Song C. Stress analysis of 3d complex geometries using the scaled boundary polyhedral finite elements. Comput Mech 2016;58(4):697–715. doi:10.1007/s00466-016-1312-0. [22] Cangiani A, Georgoulis EH, Houston P. Hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes. Math Models Methods Appl Sci 2014;24(10):2009–41. [23] Song C. Evaluation of power-logarithmic singularities, t-stresses and higher order terms of in-plane singular stress fields at cracks and multi-material corners. Eng Fract Mech 2005;72:1498–530. [24] Song C. Analysis of singular stress fields at multi-material corners under thermal loading. Int J Numer Methods Eng 2006;65:620–52. [25] Song C, Ooi ET, Natarajan S. A review of the scaled boundary finite element method for two-dimensional linear elastic fracture mechanics. Eng Fract Mech 2018;187:45–73.
[26] Huynh HD, Tran P, Zhuang X, Nguyen-Xuan H. An extended polygonal finite element method for large deformation fracture analysis. Eng Fract Mech 2019;209:344–68. [27] Huynh HD, Nguyen MN, Cusatis G, Tanaka S, Bui TQ. A polygonal XFEM with new numerical integration for linear elastic fracture mechanics. Eng Fract Mech 2019. doi:10.1016/j.engfracmech.2019.04.002. [28] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. Int J Numer Methods Eng 2016;108:1451–76. [29] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: a stable solution to varying horizons. Comput Methods Appl Mech Eng 2017;318:762–82. [30] Beirão Da Veiga L, Brezzi F, Marini L. Virtual elements for linear elasticity problems. SIAM J Numer Anal 2013;51:794–812. [31] Gain AL, Talischi C, Paulino GH. On the virtual element method for threedimensional linear elasticity problems on arbitrary polyhedra meshes. Comput Methods Appl Mech Eng 2014. doi:10.1016/j.cma.2014.05.005. [32] Wolf J, Song C. The scaled boundary finite-element method - a premier derivation. Comput Struct 2000;78:191–210. [33] Wolf J, Song C. The scaled boundary finite-element method - a fundamental solution-less boundary element method. Comput Methods Appl Mech Eng 2001;190:5551–68. [34] Deeks A, Wolf J. A virtual work derivation of the scaled boundary finite-element method for elastostatics. Comput Mech 2002;28:489–504. [35] Natarajan S, Ooi ET, Chiong I, Song C. Convergence and accuracy of displacement based finite element formulation over arbitrary polygons: laplace interpolants, strain smoothing and scaled boundary polygon formulation. Finite Elem Anal Des 2014;85:101–22. [36] Talischi C, Paulino GH, Pereira A, Menezes IF. Polytop: a matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes. Struct Multidisc Optim 2012;45:329–57. [37] Tada H, Paris P, Irwin G. The stress analysis of cracks. The Americal Society of Mechanical Engineers, New York; 2000. [38] Ooi ET, Song C, Tin-Loi F, Yang Z. Polygon scaled boundary finite elements for crack propagation modelling. Int J Numer Methods Eng 2012;91:319–42. [39] Aliabadi M, Rooke D, Cartwright D. Mixed-mode Bueckner weight functions using boundary element analysis. Int J Fract 1987;34:131–47.
16