# A Non-Convex Partition of Unity and Stress Analysis of a Cracked Elastic Medium

Authors: Won-Tak Hong
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

A stress analysis using a mesh-free method on a cracked elastic medium needs a partition of unity for a non-convex domain whether it is defined explicitly or implicitly. Constructing such partition of unity is a nontrivial task when we choose to create a partition of unity explicitly. We further extend the idea of thealmost everywhere partition of unity and apply it to linear elasticity problem. We use a special mapping to build a partition of unity on a non-convex domain. The partition of unity that we use has a unique feature: the mapped partition of unity has a curved shape in the physical coordinate system. This novel feature is especially useful when the enrichment function has polar form, f(r,θ)=rλg(θ), because we can partition the physical domain in radial and angular directions to perform a highly accurate numerical integration to deal with edge-cracked singularity. The numerical test shows that we obtain a highly accurate result without refining the background mesh.


---

## Body

## 1. Introduction

Many mesh-free methods [1–9] showed great success to solve challenging problems. Whether it has explicitly defined a partition of unity function or not, most of the mesh-free methods can be categorized as a partition of unity method. Some of them [2–5, 8, 10, 11] use partition of unity explicitly by utilizing background mesh. The partition of unity function allows researchers to use a priori knowledge of the solution and enables enriching the solution space with special functions [8, 12–14]. The main interest of such enrichment was focused on crack modeling [1, 2, 7, 8, 15, 16]. The partition of unity enrichment with discontinuous functions allows to model discontinuity in the domain such as a crack in an elastic medium. Thus, using the Heaviside enrichment function, researchers could avoid non-convex partitioning at all, but an accurate numerical integration becomes challenging because of the discontinuous enrichment [16–21]. The discontinuous enrichment, however, was an ingenious tool to model dynamic crack propagation, since it can be replaced easily without changing the given background mesh [2]. Thus, engineers and researchers could save time by enrichment and avoid remeshing the entire domain. On the other hand, there has been different effort to deal with a non-convex domain. Instead of enriching with discontinuous function, a mapping technique has been used to partition the given convex domain and coined almost everywhere partition of unity [22] because of failing partition of unity at some points at the boundary. The later effort focused more on accuracy and error estimation of elliptic boundary value problem, but elasticity equation has not yet been explored.The goal of this study is performing accurate stress analysis on the cracked elastic medium with the almost everywhere partition of unity. It has been well known that the singularity can cause poor performance of finite element approximation for the elliptic boundary value problems [23–26], including elasticity problems. The quality of enriched mesh-free method depends on the accuracy of the numerical integration [15]. The challenge for elasticity equation is the nonsmooth integrand due to the enrichment. Thus, to tackle the difficulty related to the numerical integration after enrichment, we use a mapping technique rather than using a discontinuous enrichment because the later requires integrating discontinuous function, which complicates the numerical integration. Another reason that we favor partitioning non-convex domain explicitly over discontinuous enrichment is the matrix condition number. Most of the study with discontinuous enrichment uses finite element mesh as the partition of unity function, and it could result in unexpected linear dependency with polynomial enrichment [3, 27, 28]. We use a flat-top partition of unity function to avoid linear dependency and thus keep the growth of matrix condition number reasonable with the polynomial and nonpolynomial enrichment.With the almost everywhere partition of unity as defined in [22], it still has difficulty in dealing with nonsmooth enrichment function that is given in polar form because the supports of the partition of unity functions are quadrangular. Thus, upon defining a curved annular shaped partition of unity on a non-convex domain, we can integrate the enriched function not only in radial direction but also in angular direction. As a result, with the help of mapped numerical quadrature, we obtain highly accurate computed stresses.


## 2. Preliminaries

Defining a partition of unity on a given domain is the most important part of a mesh-free method. Among many partition of unity functions, we adopt the following standard definitions and theorems introduced in [28].Definition 1.
For a point finite open covering{Ok:k∈Λ} of a domain Ω⊂Rd, a family of Lipschitz functions, {ϕk:k∈Λ}, is called a partition of unity subordinate to the covering {Ok:k∈Λ} if supp(ϕk)⊆O¯k for each k∈Λ and there exists a number C such that ϕk∞,Rd≤C for each k with the following condition:(1)∑k∈Λϕkx=1for  eachx∈Ω.Definition 2.
For integersn≥1, we define a piecewise polynomial function by(2)ϕgnx=1+xngnxx∈-1,0,1-xngn-xx∈0,1,0x≥1,where gn(x)=a0(n)+a1(n)(-x)+a2(n)(-x)2+⋯+an-1(n)(-x)n-1 whose coefficients are inductively constructed by the following recursion formula:(3)akn=1,ifk=0,∑j=0kajn-1,if0<k≤n-2,2an-2n,ifk=n-1.Let us defineϕgnL(x)=(1+x)ngn(x) and ϕgnR(x)=(1-x)ngn(-x). The superscript R and L stand for right and left simple polynomials, respectively. Then, with these simple polynomials, we construct a piecewise Cn-1 polynomial function with flat-top property in one dimension as follows:(4)ψa,bδ,n-1=ϕgnLx-a+δ2δifx∈a-δ,a+δ,1,ifx∈a+δ,b-δ,ϕgnRx-b-δ2δifx∈b-δ,b+δ,0,ifx∉a-δ,b+δ.The above partition of unity with flat-top is exactly the same as the convolution of the characteristic function χ[a,b] and the scaled window function defined by w(x)=A(1-(x/δ)2)n for |x|<δ, 0 for |x|>δ. A is a constant that makes ∫-∞∞‍wdx=1 [28]. The most obvious way to construct a higher dimension flat-top partition of unity function, although its support is rather restricted to be rectangular, is using a tensor product of one-dimensional partition of unity functions [28]:(5)Ψx,y=ψa,bδ,n-1x×ψc,dδ,n-1y.Other studies [11, 29] generalized the tensor product idea to a quadrangular shape. In this study, we use a partition of unity for annular shape background mesh to deal with non-convex domain.


## 3. The Linear Elasticity Equations

LetΩ be an elastic medium in R2 with boundary Γ=∂Ω. The displacement vectors are denoted by u={ux(x,y),uy(x,y)}T and stress tensor are denoted by {σ(u)}={σx(u),σy(u),τxy(u)}T. Let us denote the strain tensor {ε(u)}={εx(u),εy(u),γxy(u)}. Then the strain-displacement and stress-strain relations are given by(6)εu=Du,(7)σu=Eεu,respectively, where [D] is the following differential matrix:(8)
D
=
∂
∂
x
0
0
∂
∂
y
∂
∂
y
∂
∂
x
,and [E]=[Eij],1≤i,j≤3, is a symmetric positive definite matrix of material constants. For an isotropic elastic medium, the matrix [E] is given as(9)E=E1-ν21ν0ν10001-ν2when plane stress is the quantity of interest or(10)E=ζ+2μζ0ζζ+2μ000μwhen plane strain is the quantity of interest. Here, μ and ζ are defined as follows:(11)μ=E21+ν,ζ=νE1+ν1-2ν.E is the modulus of elasticity and ν(0≤ν<1/2) is Poisson’s ratio.The equilibrium equations of elasticity are given by(12)DTσux,y+fx,y=0,x,y∈Ω,where {f}={fx(x,y),fy(x,y)}T is the vector that describes body force per unit area.Using the stress-strain relation, (6) and (7), we can express (12) in terms of the displacement vector {u}. Let us consider the following system of partial differential equations in terms of the displacement vector:(13)DTEDux,y+fx,y=0,x,y∈Ω,subject to the following boundary conditions,(14)Nσus=T~s=T~xs,T~ysT,s∈ΓN,us=u~s=u~xs,u~ysT,s∈ΓD,where Γ=ΓN∪ΓD=∂Ω, {nx,ny}T is an outward unit normal vector to the traction boundary ΓN, and(15)N=nx0ny0nynx.LetHD1(Ω)={{w}={wx,wy}∈[H1(Ω)]2:{w}=0onΓD}. Then the variational equation of the linear elasticity equation, (13) and (14), becomes as follows.Find the vector{u}={ux,uy} such that ux,uy∈H1(Ω), {u}={u~} on ΓD, and(16)Bu,v=Fv∀v∈HD1Ω,where(17)Bu,v=∫ΩDvTEDudxdy,(18)Fv=∫ΩvTfdxdy+∮ΓNvTT~ds.U({u})=(1/2)B({u},{v}) is the strain energy of the displacement vector {u}.Let us denote the basis functions defined onΩ by Φix,y,i=1,2,…,n. The components of the displacement vector in terms of basis functions Φi are given as the following forms:(19)uxx,y=∑k=1nckΦkx,y,uyx,y=∑k=1ncn+kΦkx,y,where ci(i=1,2,…,2n) are called the amplitudes of the basis functions Φi. Let(20)Φk=Φkx,y0,k=1,2,…,n,Φk=0Φkx,y,k=n+1,2,…,2n.Then the displacement vector {u} can be written as (21)u=∑k=12nakΦk.Substituting (20) into (17), we have the following.Lemma 3.
The bilinear formB({u},{v}) on e, which is the common support of Φi and Φj, becomes one of the following four equations:(22)∫e∇ΦjTE11E13E31E33∇ΦidΩ,ifv=Φi,0T,u=Φj,0T,∫e∇ΦjTE33E32E23E22∇ΦidΩ,ifv=0,ΦiT,u=0,ΦjT,∫e∇ΦjTE13E12E33E32∇ΦidΩ,ifv=Φi,0T,u=0,ΦjT,∫e∇ΦjTE31E33E21E23∇ΦidΩ,ifv=0,ΦiT,u=Φj,0T.


## 4. A Mesh-Free Approximation for the Linear Elasticity Equation on an Edge-Cracked Domain

We define themesh-free approximation spaceVh,k as follows:(23)Vh,k=spanΦix,y=ΨJix,yfIikx,y,i=1,…,Nk,J=1,…,M,where ΨJ(i)(x,y)=ΨJ is the partition of unity function and fI(i)k(x,y) is the local approximation function that has polynomial reproducing order k. Standard h-version of finite element method has polynomial reproducing order 1. We refer to [28] for different local approximation functions that can be used. In this study, we use Lagrange interpolating polynomials of degree k and some singular functions for the local approximation function.Using the mesh-free approximation spaceVh,k, we state themesh-free approximation of the linear elasticity equation as follows.Find the vector{uh}={uhx,uhy} such that uhx,uhy∈Vh,k, {uh}={u~h} on ΓD, and(24)Buh,vh=Fvh∀vh∈VDh,k,where the bilinear forms B and F are defined as in (17) and (18).We define theenriched mesh-free approximation space  Venh,k as follows:(25)Venh,k=Vh,k⊕spanΦix,y=ΨJix,ysix,y,i=1,…,NJ,J=1,…,M,where si(x,y) is a special function that is enriched on the patches with supp(ΨJ(i)), which contain the crack tip.Then theenriched mesh-free approximation of the linear elasticity equation becomes as follows.Find the vector{uh}={uhx,uhy} such that uhx,uhy∈Venh,k, {uh}={u~h} on ΓD, and(26)Buh,vh=Fvh∀vh∈Ven,Dh,k,where the bilinear forms B and F are defined as in (24).Using Lemma3, we see that the bilinear form B({u},{v}) on the common support e can be obtained using the following integral:(27)∫e∇x,yΦix,yT∇x,yΦjx,ydΩ.If we use Venh,k, (27) becomes either one of the four following equations:(28)∫e∇x,yΨJifIikT∇x,yΨJjfIjkdxdy,(29)∫e∇x,yΨJifIikT∇x,yΨJjsIjdxdy,(30)∫e∇x,yΨJisIiT∇x,yΨJjfIjkdxdy,(31)∫e∇x,yΨJisIiT∇x,yΨJjsIjdxdy,where i and j are global basis number and I and J are index functions depending on the global basis number. If fIk is polynomial, then the integrands in (28) are piecewise polynomial so we can accurately evaluate the integral using numerical quadratures.IfSI(j) is an enrichment of type rg(θ), then the integrals in (29)–(31) have singular integrands of type r-1/2 or r-1. However, if we use the change of variable given in (33), these terms are changed to r2 and r, respectively. Thus, the stiffness calculation given in (27) is regular on the (r,θ)-coordinate and we could perform highly accurate numerical integration.


## 5. Numerical Examples

We consider two examples. We useC2 partition of unity function defined in (4) throughout this section. Also, we fix the polynomial reproducing order, k in (23), when comparing the enriched mesh-free solution with the finite element solution without enrichment.Example 1.
Let us consider the equation of elasticity on a domainΩ={(x,y):-2≤x,y≤2} with a crack along the negative x-axis. Assume that Young’s modulus E=1000 and Poisson’s ratio ν=0.3. We impose the following stresses along all boundaries of the given domain:(32)σx=14r3cos⁡θ2+cos⁡5θ2,σy=14r5cos⁡θ2-cos⁡5θ2,τxy=14rsin⁡5θ2-sin⁡θ2.The displacement vector {u} is fixed at the crack tip (0,0) and y-component of the displacement vector, uy(x,y), is fixed at (2,2) to prevent rigid body motion.Let us illustrate how to construct a partition of unity on a non-convex domain. We consider the background mesh pictured in Figure1(a). We use eight patches QJ,J=I,II,…,VIII as shown in Figure 1(a). For each patch QJ, we construct a partition of unity function ΨJ that corresponds to the patch QJ. The support of the partition of unity functions ΨJ is shown in Figure 1(b). Note that the support overlaps. The thin strips that are enclosed by dotted lines are the overlapping region of the partition of unity functions. For example, the support of the third partition of unity functions, ΨIII, is the shaded region in Figure 2. Note that the partition of unity functions does not overlap along the negative x-axis where the crack is located.Figure 1
An edge-cracked domain: (a) the background mesh and eight patches and (b) the overlapped supports of partition of unity functions.
(a) (b)Figure 2
Polar mapping, given in (33), maps the reference patch to the physical patch. The third reference patch is shaded.LetP:(r,θ)→(x,y) be the coordinate transformation:(33)x=r2cos⁡θ;y=r2sin⁡θ.Then, as shown in Figure 2, through the mapping P, the partition of unity on the cracked domain in the (r,θ)-coordinate system is generated by using the partition of unity of the reference domain in the rectangular coordinate system. The four outer patches QJ,J∈{V,…,VIII}, are mapped to quadrangles Q^J,J∈{V,…,VIII}, with one curved side.Let the partition of unity function on the reference coordinate system beΨ^J(r,θ), where J∈{I,II,…,VIII}. This reference partition of unity functions Ψ^J(r,θ) can be easily constructed by using the tensor product of one-dimensional flat-top partition of unity function (4); (34)Ψ^Jr,θ=ψδ,n-1r×ψδ,n-1θ.Then the partition of unity functions ΨJ(x,y) on the physical coordinate system is given by(35)ΨJx,y=Ψ^J∘P-1x,y,where J∈{I,II,…,VIII}.Next, we enrich the following four singular functions in the patch that contains the crack tip.(36)rcos⁡θ2,rsin⁡θ2,rcos⁡5θ2,rsin⁡5θ2.A prior knowledge to the solution behavior or a numerical experiment could provide an optimal or nearly optimal choice for the enrichment functions.In order to show the effectiveness of the enrichment, we compare the enriched mesh-free solution with plain mesh-free solution with polynomial reproducing order four. The plain mesh-free solution with polynomial reproducing order four is comparable to the standardp-version finite element methods with p=4. Figure 3 shows the mesh-free solution with and without enrichment along with true stress τxy near the crack tip, [0.001,0.01]×[-π,π]. We clearly see the effectiveness of the enrichment in handling crack singularities. There is virtually no difference between the true stress and the computed stress obtained by enriched mesh-free method.Figure 3
The stressτxy comparison: (a) the enriched mesh-free method with crack-tip enrichment; (b) the mesh-free method without crack-tip enrichment; and (c) the true stress.
(a) (b) (c)Example 2.
Let us consider the equations of elasticity on a domain shown in Figure4, which is isotropic with material constants: E=1000 and ν=0.3. The boundary conditions are given as follows: un=0, ut=0 along Γ1∪Γ2, Tn=10, Tt=2 along Γ5, Tn=0, Tt=0 along ∂Ω∩(Γ1∪Γ2∪Γ5)c.Figure 4
Isotropic cracked elastic domain with a weak singularity at(2,2) and a strong singularity at the crack tip, (0,0).We use Mesh I and Mesh II, shown in Figures5(a) and 5(b), for the p-version of finite element method. MAM (Method of Auxiliary Mapping) [30] is powerful tool that can be used with p-version finite element method when the problem has known singularity. MAM uses the same Mesh I but special mappings are used on the shaded regions; see Figure 5(c). The given domain, Figure 4, is essentially the same as Example 1, except the crack lying in the positive x-axis instead of negative x-axis. We use the background mesh given in Figure 5(d) to construct partition of unity functions. The eight enrichment patches are shaded. For this example, we use only one singular enrichment function, r1/2sin⁡(θ/2). The degrees of freedom are compared in Table 1 for four different methods. We see that the degree of freedom for the enriched mesh-free method is the smallest, except p=2, and has the simplest mesh setup.Table 1
Degrees of freedom comparison.
p p-FEM with Mesh I & MAM p-FEM with Mesh II Meshfree enrichment 2 120 280 128 4 376 792 328 6 808 1688 624 8 1416 2968 —Figure 5
Meshes forp-FEM and background mesh for the mesh-free method: (a) Mesh I (22 elements); (b) Mesh II (48 elements); (c) MAM with Mesh I and auxiliary mapping zone is shaded; and (d) background mesh for enriched mesh-free method.
(a) (b) (c) (d)We calculate the strain energy of computed solutionuh by U(uh)=(1/2)B(uh,uh), where B is the bilinear form defined in (17). The computed energy for four different methods is summarized in Table 2. The true energy Uexact for this problem is 2.113815563245032 and can be obtained by the extrapolation technique [24]. We calculate the relative error in energy norm as follows:(37)eE,rel=Uexact-UuhUexact1/2.Table 3 lists the computed relative errors in percent. The enriched mesh-free method is not only robust but also highly efficient at handling singularity arising in elasticity problem on a cracked domain. With the right enrichment function, we see that the approximate solution of the enriched mesh-free method converges much faster than MAM with fewer degrees of freedom.Table 2
Total strain energy comparison for Example2.
p p-FEM with Mesh I p-FEM with Mesh II p-FEM with MAM Meshfree enrichment 2 1.60551078343838 2.01871772810441 1.70213202981048 2.019413351175610 4 1.81518650351132 2.08274831820450 2.11004639654621 2.112748709582020 6 1.89915218378087 2.09344129702200 2.11356074177201 2.113806861694690 8 1.94532926077425 2.09816367972923 2.11378584006680 —Table 3
Relative error comparison in the energy norm (%).
p p-FEM with Mesh I p-FEM with Mesh II p-FEM with MAM Meshfree enrichment 2 49.04 21.21 15.81 21.3 4 37.59 12.12 4.22 2.25 6 31.87 9.82 1.10 0.20 8 28.23 8.60 0.38 —


## 6. Concluding Remarks

We have used a mapping technique, given in (33), to deal with the crack singularity. In this study, we map the entire computational domain to the (r,θ)-coordinate system to remove the r-singularity. A domain that has more general shape can be mapped using conformal mapping technique [17]; however, essential boundary condition imposition becomes nontrivial with the current almost everywhere partition of unity. Thus, the idea of the mapping method in this study can be more powerful if we can avoid global coordinate mapping so that the mapping can be used only locally where it is needed. To that end, a coupling method [31] or generalized product partition of unity [29] can be utilized together to localize the mapping technique. Also the same mapping technique can be used to improve singular integral given in [14].



---

*Source: 9574341-2017-02-07.xml*