Log In Sign Up

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

In this paper, we consider the stress of a hole with the given fourfold shape (with corner) on an infinite plane under uniaxial tension. Complex Goursat functions formulation by Muskhelishvili (1953) gives a set of singular integral equations on the boundary to solve this problem. We develope a numerical method using a set of Chebyshev polynomial with some constraints to represent the Goursat function on the boundary and apply the collocation method on roots of Legendre polynomial to solve integral equations. Our results show that the numerical method spectrally converges to the known exact solution when boundary shape is a circle, ellipse. We also applied our numerical method on two overlapped circle shape (with corner) and find the result also converge to the exact solution on Ling (1948).


page 1

page 2

page 3

page 4


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

In this work, a new approach has been developed to obtain numerical solu...

Spectral formulation of the boundary integral equation method for antiplane problems

A spectral formulation of the boundary integral equation method for anti...

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

We propose a numerical method to solve the three-dimensional static Maxw...

Tiling Rectangles and the Plane using Squares of Integral Sides

We study the problem of perfect tiling in the plane and exploring the po...

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

A general numerical method using sum of squares programming is proposed ...

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: , where and

are unit vector on x,y direction. Let

denote elastic solid region, denote void region,

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


Infinitesimal strain tensor is defined by , First Piola-Kirchhoff stress tensor of elastic solid is given by , where is trace of matrix and where and are Lamᅵ coefficients . Thus,

on xy-plane with far field condition:


where is parametre after nondimensionalization. Apply Law of Linear momentum in Lagrangian Form (here are ):


Let be a smooth function called stress function defined on , such that , , . Then equation (1) hold automatically. This part of mechanics from Gonzalez and Stuart (2008). Combine the symmetry of strain tensor . satisfies the following biharmornic equation:


Then using two function and (called Goursat function) which is holomorphic on to represent with complex variable (different from z-direction here): , and . the relation between stress tensor and Goursat functions are:


Boundary condition at infinity in terms of and becomes:


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 , the constraint on inner boundary is given by


where Goursat functions , can be writen as and . Take conjugate on both side, then farfield condition (6),(7) and boundary conditions (8) are equivalent to:


Purpose of making the substitution is to remove the singularity of , at . and are analytic on the region . Let us multiply bothside of (11) by the factor , where is an arbitrary point on , and integrate along boundary , is with counterclockwise direction. Since and are analytic on the region , is analytic on , by Cauchy integral formular, value of Cauchy integrals given by:

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


Now let , where is a point on boundary , the limits has following property:


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


Then the stress at any point inside solid can be derived by . Since there is no body force in solid, equation (12) holds inside solid. Then we can apply equation (12) and find , , stress tensor 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 can change up to a complex constant even for the same stress solution. To make the solution unique, we assume is fourfold-symmetry and make some restrictions on it. Because of the fourfold symmetry, we only need to consider the 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 , where is radius from surface and denotes polar angle. If surface is given analytically, e.g when surface is a circle or for a ellipse with eccentricity and semi-minor axes . In this case, and can be obtained directly without error. If surfac is not given explicitly. We use a series of Chebyshev polynomial to represent where are th Chebyshev polynomial on . Coefficients are derived from value of at Chebyshev nodes on

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 is large.

Boundary value of complex function with complex variable determines the stress distribution. Goursat function is holomorphic on , 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:


where , are unknowns. can be extended to by:

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


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 on . We pick as collocation points, where and is the i-th root of degree () Legendre polynomial (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 ( equations) together with two equations (18) and (19) to solve total unknowns , . Gauss–Legendre quadrature by Golub and Welsch (1969) gives good approximation to nonsingular integral by using value of function at Legendre-Gauss quadrature points:


where are weight of Gauss–Legendre quadrature on . Discretization of integral equations in our model is based on Gauss–Legendre quadrature. can be obtained by using Clenshaw’s recurrence formula from Clenshaw (1955) to minimize truncation error when make substitution .

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


First term is


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


where . Notice the singularity at is removable singularity canceled by taking limit at


We find by using algorithm on Press et al. (1992) based on the relation between the Chebyshev coefficient and the Chebyshev coefficient of the derivative:


where are the Chebyshev coefficients and are Chebyshev coefficient of the derivatives. The Chebyshev coefficient of is linear to unknowns , from equation (25). Then integral (23) can be evaluated using Gauss–Legendre quadrature and implement symmetry of and . We can evaluate by letting .

Then we consider integrodifferential term in equation (16):


Follow the same logic, first term is


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


the singularity at is removable singularity with


Substitute our results removing singularity into (16), and we discretize equation in polar coordinate using Gauss-Legendre quadrature. Collocation points are. Notice coefficients of derivative of a Chebyshev approximated function are linear to the coefficients of original Chebyshev approximated function, so and are linear in and . Gauss-Legendre quadrature gives linear combinations of function value at collocation points. Finally, we get a linear system in and .

3.3 Analyticity equations

In section 3.1, we assume is a smooth function of on first quadrant of boundary from the definition (17). To make an analytic function in we define in by Cauchy integral: .


because is bounded on which satisfy condition (9) at infinity. This definition of guarantee the analyticity in . To make analytic on , should be continuous to any point on boundary which give us the continuous condition at each collocation points :

take limit inside and write in Cauchy principle value,


The integral can be evaluated in the same way we find the Cauchy integral of in section 3.2. Complex function is analytic on is equivalent to equation (31) hold on every point on . 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 and , and so does .

Finally, we combine equations from (31) at collocation points, 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 unknowns and . 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 and 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 norm:


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

where is Young’s modulus and is Poission’s ratio of the material and the exact solution for is given by Ling (1948). Comparing is a good mesurement for our method. Define norm error is:


4.2 Test for circle

The solution of exterior problem when interface is a unit circle in the case of uniaxal tension by using conformal map to interior problem. and . Exact solution of is given by . Absolute value of and are presented in Fig.1a in logy plot. Coefficient of reach magnitude near truncation error at which ensure the accuracy of our collocation method. The convergence rate is represented by the error vesus plot in Fig.1b. Results illustrate the algorithm is spectrally convergence in number of collocation method when the boundary shape is a unit circle.

Fig. 1. (a) absolute value of coefficient of or versus index with the number of collocation points . (b) error of versus number of collocation points . Here error is defined in Section 4.1.

4.3 Test for ellipse

If the ellipse is given on complex plane by conformal map to unit circle interior problem, where then exact solution is . Eccentricity or ellipse is , and . We can get and by taking derivative:

Then apply our solver to find approximation of , numerical results when are shown below:

Fig. 2. (a) real part and imaginary part of exact solution versus . (b) absolute value of coefficient of or versus when number of collocation points (c) error of versus number of collocation points .

4.4 Test for overlapped circles ()

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 then the shape has two corners on y-axis.

Fig. 3(a) Overlapped circle shape when . Fig. 3(b) vesus when

The equation of this shape in polar coordinate is on first quadrant. Exact solution of is given by an integral equation:


where defined by:


where is center angle, is tension parallel to the x-axis and is tension parallel to the y-axis. The value of is evaluated by MATLAB numerical integration function ’integral’ (see Shampine (2008)). We test the case of longitudinal tension when and with . Notice as , from (36), in this case. Integral (34) contains highly oscilllating term which cause significant error when taking numerical integration. For this reason, we drop a small interval when showing error and evaluating error.

Fig. 3(c) absolute value of coefficient of or versus number of collocation points . (d) Error of vesus when .

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 when . The coefficients of is not decay so fast so the rate of convergence for overlapped circle with is slower (see Fig. 3(e))

Fig. 3(d) error on versus number of collocation points .

4.5 Test for overlapped circles ()

The shape of overlapped circles when 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 at corner.

Fig. 3(a) Overlapped circle shape when . Fig. 3(b) Numerical solution and exact solution of vesus when

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 and for large , the error comes from taking derivative of , truncation error for large . 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 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

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.


  • Boyd (2001) Boyd, J.P., 2001. Chebyshev and Fourier spectral methods. Courier Corporation.
  • 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.