Elasticity with Arbitrarily Shaped Inhomogeneity

TL;DR

This paper introduces a semi-analytic method using conformal maps to determine stress tensors in elastic materials with arbitrarily shaped inhomogeneities, improving understanding of complex inhomogeneity effects.

Contribution

The paper presents a novel semi-analytic approach employing dual conformal maps to solve elasticity problems with arbitrary inhomogeneity shapes, expanding analytical capabilities.

Findings

Method converges for mild shape distortions

Provides accurate stress tensor calculations

Compared successfully with previous results

Abstract

A classical problem in elasticity theory involves an inhomogeneity embedded in a material of given stress and shear moduli. The inhomogeneity is a region of arbitrary shape whose stress and shear moduli differ from those of the surrounding medium. In this paper we present a new, semi-analytic method for finding the stress tensor for an infinite plate with such an inhomogeneity. The solution involves two conformal maps, one from the inside and the second from the outside of the unit circle to the inside, and respectively outside, of the inhomogeneity. The method provides a solution by matching the conformal maps on the boundary between the inhomogeneity and the surrounding material. This matching converges well only for relatively mild distortions of the unit circle due to reasons which will be discussed in the article. We provide a comparison of the present result to known previous results.

Figures

Equations

$$\mathbf{u}(\mathbf{r})\equiv {\mathbf{r}}^{\prime }-\mathbf{r}.$$

$${ϵ}_{ij}\equiv [{\partial }_i{\mathbf{u}}_j+{\partial }_j{\mathbf{u}}_i]/2.$$

$${\sigma }_{ij}=2{\mu }_i{ϵ}_{ij}+{\lambda }_i{\delta }_{ij}{ϵ}_{kk}$$

$${\sigma }_{xx}=\frac{{\partial }^{2}U}{\partial y^2};{\sigma }_{xy}=-\frac{{\partial }^{2}U}{\partial x\partial y};{\sigma }_{yy}=\frac{{\partial }^{2}U}{\partial x^2},$$

$${\nabla }^2{\nabla }^{2}U(x,y)=0.$$

$$U(x,y)=\Re [\bar{z}\tilde{\phi }(z)+\eta (z)],$$

$${\sigma }_{yy}(x,y)$$

$${\sigma }_{xx}(x,y)$$

$${\sigma }_{xy}(x,y)$$

$$\tilde{\psi }(z)={\eta }^{\prime }(z)$$

$${\sigma }_{yy}(x,y)$$

$${\sigma }_{xx}(x,y)$$

$${\sigma }_{xy}(x,y)$$

$${\tilde{\phi }}_i$$

$${\tilde{\psi }}_i$$

$$\mathbf{\sigma }\cdot \text{n}={\sigma }_{xn}+i{\sigma }_{yn}$$

$${\sigma }_{xn}^{(1)}+i{\sigma }_{yn}^{(1)}={\sigma }_{xn}^{(2)}+i{\sigma }_{yn}^{(2)},$$

$$u_x^{(1)}+i{u}_y^{(1)}=u_x^{(2)}+i{u}_y^{(2)}.$$

$$\frac{d}{ds}\left(\frac{\partial U_1}{\partial x}+i\frac{\partial U_1}{\partial y}\right)=\frac{d}{ds}\left(\frac{\partial U_2}{\partial x}+i\frac{\partial U_2}{\partial y}\right)$$

$$\frac{\partial U_1}{\partial x}+i\frac{\partial U_1}{\partial y}=\frac{\partial U_2}{\partial x}+i\frac{\partial U_2}{\partial y}+\mathcal{C},$$

$${\tilde{\phi }}^{(1)}(z)+z\overline{{\tilde{\phi }}^{\prime (1)}(z)}+\overline{{\tilde{\psi }}^{(1)}(z)}={\tilde{\phi }}^{(2)}(z)+z\overline{{\tilde{\phi }}^{\prime (2)}(z)}+\overline{{\tilde{\psi }}^{(2)}(z)}+\mathcal{C},$$

$$\frac{[{\kappa }_1{\tilde{\phi }}^{(1)}(z)-z\overline{{\tilde{\phi }}^{\prime (1)}(z)}-\overline{{\tilde{\psi }}^{(1)}(z)}]}{{\mu }_1}=\frac{[{\kappa }_2{\tilde{\phi }}^{(2)}(z)-z\overline{{\tilde{\phi }}^{\prime (2)}(z)}-\overline{{\tilde{\psi }}^{(2)}(z)}]}{{\mu }_2}.$$

$${\sigma }_{xx}(\infty )=0,{\sigma }_{yy}(\infty )={\sigma }_{\infty },{\sigma }_{xy}(\infty )=0.$$

$${\tilde{\phi }}^{(1)}(z)$$

$${\tilde{\psi }}^{(1)}(z)$$

$${\tilde{\phi }}^{(2)}(z)$$

$${\tilde{\psi }}^{(2)}(z)$$

$${\tilde{\phi }}_1^{(2)}=\frac{{\sigma }_{\infty }}{4};{\tilde{\psi }}_1^{(2)}=\frac{{\sigma }_{\infty }}{2}.$$

$$z=\Phi (\omega )$$

$$z=\Lambda (\zeta )$$

$$\zeta ={\Lambda }^{-1}(z)$$

$$\omega ={\Phi }^{-1}(z).$$

$$\frac{1}{2\pi }{\int }_0^{2\pi }{e}^{ni\theta }{e}^{-mi\theta }={\delta }_{n,m}.$$

$${\tilde{\phi }}^{(2)}(z)\equiv {\phi }^{(2)}\left({\Phi }^{-1}(z)\right),\tilde{\psi }(z{)}^{(2)}\equiv {\psi }^{(2)}\left({\Phi }^{-1}(z)\right).$$

Full text

I Introduction

Elasticity theory in homogeneous materials is a well developed subject. Much less is known about inhomogeneous materials where the solution of the basic equations of elasticity becomes very involved. In this paper we focus on a material which consists of one finite area (of arbitrary shape) in which a material with given elastic properties is embedded in an infinite sheet of material of different elastic properties. This situation is known as “elastic inhomogeneity” and it appears in a variety of solid mechanical contexts. Previous studies have concentrated on solving this problem for the relatively symmetric case of an ellipse 54Har where it was solved analytically. In other cases the problem was solved for small perturbations to the circle 90Gao .

Mathematically the problem is set as follows, and see Fig. 1. A patch of material of type 1 occupies an area $\Omega$ and is delineated by a sharp boundary which will be denoted $\partial\Omega$. The rest of the infinite plane is made of material of type 2. The material is subjected to forces at infinity (and see below the precise boundary conditions), and is therefore deformed. Before the deformation each point of the material is assigned a point ${\bm{r}}$ in the two-dimensional plane. The forces at infinity result in a displacement of the material points to a new equilibrium position ${\bm{r}}^{\prime}$. The displacement field ${\bm{u}}({\bm{r}})$ is defined as

[TABLE]

The strain field is defined accordingly as

[TABLE]

In the context of linear elasticity in isotropic materials one then introduces the stress field according to Hook’s law

[TABLE]

Where $\lambda_{i}=\mu_{i}/(\frac{1}{2\nu_{i}}-1)$. $\mu_{i}$ and $\nu_{i}$ take on different values $\mu_{1}$, $\nu_{1}$ in $\Omega$ and $\mu_{2}$, $\nu_{2}$ in the rest of the material. In equilibrium the stress tensor should be divergenceless $\frac{\partial\sigma_{ij}}{\partial x_{j}}=0$ at each point in the sheet. By defining the stress (or Airy) potential U:

[TABLE]

the former equation for the stress tensor becomes a partial differential equation for the stress potential:

[TABLE]

This equation, which is known as the bi-Laplace or the bi-harmonic equation is conveniently solved as a non-analytic combination of analytic funcitons. To this aim we introduce the complex notation $z\equiv x+iy$, and note the general solutions of Eq. (5) in the form

[TABLE]

where $\bar{z}\equiv x-iy$ and $\tilde{\varphi}(z)$ and $\eta(z)$ are any two holomorphic functions. What remains to do in any particular problem is to find the unique analytic functions such that the stress tensor satisfies the boundary conditions. This stress tensor is determined by the two holomorphic functions as:

[TABLE]

We define for convenience:

[TABLE]

And then:

[TABLE]

Note that the stress tensor is determined by derivatives of the holomorphic functions, and not by the functions themselves. This leaves us with some freedom, since the functions can be changed with the following gauge:

[TABLE]

As we shall see below, not all these gauge freedoms are true freedoms once we introduce the boundary and continuity conditions. Since the elastic properties are different inside and outside $\Omega$, the potential functions will be different in the two regions: $\tilde{\varphi}_{1}$ and $\tilde{\psi}_{1}$ which are defined on $\Omega$ and $\tilde{\varphi}_{2}$ and $\tilde{\psi}_{2}$ which are defined on $\mathbb{C}\backslash\Omega$. Nevertheless we will demand continuity of the physical fields. In particular the normal force

[TABLE]

and the displacement $\textbf{u}(\textbf{r})$ must be continues across the interface (by Newton’s third law) in the absence of surface tension. Therefore the continuity conditions are:

[TABLE]

[TABLE]

The continuity conditions for the stress, can be rewritten as:

[TABLE]

or, after integrating:

[TABLE]

where ${\mathcal{C}}$ is a complex constant of integration. In terms of the holomorphic functions, the condition (12) translates to 53Mus :

[TABLE]

and the condition (13) becomes:

[TABLE]

where $\kappa_{i}=(3-\nu_{i})/(1+\nu_{i})$.

In addition to these continuity conditions on $\partial\Omega$ we need to specify boundary conditions at infinity. We choose

[TABLE]

II Solution in terms of conformal maps

Solutions to the problem of finding the stress field outside a given domain using conformal maps were described for example in Barra . Here we need to solve for the stress field both inside and outside the given domain. In the following we assume that the center of coordinates is inside $\Omega$ and the point at infinity is outside $\Omega$. Since the stress functions are holomorphic in their domains of definition, we can expand them in the appropriate laurent series, which for the functions with superscript (1) is of the form

[TABLE]

i.e. we have no poles at the origin. For the outside domain (functions with superscript (2)) the most general expansions in agreement with the boundary conditions (18) are of the form

[TABLE]

i.e we have a pole of order 1 at infinity. Accordingly, the leading terms of Eqs. (20) are determined by the boundary conditions. We now use one of the gauge freedoms to eliminate the imaginary part of $\varphi^{(2)}$ and write:

[TABLE]

The standard way to proceed 53Mus would be to substitute the series expansions in the continuity conditions and find the linear equations that determine all the coefficients by equating terms of the same order in $z$. However this cannot be done in general since the functions $z^{n}$ are not orthogonal on arbitrary contours $\partial\Omega$. To overcome this, one maps the regions $\Omega$ and $\mathbb{C}\backslash\Omega$ into the interior and exterior of the unit circle, respectively. That is, we need two holomorphic, invertible (and thus conformal) functions, one is

[TABLE]

which maps the exterior of the unit circle into $\mathbb{C}\backslash\Omega$, and the other is

[TABLE]

which maps the unit disk into $\Omega$. Since they are both invertible they have inverse functions which we denote

[TABLE]

and

[TABLE]

Now we express the functions $\tilde{\varphi}^{(i)}$ and $\tilde{\psi}^{(i)}$ in terms of $\omega$ and $\zeta$ and then expand them on the boundary of the unit circle. This expansion will be a Fourier series where the powers of $\omega$ or $\zeta$ satisfies the orthogonality relation:

[TABLE]

We have here used $e^{i\theta}$ to represent either $\omega$ or $\zeta$. The orthogonality allows us to equate the coefficients of the series term by term. We define:

[TABLE]

and:

[TABLE]

We can expand $\varphi_{i}$ and $\psi_{i}$ in terms of $\omega$ and $\zeta$ on the unit circle. Since the original functions were holomorphic and meromorphic in the original domains, the functions:

[TABLE]

and:

[TABLE]

are holomorphic inside and outside the unit disc, respectively. Therefore we can expand in terms of $\omega$ and $\zeta$:

[TABLE]

[TABLE]

We now assume that the map of the exterior domain, $\Phi$, maps the point at infinity to infinity. That is, it will have a Laurent series on the form

[TABLE]

From this we get the following relations (after substituting and taking the limit $\omega\rightarrow\infty$):

[TABLE]

we can also use the last two freedoms to choose $\tilde{\varphi}^{(2)}_{0}=-F_{0}\tilde{\varphi}^{(2)}_{1}$ such that $\varphi^{(2)}_{0}=0$. In the interior domain, the functions $\tilde{\varphi}^{(1)}$ and $\tilde{\psi}^{(1)}$ also have five freedom. However, the requirement of continuity of the displacement field ${\bm{u}}$ across the boundary $\partial\Omega$ removes three of these freedoms. This continuity was expressed by Eq. (13). Applying the apparent gauge freedoms on the LHS of that equation and then subtracting the resulting equation from Eq. (13) we find the three conditions

[TABLE]

Using the remaining two freedoms, we can eliminate the constant term in the expansion of $\psi^{(1)}$ by setting $\psi^{(1)}_{0}=0$. Note that this is possible only when we choose $\Lambda(\zeta)$ such that $\Lambda(0)=0$. We may always define our mapping $\Lambda$ such that this is satisfied. In terms of the conformal maps we transform the boundary conditions into

[TABLE]

[TABLE]

III Method of Solution

At this point we need to substitute the expansions (31),(32),(33) and an expansion similar to (33) for $\Lambda(\zeta)$ into the equations (36) and (37) and solve for the coefficients $\varphi^{(i)}_{k}$ and $\psi^{(i)}_{k}$. To understand how to do this in principle we write the expanded equations (31) and (32) in an abstract form

[TABLE]

where $p_{k}$ are linear combinations of the coefficients $\varphi^{(1)}_{n}$ and $\psi^{(1)}_{n}$ whereas $q_{m}$ are linear combinations of $\varphi^{(2)}_{n}$ and $\psi^{(2)}_{n}$. As this equation stands we cannot use the orthogonality relation Eq. (26). Therefore we expand moments of $\omega$ in terms of $\zeta$ in the form

[TABLE]

We now insert this expression in Eq. (38),

[TABLE]

and equate the coefficients of same powers to achieve a set of linear equations for the coefficients of $\varphi^{(i)}$ and $\psi^{(i)}$. The actual algebraic manipulations that are involved in reaching a finite set of linear equations are presented in the appendix.

Needless to say, when we derive a finite set of equations we lose precision. To see this we note that to get the right number of equations for the number of unknowns (see Appendix) we need to truncate the summations on the LHS and the RHS of Eq. (40) at the same finite N, i.e.

[TABLE]

For a precisely circular inclusion this truncation introduces no loss of information. For this particular shape the expansion Eq. (39) has only one term with $n=m$, i.e. $a_{n,m}=\delta_{n,m}$. Obviously, when the inclusion shape deviates from the circle, the representation of $\omega^{m}$ in Fourier space deviates from a delta function and it becomes more spread. An example of this phenomenon is presented in Fig.2 for an inclusion in the form of an ellipse with aspect ratio of about 1.5. The upper panel shows the paramterization of the outer mapping as function of that of the inner, $\arg(\Phi^{-1}(\Lambda(\omega)))$. In the lower panel we show the power spectrum $|a_{n,m}|^{2}$ of the moments $m=1$ and $m=25$ of the function in the upper panel. If we truncate the expansion at the dashed line $N=n=25$, we lose high frequency information for the higher moments. This loss of informaion will lead to stress field calculations which are less accurate.

To see the difficulty in a pictorial way we can consider the field lines of the conformal mappings for an inclusion that is elongated in shape, see for example Fig. 3. The external field lines concentrate at the convex parts of the inclusion, whereas the internal field lines concentrate on the concave parts. It becomes increasingly difficult to match field lines since they make a large discontinues jump when we go from the interior to the exterior domain. Similarly, for the ellipse in Fig. 2, when we increase the aspect ratio of the ellipse, the slope in the steep parts of $\arg(\Phi^{-1}(\Lambda(\omega)))$ become even larger, requiring higher order frequencies in our expansions. Eventually for large aspect ratios, our method will break down.

IV Obtaining the conformal maps

In all the calculations we assumed that the conformal maps $\Phi(\omega)$ and $\Lambda(\zeta)$ are available. For arbitrary inclusion shapes this is far from obvious, and special methods are necessary to obtain these maps. An efficient method to obtain the conformal map from the exterior of the unit circle to the exterior of an arbitrary given shape had been discussed in great detail in 06MPST . In the present case we use a slightly different method namely the geodesic algorithm 06MR . This method, like the former one, is based on the iterations of a generic conformal map $\gamma_{\mathbf{\xi}}$ defined by a set of parameters $\mathbf{\xi}$. We then construct the conformal map to an arbitrary shape by an appropriate choice of parameters $\mathbf{\xi}$. In the geodesic algorithm, we discretize the interface of the inclusion by a sequence of points $\{z_{k}\}_{k=0}^{n}$. The points appear sequentially in the positive direction of the interface.

We now briefly summarize how to construct the conformal map (see 06MR ). First we construct iteratively the inverse map that brings the interior of the inclusion to the upper half plane and the interface to the real axis. The conformal map to the shape then follows directly from the inverse. The construction is done in three steps. In the first step, we move one point to infinity and another to the center of coordinates, e.g. $z_{0}$ and $z_{1}$, respectively. For that purpose we use the mapping,

[TABLE]

In the next step we find a map that connects $z_{2}$ to the real axis by a semi circular arc. The inverse of this mapping, $\gamma_{\xi_{2}}$, brings $z_{2}$ to the real axis.

[TABLE]

where $\xi_{2}=\gamma_{1}(z_{2})$ and $a=|\xi_{2}|^{2}/\Re\xi_{2}$ and $b=|\xi_{2}|^{2}/\Im\xi_{2}$. Iteratively, we apply this mapping to all the points $z_{3},\ldots z_{n}$ where in general

[TABLE]

In the third and last step we unfold the remaining part of the interior to the whole upper half-plane by the map

[TABLE]

with

[TABLE]

The conformal map $\Psi$ from the upper half-plane to the interior domain and from the lower half-plane to the exterior domain is then given by

[TABLE]

From the conformal map $\Psi$ we easily construct the map $\Phi$ from the unit circle to the inclusion. In Fig. 6 we illustrate how this is done.

V Examples

In order to check the validity of our method we calculated the stress fields created by inhomogeneities with 2 different geometries: an ellipse with semi axes $0.9$ and $\sqrt{1+0.9^{2}}$ (aspect ratio of about 1.5) and a smoothed triangular curve $1.75z+\frac{0.2}{z^{2}}$ (see fig. 4). In the case of the elliptical inhomogeneity we compared our method to the known analytical solution which was first obtained by Hardiman 54Har in 1954. In the example below, the boundary conditions at infinity were set to $\sigma_{\infty}=-1\frac{N}{m^{2}}$. The shear moduli used were $\mu_{1}=1\frac{N}{m^{2}}$ for the inhomogeneity and $\mu_{2}=1.2\frac{N}{m^{2}}$ for the matrix. The Poisson ratio was taken to be $\nu=1$ for both inhomogeneity and matrix. In figs. 5 and 7 we can see the components of the stress field calculated outside the ellipse.

The blue line is the stress calculated using Hardiman’s solution and the red spots corresponds to the values obtained by our method. Similarly, we have calculated the stress field outside the triangular-like inhomogeneity (figs 8 and 9).

VI Concluding Remarks

In comparing our approach to other available algorithms, for example finite elements approximations to the equations of linear elasticity, we should stress that our approach works equally well for compressible and incompressible materials, There is no problem in taking the incompressible limit as the Poisson ratio approaches 1/2. This is not the case for finite elements methods. While the examples shown above worked out very well, indicating that the proposed algorithm is both elegant and numerically feasible, unfortunately it deteriorates very quickly when the shape of the inhomogeneity deviates strongly from circular symmetry. The difficulty in matching the two conformal maps is significant, as can be gleaned from from Figs. 2 and 3. One could think that the problem could be overcome in principle by increasing the numerical accuracy, but in practice, when the inhomogeneity has horns, spikes or deep fjords, the difficulties becomes insurmountable. Similar difficulties in another guise are however expected when any other analytic or semi-analytic method is used, leaving very contorted inhomogeneities as a remaining challenge for elasticity theory.

Acknowledgements.

This work had been supported in part by the Israel Science Foundation, the German Israeli Foundation and the Minerva Foundation. We thank Eran Bouchbinder, Felipe Barra, Anna Pomyalov and Charles Tresser for some very useful discussions.

Appendix A Explicit calculation

Starting from (36) and (37) we write the series expansions in a compact form:

[TABLE]

[TABLE]

were we have eliminated the zero terms in $\psi^{1}$ because of the gauge freedom. Substituting:

[TABLE]

[TABLE]

We can also expand:

[TABLE]

[TABLE]

Substituting and changing the indices of summation: $m=1-k+n$ which leads to $n=m+k-1$ we get :

[TABLE]

[TABLE]

In order to find the relations between the coefficients we need that both sides will be expressed in the same Fourier harmonics. Therefore, we need to expand:

[TABLE]

In a Fourier series. The condition for this to be expanded in a Fourier series is that $\omega(\zeta)$ is $L^{2}$ (i.e. $\int_{-\pi}^{\pi}|\omega(e^{\imath\theta})|^{2}d\theta<\infty$ on the segment $[-\pi,\pi]$). If that is the case, we can expand:

[TABLE]

Actually, it is found to be more convenient to expand powers of $\zeta$ in a Fourier series:

[TABLE]

were $\zeta=e^{\imath\theta}$ and $\int_{-\pi}^{\pi}|\omega(e^{\imath\theta})^{m}|^{2}d\theta<\infty$ on the segment $[-\pi,\pi]$.

[TABLE]

[TABLE]

Define:

[TABLE]

and

[TABLE]

Substituting:

[TABLE]

[TABLE]

Using the linear independence of the $\zeta^{n}$’s with respect to the Fourier integral, and ’cutting’ the infinite series at some number N, we get a set of linear equations which is of the form:

[TABLE]

Where $\mathbf{v}$ is the vector of coefficients ($\varphi$’s and $\psi$’s), $\hat{\mathbf{M}}$ is a $(4N+2)\times(2N+1)$ matrix of constants and $\mathbf{c}$ is the $C_{n}$’s.

Next, we substitute the expansions in the continuity equation for the displacement:

[TABLE]

Substituting:

[TABLE]

[TABLE]

We can also expand:

[TABLE]

[TABLE]

Substituting and changing the indices of summation: $m=1-k+n$ which leads to $n=m+k-1$:

[TABLE]

[TABLE]

Expanding $\omega(\zeta)$ as before and defining:

[TABLE]

we get:

[TABLE]

[TABLE]

When ’cutting’ the infinite series in the same way, we get again a matrix equation (with the same dimensions) of the form:

[TABLE]

Combining equations (51) and (58) we get a $(4N+2)\times(4N+1)$ matrix equation:

[TABLE]

where

[TABLE]

and

[TABLE]