Numerical solution to stress distribution of a hole with corners on infinite plane

1 Introduction

Determining equilibrium shape with anisotropy surface energy is an important problem in material science. Wulff’s construction by Wulff (1901) gives a theorem describe equilibrium of an elastic solid with anisotropy surface. Many researchers also find analytically presentations of the equilibrium (Herring (1951), Cabrera (1964) and Pimpinelli and Villain (1999)). Their results show that when anisotropy is large, the equilibrium shape has corners.

Some researchers are focusing on the equilibrium shape if stress applied. Stress introduces the elastic strain energy in this case, which can change the equilibrium shape with anisotropy surface energy. Our goal is to solve the equilibrium shape of the void inside an elastically-stressed solid with anisotropic surface energy. There are several ways to solve the stress numerically on a shape in this corner case, such as adding a regularization term to smooth the corner and solve a smooth shape by Siegel et al. (2004), or applying a conformal transformation to a unit circle Soutas-Little (2010) approach the shape with corners e.g. rectangular hole Savin (1961), Pan et al. (2013) and solve the stress utilize the solution for infinite plane with a unit circle (Muskhelishvili (1953b)), where the shape is smooth as well.

However, these results are trying to avoid corner by smooth approach may not be so accurate on stress and elastic energy, especially near corners. Srolovitz and Davis (2001) discuss the elastic energy locally near the corner and shows stresses do not modify the corner angle. Siegel et al. (2004)gives a conflicting result by analysis the shape with regularization showing the corner angle changed by stress.

To find out if stress changes the corner angle, the solution of stress distribution on a given shape is essential. The stress determines the elastic energy on the boundary, and the shape with minimum energy is the exact shape. In this paper, we find a reliable solver with high accuracy to solve the stress distribution of the shape consider the effects of corners as the first step.

2 Mathematical formulation

Our formulation follows Muskhelishvili’s complex variable formulation on 2D-elasticity by Mikhlin (1957). We speak of a simply connected void inside elastic solid with biaxial stress applied far form the void. Assume the biaxial stress is on x and y direction and cause planar deformation. The displacement of elastic solid is: $u (x, y, z) = u_{1} (x, y) e_{1} + {u %}_{2} (x, y) e_{2}$ , where $e_{1}$ and $e_{2}$

are unit vector on x,y direction. Let

D

denote elastic solid region,

D^{'}

denote void region,

\partial D

represent the interface between void and solid. Displacement, stress tensor are defined on

D

Infinitesimal strain tensor is defined by $E = \frac{1}{2} (\nabla u + \nabla u^{T})$ , First Piola-Kirchhoff stress tensor of elastic solid is given by $P = λ (t r (E)) I + μ (E + E^{T})$ , where $t r (E)$ is trace of matrix $E$ and where $λ$ and $μ$ are Lamï¿œ coefficients . Thus,

P = (\begin{matrix} σ_{x} & τ_{x y} τ_{x y} & σ_{y} \end{matrix})

on xy-plane with far field condition:

(\begin{matrix} σ_{x} & τ_{x y} τ_{x y} & σ_{y} \end{matrix}) \to (\begin{matrix} 1 & 0 0 & χ \end{matrix}) a s \sqrt{x^{2} + y^{2}} \to \infty,

(1)

where $χ = σ_{1} / σ_{2}$ is parametre after nondimensionalization. Apply Law of Linear momentum in Lagrangian Form (here $x_{1}, x_{2}$ are $x, y$ ):

\nabla \cdot P = \sum i, j = 1, 2 \frac{\partial P_{i j}}{\partial x_{j}} e_{i} = 0.

(2)

Let $W (x, y)$ be a smooth function called stress function defined on $D$ , such that $σ_{x} = \partial^{2} W / \partial y^{2}$ , $τ_{x y} = - \partial^{2} W / \partial x \partial y$ , $σ_{y} = \partial^{2} W / \partial x^{2}$ . Then equation (1) hold automatically. This part of mechanics from Gonzalez and Stuart (2008). Combine the symmetry of strain tensor $E$ . $W (x, y)$ satisfies the following biharmornic equation:

\frac{\partial^{4} W}{\partial x^{4}} + 2 \frac{\partial^{4} W}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} W}{\partial y^{4}} = 0

(3)

Then using two function $ϕ (z)$ and $ψ (z)$ (called Goursat function) which is holomorphic on $D$ to represent $W (x, y)$ with complex variable $z = x + i y$ (different from z-direction here): $W (x, y) = R e {¯ z ϕ (z) + ς (z)}$ , and $ψ (z) = ς^{'} (z)$ . the relation between stress tensor and Goursat functions are:

σ_{x} + σ_{y} = 4 R e {ϕ^{'} (z)}

(4)

σ_{y} - σ_{x} + 2 i τ_{x y} = 2 [¯ z ϕ^{''} (z) + ψ^{'} (z)]

(5)

Boundary condition at infinity in terms of $φ$ and $ψ$ becomes:

ϕ (z) = \frac{1 + χ}{4} z + O (\frac{1}{z}) a s ∣ z ∣\to \infty

(6)

ϕ (z) = \frac{χ - 1}{2} z + O (\frac{1}{z}) a s ∣ z ∣\to \infty .

(7)

To make the solution unique also for convenience, the boundary condition at infinity in this paper drop some arbitrary constant on $φ$ and $ψ$ , which does not affect stresses (see Muskhelishvili (1953a) and Siegel et al. (2004)). Since no external force applied to $\partial D$ , the constraint on inner boundary $\partial D$ is given by

φ (z) + z ¯ ¯¯¯¯¯¯¯¯¯¯ ¯ φ^{'} (z) + ¯ ¯¯¯¯¯¯¯¯ ¯ ψ (z) = 0 o n z \in \partial D,

(8)

where Goursat functions $φ$ , $ψ$ can be writen as $ϕ (z) = (1 + χ) z / 4 + φ (z)$ and $ψ (z) = (χ - 1) z / 2 + h (z)$ . Take conjugate on both side, then farfield condition (6),(7) and boundary conditions (8) are equivalent to:

φ (z) \to 0 a s ∣ z ∣\to \infty

(9)

h (z) \to 0 a s ∣ z ∣\to \infty

(10)

¯ ¯¯¯¯¯¯¯¯ ¯ φ (z) + \frac{1 + χ}{4} ¯ ¯ ¯ z + ¯ ¯ ¯ z (φ^{'} (z) + \frac{1 + χ}{4}) + h (z) + \frac{χ - 1}{2} z = 0 o n z \in \partial D .

(11)

Purpose of making the substitution is to remove the singularity of $φ$ , $ψ$ at $\infty$ . $ϕ$ and $h$ are analytic on the region $D \cup \infty$ . Let us multiply bothside of (11) by the factor $1 / 2 π i \cdot d z / (z - t)$ , where $t$ is an arbitrary point on $D^{'}$ , and integrate along boundary $\partial D$ , $L$ is $\partial D$ with counterclockwise direction. Since $ϕ$ and $h$ are analytic on the region $D \cup \infty$ , $z$ is analytic on $D^{'}$ , by Cauchy integral formular, value of Cauchy integrals given by:

\frac{1}{2 π i} \int L \frac{φ (z)}{z - t} d z = 0 \frac{1}{2 π i} \int L \frac{h (z)}{z - t} d z = 0 \frac{1}{2 π i} L \int \frac{z}{z - t} d z = t .

Equation (11) becomes an integral equation does not involve $h (z)$ :

\frac{1}{2 π i} \int L \frac{¯ ¯¯¯¯¯¯¯¯ ¯ φ (z)}{z - t} d z + \frac{1 + χ}{2 π i \cdot 2} \int L \frac{¯ ¯ ¯ z}{z - t} d z + \frac{1}{2 π i} \int L \frac{¯ ¯ ¯ z φ^{'} (z)}{z - t} d z + \frac{χ - 1}{} 2 t = 0

(12)

Now let $t \to z_{0}$ , where $z_{0}$ is a point on boundary $\partial D$ , the limits has following property:

lim t \to z_{0} \frac{1}{2 π i} \int L \frac{¯ ¯¯¯¯¯¯¯¯ ¯ φ (z)}{z - t} d z = \frac{1}{2} ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ φ (z_{0}) + \frac{1}{2 π i} \int L \frac{¯ ¯¯¯¯¯¯¯¯ ¯ φ (z)}{z - z_{0}} d z

(13)

lim t \to z_{0} \frac{1}{2 π i} \int L \frac{¯ ¯ ¯ z}{z - t} d z = \frac{1}{2}_{0} + \frac{1}{2 π i} \int L \frac{¯ ¯ ¯ z}{z - z_{0}} d z

(14)

lim t \to z_{0} \frac{1}{2 π i} \int L \frac{¯ ¯ ¯ z φ^{'} (z)}{z - t} d z = \frac{1}{2}_{0} φ^{'} (z_{0}) + \frac{1}{2 π i} \int L \frac{¯ ¯ ¯ z φ^{'} (z)}{z - z_{0}} d z .

(15)

All integrals on the right are in sence of cauchy principle value. We arrived at an singular boundary integrodifferential equation on $\partial D$ :

\frac{1}{2} ¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ φ (z_{0}) + \frac{1}{2 π i} \int L \frac{¯ ¯¯¯¯¯¯¯¯ ¯ φ (z)}{z - z_{0}} d z + \frac{1 + χ}{4}_{0} + \frac{1 + χ}{4 π i} \int L \frac{¯ ¯ ¯ z}{z - z_{0}} d z

+ \frac{1}{2}_{0} φ^{'} (z_{0}) + \frac{1}{2 π i} \int L \frac{¯ ¯ ¯ z φ^{'} (z)}{z - z_{0}} d z + \frac{χ - 1}{2} z_{0} = 0.

(16)

Then the stress at any point $ζ$ inside solid can be derived by $φ (ζ) = 1 / 2 π i \cdot \int - L φ (z) d z / (z - ζ)$ . Since there is no body force in solid, equation (12) holds inside solid. Then we can apply equation (12) and find $ψ$ , $W$ , stress tensor $P$ at $ζ$ in solid.

3 Numerical method

The surface of the elastic solid is described by a closed continuous curve on xy-plane. For instance, we consider the prototype model for surface energy with fourfold anisotropy. Fourfold-symmetry anisotropy surface energy function form a fourfold symmetry surface. Surface is smooth and convex on each piece and has corners when anisotropy is significant (see Herring (1951)). Let the center of the void is origin, and corners are located on the intercepts between surface and axis. As we knew, Goursat function $ϕ (z)$ can change up to a complex constant even for the same stress solution. To make the solution unique, we assume $φ (z)$ is fourfold-symmetry and make some restrictions on it. Because of the fourfold symmetry, we only need to consider the $z_{0}$ on the first quadrant and make the extension to the entire surface when evaluating the integrals on (16).

3.1 Modelling the surface

Surface is represent in polar coordinate $r (θ)$ , where $r$ is radius from surface and $θ \in [0, \frac{π}{2}]$ denotes polar angle. If surface is given analytically, e.g $r (θ) = 1$ when surface is a circle or $r (θ) = b / \sqrt{1 - (e cos θ)^{2}}$ for a ellipse with eccentricity $e$ and semi-minor axes $b$ . In this case, $r^{'} (θ)$ and $r^{''} (θ)$ can be obtained directly without error. If surfac is not given explicitly. We use a series of Chebyshev polynomial to represent $r (θ) = [k = 1] N \sum c_{k} T_{k} (θ),$ where $T_{k} (θ)$ are $k$ th Chebyshev polynomial on $[0, \frac{π}{2}]$ . Coefficients $c_{k}$ are derived from value of $r$ at Chebyshev nodes on $θ \in [0, \frac{π}{2}]$

using Chebyshev interpolation in

Press et al. (1992). More details about error analysis of Chebyshev interpolation can be found in Boyd (2001). Since our surface shape is smooth on first quadrant, the error of interpolation are negligible when

N

is large.

Boundary value of complex function $φ (z)$ with complex variable $z (θ) = r (θ) cos θ + i \cdot r (θ) cos θ$ determines the stress distribution. Goursat function $φ$ is holomorphic on $D$ , which leads to $φ$ is a continuous function of $θ$ on boundary. Chebyshev basis also works in representing Goursat function on the first quadrant due to the smoothness of the shape on the first quadrant:

φ (θ) = [k = 1] N \sum a_{k} T_{k} (θ) + i \cdot [k = 1] N \sum b_{k} T_{k} (θ),

(17)

where $a_{k}$ , $b_{k}$ are unknowns. $φ (θ)$ can be extended to $θ \in [0, 2 π]$ by:

φ (θ) = - R e {φ (π - θ)} + i \cdot I m {φ (π - θ)} o n θ \in [π / 2, π]

φ (θ) = - R e {φ (θ - π)} - i \cdot I m {φ (θ - π)} o n θ \in [π, 3 π / 2]

φ (θ) = R e {φ (2 π - θ)} - i \cdot I m {φ (2 π - θ)} o n θ \in [3 π / 2, 2 π]

and $φ$ still a continuous function of $θ$ , implys continuity conditions on ends of interval:

R e {φ (θ)} = 0 a t θ = π / 2.

(18)

I m {φ (θ)} = 0 a t θ = 0 .

(19)

Some more constraints will be discussed in section 3.3 to ensure $φ (θ)$ is boundary value of an analytic function.

3.2 Discretization of integral equations

Equation (16) hold for any $z_{0}$ on $L$ . We pick $z_{i} = z (θ_{i})$ as collocation points, where $θ_{i} = π (x_{i} + 1) / 4$ and $x_{i}$ is the i-th root of degree ( $N - 1$ ) Legendre polynomial $P_{n - 1}$ (Legendre-Gauss quadrature points). The choice of number of collocation points is determined by number of unknowns. We implement both real part and imagine part of equation (16) on each $z_{i}$ ( $2 (N - 1)$ equations) together with two equations (18) and (19) to solve total $2 N$ unknowns $a_{k}$ , $b_{k}$ . Gauss–Legendre quadrature by Golub and Welsch (1969) gives good approximation to nonsingular integral by using value of function at Legendre-Gauss quadrature points:

\int_{0}^{π / 2} f (θ) d θ = [i = 1] N - 1 \sum ω_{i} f (θ_{i}),

(20)

where $ω_{i}$ are weight of Gauss–Legendre quadrature on $[0, π / 2]$ . Discretization of integral equations in our model is based on Gauss–Legendre quadrature. $φ (θ_{i})$ can be obtained by using Clenshaw’s recurrence formula from Clenshaw (1955) to minimize truncation error when make substitution $θ = θ_{i}$ .

Now we consider integrals in equation (12) term by term. Singularity of Cauchy integrals on equation (16) is separated by:

\int L \frac{¯ ¯¯¯¯¯¯¯¯ ¯ φ (z)}{z - z_{i}} d z = \int L \frac{¯ ¯¯¯¯¯¯¯¯¯¯ ¯ φ (z_{i})}{z - z_{i}} d z + \int L \frac{¯ ¯¯¯¯¯¯¯¯ ¯ φ (z) - ¯ ¯¯¯¯¯¯¯¯¯¯ ¯ φ (z_{i})}{z - z_{i}} d z

(21)

First term is

\int L \frac{¯ ¯¯¯¯¯¯¯¯¯¯ ¯ φ (z_{i})}{z - z_{i}} d z = ¯ ¯¯¯¯¯¯¯¯¯¯ ¯ φ (z_{i}) \cdot \int L \frac{1}{z - z_{i}} d z = π i ¯ ¯¯¯¯¯¯¯¯¯¯ ¯ φ (z_{i})

(22)

from boundary version of Cauchy integral formula. Second term takes the form in $θ$

\int L \frac{¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ φ (z) - φ (z_{i})}{z - z_{i}} d z = \int_{0}^{2 π} \frac{¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ φ (θ) - φ (θ_{i})}{z (θ) - z (θ_{i})} \cdot z^{'} (θ) d θ,

(23)

where $z^{'} (θ) = r^{'} (θ) cos θ - r (θ) sin θ + i (r^{'} (θ) sin θ + r (θ) cos θ)$ . Notice the singularity at $θ = θ_{i}$ is removable singularity canceled by taking limit at $θ_{i}$

lim θ \to θ_{i} \frac{¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯ φ (θ) - φ (θ_{i})}{z (θ) - z (θ_{i})} = \frac{¯ ¯¯¯¯¯¯¯¯¯¯¯ ¯ φ^{'} (θ_{i})}{z^{'} (θ_{i})} .

(24)

We find $φ^{'} (θ_{i})$ by using algorithm on Press et al. (1992) based on the relation between the Chebyshev coefficient and the Chebyshev coefficient of the derivative:

c_{i - 1}^{'} = c_{i + 1}^{'} + 2 (i - 1) c_{i} i \geq 1,

(25)

where $c_{i}$ are the Chebyshev coefficients and $c_{i}^{'}$ are Chebyshev coefficient of the derivatives. The Chebyshev coefficient of $φ^{'}$ is linear to unknowns $a_{k}$ , $b_{k}$ from equation (25). Then integral (23) can be evaluated using Gauss–Legendre quadrature and implement symmetry of $φ (θ)$ and $r (θ)$ . We can evaluate $\int L ¯ ¯ ¯ z d z / (z - z_{0})$ by letting $φ (θ) = z (θ)$ .

Then we consider integrodifferential term $\int L ¯ ¯ ¯ z φ^{'} (z) d z / (z - z_{0})$ in equation (16):

\int L \frac{¯ ¯ ¯ z φ^{'} (z)}{z - z_{i}} d z = \int L \frac{_{i} φ^{'} (z_{i})}{z - z_{i}} d z + L \int \frac{¯ ¯ ¯ z φ^{'} (z) -_{i} φ^{'} (z_{i})}{z - z_{i}} d z .

(26)

Follow the same logic, first term is

\int L \frac{_{i} φ^{'} (z_{i})}{z - z_{i}} d z = π i_{i} φ^{'} (z_{i}) .

(27)

Apply Gauss–Legendre quadrature on second integral in terms of $θ$ :

\int L \frac{¯ ¯ ¯ z φ^{'} (z) -_{i} φ^{'} (z_{i})}{z - z_{i}} d z = \int_{0}^{2 π} \frac{¯ ¯¯¯¯¯¯¯ ¯ z (θ) φ^{'} (θ) - ¯ ¯¯¯¯¯¯¯¯¯ ¯ z (θ_{i}) φ^{'} (θ_{i})}{} z (θ) - z (θ_{i}) \cdot z^{'} (θ) d θ,

(28)

the singularity at $θ = θ_{i}$ is removable singularity with

lim θ \to θ_{i} \frac{¯ ¯¯¯¯¯¯¯ ¯ z (θ) φ^{'} (θ) - ¯ ¯¯¯¯¯¯¯¯¯ ¯ z (θ_{i}) φ^{'} (θ_{i})}{z (θ) - z (θ_{i})} = \frac{¯ ¯¯¯¯¯¯¯¯¯¯ ¯ z^{'} (θ_{i}) φ^{'} (θ_{i}) + ¯ ¯¯¯¯¯¯¯¯¯ ¯ z (θ_{i}) φ^{''} (θ_{i})}{z^{'} (θ_{i})} .

(29)

Substitute our results removing singularity into (16), and we discretize equation in polar coordinate using Gauss-Legendre quadrature. Collocation points are $θ_{i}$ . Notice coefficients of derivative of a Chebyshev approximated function are linear to the coefficients of original Chebyshev approximated function, so $φ^{'} (θ_{i})$ and $φ^{''} (θ_{i})$ are linear in $a_{k}$ and $b_{k}$ . Gauss-Legendre quadrature gives linear combinations of function value at collocation points. Finally, we get a linear system in $a_{k}$ and $b_{k}$ .

3.3 Analyticity equations

In section 3.1, we assume $φ$ is a smooth function of $θ$ on first quadrant of boundary $\partial D$ from the definition (17). To make $φ$ an analytic function in $D,$ we define $φ (ζ)$ in $D$ by Cauchy integral: $φ (ζ) = 1 / 2 π i \int - L φ (z) d z / (z - ζ)$ .

φ (ζ) = \frac{1}{2 π i} \int - L \frac{φ (z)}{z - ζ} d z \to 0 a s ∣ ζ ∣\to \infty,

(30)

because $φ (z)$ is bounded on $\partial D$ which satisfy condition (9) at infinity. This definition of $φ (ζ)$ guarantee the analyticity in $D$ . To make $φ (ζ)$ analytic on $D \cup \partial D$ , $φ$ should be continuous to any point on boundary $\partial D$ which give us the continuous condition at each collocation points $z_{i}$ :

lim ζ \to z_{i} \frac{1}{2 π i} \int - L \frac{φ (z)}{z - ζ} d z = φ (z_{i})

take limit inside and write in Cauchy principle value,

\frac{1}{2 π i} \int L \frac{φ (z)}{z - z_{i}} d z + \frac{1}{2} φ (z_{i}) = 0.

(31)

The integral can be evaluated in the same way we find the Cauchy integral of $¯ ¯¯¯¯¯¯¯¯ ¯ φ (z)$ in section 3.2. Complex function $φ$ is analytic on $D \cup \partial D$ is equivalent to equation (31) hold on every point on $\partial D$ . Then equation (31) at each collocation points is a necessary condition to $φ (θ)$ is boundary value of an analytic function. Equation (31) at collocation points give analyticity constraints on coefficients. This integral is a linear combination of $a_{k}$ and $b_{k}$ , and so does $φ (z_{i}) = φ (θ_{i})$ .

Finally, we combine $2 (N - 1)$ equations from (31) at collocation points, $2 (N - 1)$ equations from discretized integral equation (16), and two equations (18), (19) from boundary conditions on both end (consider both of real and imaginary part of equations). Then solve $2 N$ unknowns $a_{k}$ and $b_{k}$ . The theorem of the uniqueness of 2D linear elastic problem is given on Knops and Payne (2012). Our model is a numerical discretization on elastic equation (16) with analyticity constraint, and drop constants cause nonuniqueness. The solution of unknows $a_{k}$ and $b_{k}$ is unique following the uniqueness of the original problem. We solve the non-square linear system by MATLAB matrix left division operator in Section 4.

4 Numerical results

Now we have numerical method solving Goursat function on boundary. In this section, we test our numerical method using boundary shape is circle, ellipse, overlapped circles.

4.1 Measurements of error

In case of circle or ellipse, exact solution of Goursat function on boundary are given by Muskhelishvili (1953a). We can compare $φ$ we get from our numerical method with exact one in sense of $L^{2}$ norm:

E r r o r_{L^{2}} = [\int_{0}^{π / 2} | φ (θ) - φ_{e x a c t} (θ) |^{2} d θ]^{1 / 2}

(32)

In case of overlapped circle, there is no exist result give exact solution of $φ$ . But in Spencer and Meiron (1994)the strain energy density takes the form

ω = \frac{1 - ν^{2}}{E} (σ_{x} + σ_{y})^{2}

where $E$ is Young’s modulus and $ν$ is Poission’s ratio of the material and the exact solution for $σ_{x} + σ_{y}$ is given by Ling (1948). Comparing $σ_{x} + σ_{y} = 4 R e {φ^{'} (z)}$ is a good mesurement for our method. Define $L^{2}$ norm error is:

E r r o r_{L^{2}} = [\int_{0}^{π / 2} ((σ_{x} + σ_{y}) - (σ_{x} + σ_{y})_{e x a c t})^{2} d θ]^{1 / 2}

(33)

4.2 Test for circle

The solution of exterior problem when interface is a unit circle in the case of uniaxal tension $χ = 0$ by using conformal map $z (ζ) = 1 / ζ$ to interior problem. $r (θ) = 1$ and $r^{'} (θ) = r^{''} (θ) = 0$ . Exact solution of $φ$ is given by $φ_{e x a c t} (z) = 1 / 2 / z$ . Absolute value of $a_{k}$ and $b_{k}$ are presented in Fig.1a in logy plot. Coefficient of $φ (z)$ reach magnitude near truncation error at $N = 15$ which ensure the accuracy of our collocation method. The convergence rate is represented by the $L^{2}$ error vesus $N$ plot in Fig.1b. Results illustrate the algorithm is spectrally convergence in number of collocation method $N$ when the boundary shape is a unit circle.

Fig. 1. (a) absolute value of coefficient of $a_{k}$ or $b_{k}$ versus index $k$ with the number of collocation points $N = 32$ . (b) $L^{2}$ error of $φ (z)$ versus number of collocation points $N$ . Here $L^{2}$ error is defined in Section 4.1.

4.3 Test for ellipse

If the ellipse is given on complex plane by conformal map $z (ζ) = 1 / ζ + m ζ$ to unit circle interior problem, where then exact solution is $φ (z) = (1 - m) / ζ$ . Eccentricity or ellipse is $e = \sqrt{1 - (1 - m)^{2} / (1 + m)^{2}}$ , and $r (θ) = (1 - m) / \sqrt{1 - (e cos θ)^{2}}$ . We can get $r^{'} (θ)$ and $r^{''} (θ)$ by taking derivative:

r^{'} (θ) = - \frac{1 - m}{2} (1 - e^{2} {cos}^{2} θ)^{- 3 / 2} e^{2} sin 2 θ

r^{''} (θ) = - \frac{1 - m}{2} [- 3 / 2 \cdot (1 - e^{2} {cos}^{2} θ)^{- 5 / 2} \cdot (e^{2} sin 2 θ)^{2} + 2 (1 - e^{2} {cos}^{2} θ)^{- 3 / 2} e^{2} cos 2 θ]

Then apply our solver to find approximation of $φ (z)$ , numerical results when $m = 0.5$ are shown below:

Fig. 2. (a) real part and imaginary part of exact solution $φ_{e x a c t}$ versus $θ \in [0, π / 2]$ . (b) absolute value of coefficient of $a_{k}$ or $b_{k}$ versus $k$ when number of collocation points $N = 96$ (c) $L^{2}$ error of $φ (z)$ versus number of collocation points $N$ .

4.4 Test for overlapped circles ( $α = π / 3$ )

When shape is given by two overlapped circles with same radius, parameter $α$ defined as on fig. 3(a) give the radian of circle overlapped to the other one. In this test, we let both circles are unit circles with $α = π / 3$ then the shape has two corners on y-axis.

Fig. 3(a) Overlapped circle shape when $α = π / 3$ . Fig. 3(b) $σ_{x} + σ_{y}$ vesus $θ \in [0, π / 2)$ when $α = π / 3$

The equation of this shape in polar coordinate is $r (θ) = cos α cos θ + \sqrt{1 - {sin}^{2} θ {cos}^{2} α}$ on first quadrant. Exact solution of $σ_{x} + σ_{y}$ is given by an integral equation:

4 K \int_{0}^{\infty} \frac{{sinh}^{2} s α - s^{2} {sin}^{2} α}{s (s^{2} + 1) (sinh 2 s α + s sin 2 α)} d s + 2 (N_{1} - N_{2}) \int_{0}^{\infty} \frac{s {sin}^{2} α}{sinh 2 s α + s sin 2 α} d u = N_{1}

(34)

σ_{x} + σ_{y} = 4 (cosh ξ - cos α) sin α \int_{0}^{\infty} \frac{2 K - (N_{1} - N_{2}) s (s - cot α coth s α)}{sinh 2 s α + s sin 2 α} sinh s α cos s ξ d s

(35)

where $ξ$ defined by:

cosh ξ = \frac{1 + cos α cos β}{cos α + cos β}

(36)

where $β = θ + arcsin (sin θ cos α)$ is center angle, $N_{1}$ is tension parallel to the x-axis and $N_{2}$ is tension parallel to the y-axis. The value of $σ_{x} + σ_{y}$ is evaluated by MATLAB numerical integration function ’integral’ (see Shampine (2008)). We test the case of longitudinal tension when $N_{1} = 1$ and $N_{2} = 0$ with $α = π / 3$ . Notice as $θ \to π / 2$ , $β \to π - α$ from (36), $ξ \to \infty$ in this case. Integral (34) contains highly oscilllating term $cos s ξ$ which cause significant error when taking numerical integration. For this reason, we drop a small interval $θ \in [π / 2 - 0.01 π, π / 2]$ when showing error and evaluating $L^{2}$ error.

Fig. 3(c) absolute value of coefficient of $a_{k}$ or $b_{k}$ versus number of collocation points $N$ . (d) Error of $σ_{x} + σ_{y}$ vesus $θ \in [0, π / 2 - 0.01 π]$ when $N = 256$ .

Fig. 3(b) and 3(c) demonstrate that our numrical method work even when the shape has corners. The error of our numerical method is at around $10^{- 8}$ when $N = 256$ . The coefficients of $φ$ is not decay so fast so the rate of convergence for overlapped circle with $α = π / 3$ is slower (see Fig. 3(e))

Fig. 3(d) $L^{2}$ error on $θ \in [0, π / 2 - 0.01 π]$ versus number of collocation points $N$ .

4.5 Test for overlapped circles ( $α = 2 π / 3$ )

The shape of overlapped circles when $α = 2 π / 3$ is on Fig. 4(a). When the shape has corners with acute angle, the stress has singularity at corner (Williams (1952)). Our result approach the exact solution but with large error on the end of intervals (Runge’s phenomenon) because $φ \to \infty$ at corner.

Fig. 3(a) Overlapped circle shape when $α = 2 π / 3$ . Fig. 3(b) Numerical solution and exact solution of $σ_{x} + σ_{y}$ vesus $θ \in [0, π / 2)$ when $α = π / 3$

5 Summary and discussion

We developed a numerical method solving elasticity equations based on the Muskhelishvili’s theorem in this paper. From the results in section 4, our numerical method can handle the smooth shape and approach the exact solution in high accuracy with a small number of collocation points. Even when shape has corners, our solver also converges to the exact solution spectrally. Error on section 4.4 comes from two possible resources:

1. Error from our numerical method. e.g., the contribution of $a_{k}$ and $b_{k}$ for large $N$ , the error comes from taking derivative of $φ$ , truncation error for large $N$ . Constraints on analyticity may not adequate.

2. Error from ’Exact solution.’ We know the exact solution of overlapped circles from (35). The exact solution has some integral with infinite bounds and oscillating function. We can not estimate the error when using numerical integration to find the value because we can not find the correct answer explicitly. However, there must be errors on numerical integral. Section 4.5 is an example of the corner with an acute angle, which is a common shape with anisotropic surface energy. In this case, our model could not find the correct solution of stresses because of the singularity of

φ

at corner. We will investigate the singularity at corners and separate them out. Then we can solve this case using our algorithm.

Our collocation method can give stress distribution in high accuracy from given fourfold shape with corners. It is possible to implement our method in a non-fourfold shape by change the boundary condition (18), (19) to some other condition like periodic condition determined by the shape. Since our final goal is to find a convex equilibrium shape with anisotropy surface energy, we use $θ$ as the parameter to simplify some calculation when evaluating the area of the shape. Using arclength $s$ as parameter is an option if the shape can not be represented on the polar coordinate.

We still need some further work on how constraints on analyticity work: The difference between the smooth function and boundary value of an analytic function. Constraints on analyticity give us $2 (N - 1)$

equations but contain too many degrees of freedom. Some investigation on space span by smooth function, space of boundary value of an analytic function, space of complex function span by Chebyshev basis.

References

Boyd (2001) Boyd, J.P., 2001. Chebyshev and Fourier spectral methods. Courier Corporation.
Cabrera (1964) Cabrera, N., 1964. THE EQUILIBRIUM OF CRYSTAL SURFACES. Technical Report. VIRGINIA UNIV CHARLOTTESVILLE.
Clenshaw (1955) Clenshaw, C.W., 1955. A note on the summation of chebyshev series. Mathematics of Computation 9, 118–120.
Golub and Welsch (1969) Golub, G.H., Welsch, J.H., 1969. Calculation of gauss quadrature rules. Mathematics of computation 23, 221–230.
Gonzalez and Stuart (2008) Gonzalez, O., Stuart, A.M., 2008. A first course in continuum mechanics. Cambridge University Press. pp. 311–312.
Herring (1951) Herring, C., 1951. Some theorems on the free energies of crystal surfaces. Physical review 82, 87.
Knops and Payne (2012) Knops, R.J., Payne, L.E., 2012. Uniqueness theorems in linear elasticity. volume 19. Springer Science & Business Media.
Ling (1948) Ling, C.B., 1948. The stresses in a plate containing an overlapped circular hole. Journal of Applied Physics 19, 405–411.
Mikhlin (1957) Mikhlin, S., 1957. Integral equations, translated from russian by ah armstrong.
Muskhelishvili (1953a) Muskhelishvili, N.I., 1953a. Some basic problems of the mathematical theory of elasticity. P. Noordhoff, Groningen.
Muskhelishvili (1953b) Muskhelishvili, N.I., 1953b. Some basic problems of the mathematical theory of elasticity. P. Noordhoff, Groningen. pp. 320–331.
Pan et al. (2013) Pan, Z., Cheng, Y., Liu, J., 2013. Stress analysis of a finite plate with a rectangular hole subjected to uniaxial tension using modified stress functions. International Journal of Mechanical Sciences 75, 265–277.
Pimpinelli and Villain (1999) Pimpinelli, A., Villain, J., 1999. Physics of crystal growth. Physics of Crystal Growth, by Alberto Pimpinelli and Jacques Villain, pp. 400. ISBN 0521551986. Cambridge, UK: Cambridge University Press, February 1999. , 43–59.
Press et al. (1992) Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical recipes in Fortran 77: the art of scientific computing. Cambridge university press Cambridge. volume 2. pp. 150–153.
Savin (1961) Savin, G.N., 1961. Stress distribution around holes.
Shampine (2008) Shampine, L.F., 2008. Vectorized adaptive quadrature in matlab. Journal of Computational and Applied Mathematics 211, 131–140.
Siegel et al. (2004) Siegel, M., Miksis, M., Voorhees, P., 2004. Evolution of material voids for highly anisotropic surface energy. Journal of the Mechanics and Physics of Solids 52, 1319–1353.
Soutas-Little (2010) Soutas-Little, R.W., 2010. Elasticity. Dover Publication. pp. 195–226.
Spencer and Meiron (1994) Spencer, B., Meiron, D., 1994. Nonlinear evolution of the stress-driven morphological instability in a two-dimensional semi-infinite solid. Acta Metallurgica et Materialia 42, 3629–3641.
Srolovitz and Davis (2001) Srolovitz, D., Davis, S.H., 2001. Do stresses modify wetting angles? Acta materialia 49, 1005–1007.
Williams (1952) Williams, M., 1952. Stress singularities resulting from various boundary conditions in angular corners of plates in extension. Journal of applied mechanics 19, 526–528.
Wulff (1901) Wulff, G., 1901. Xxv. zur frage der geschwindigkeit des wachsthums und der auflösung der krystallflächen. Zeitschrift für Kristallographie-Crystalline Materials 34, 449–530.

Numerical solution to stress distribution of a hole with corners on infinite plane

Application of orthonormal Bernoulli polynomials for approximate solution of some Volterra integral equations

Spectral formulation of the boundary integral equation method for antiplane problems

Numerical solution to the 3D Static Maxwell equations in axisymmetric singular domains with arbitrary data

Tiling Rectangles and the Plane using Squares of Integral Sides

Estimation of non-symmetric and unbounded region of attraction using shifted shape function and R-composition

A continuous fracture front tracking algorithm with multi layer tip elements (MuLTipEl) for a plane strain hydraulic fracture

A projected Nesterov-Kaczmarz approach to stellar population-kinematic distribution reconstruction in Extragalactic Archaeology

1 Introduction

2 Mathematical formulation