Preprint
Article

This version is not peer-reviewed.

Theory and Applications of the Eshelby Ellipsoidal Elastic Inclusion Problem

Submitted:

07 December 2024

Posted:

11 December 2024

You are already at the latest version

Abstract
In this document, we fully review the theory and applications of the Eshelby Ellipsoidal Elastic Inclusion Problem. We rigorously derive all the equations related to the Eshelby Ellipsoidal Elastic Inclusion Problem and its applications to various Micro-mechanics problems like Ellipsodial Inhomogenity, Cracks, and Dislocations.
Keywords: 
;  ;  ;  ;  

1. Summary of Eshelby’s Inclusion and Inhomogeneity Problem

Eshelby’s Inclusion and Inhomogeneity problem is a fundamental concept in the field of continuum mechanics, particularly in the study of elastic fields and material science. It deals with the elastic behavior of a region within a homogeneous material (the matrix) that has different material properties or experiences different strains compared to the surrounding matrix. This problem is crucial in understanding how inhomogeneities such as voids, inclusions, or other defects within a material affect its overall mechanical properties.

1.1. Basic Definitions

  • Inclusion: An inclusion is a region within a material that has different elastic properties (stiffness, for example) from the surrounding material (the matrix). The inclusion is assumed to be embedded in the matrix and can have its own distinct material properties.
  • Inhomogeneity: An inhomogeneity refers to a region within a material where the material properties differ from the surrounding matrix. It is a broader term that includes inclusions but also refers to regions where properties such as density, thermal expansion, or other physical characteristics differ.

1.2. Eshelby’s Inclusion Problem

  • Eshelby’s Tensor: At the core of Eshelby’s inclusion problem is the Eshelby tensor, which describes the elastic field inside and around an inclusion when it is subjected to an external stress or strain. This tensor is a fourth-order tensor that relates the applied strain to the strain inside the inclusion.
  • Ellipsoidal Inclusions: Eshelby’s work showed that for ellipsoidal inclusions, the strain inside the inclusion is uniform and can be related to the external strain through the Eshelby tensor. This remarkable result simplifies the analysis of inclusions significantly, as it reduces the complexity of the problem.
  • Inclusion vs. Matrix: The key idea is that when an inclusion is subjected to a uniform external stress or strain, the strain field inside the inclusion remains uniform, although different from the strain field in the surrounding matrix. The specific relationship between these strains is governed by the shape of the inclusion and the Eshelby tensor.

1.3. Mathematical Formulation

  • Eigenstrain: The concept of eigenstrain (or transformation strain) is central to Eshelby’s analysis. Eigenstrain refers to a strain that would exist in the inclusion if it were isolated from the matrix and allowed to undergo a strain freely. When the inclusion is embedded in the matrix, the surrounding material restricts this strain, leading to an interaction between the inclusion and the matrix.
  • Elastic Field Equations: The elastic field due to an inclusion is governed by the equations of elasticity. For a linear elastic material, these equations are linear partial differential equations (PDEs) involving the stress and strain fields, which are solved subject to boundary conditions at the inclusion-matrix interface.
  • Eshelby’s Solution: Eshelby provided an analytical solution for the elastic field both inside and outside an ellipsoidal inclusion in an infinite medium. His solution showed that the strain inside the inclusion is constant and can be calculated using the Eshelby tensor.

1.4. Inhomogeneity Problem

  • Difference from Inclusion: In the case of an inhomogeneity, the material properties of the region differ from those of the matrix, leading to a more complex interaction between the region and the surrounding material. Unlike an inclusion, where the material inside the inclusion can be imagined as having the same properties as the matrix, an inhomogeneity represents a real difference in material properties.
  • Complexity: The solution to the inhomogeneity problem is more complex than the inclusion problem because the contrast in material properties must be accounted for. This typically requires solving the elasticity equations with variable material coefficients.
  • Perturbation Techniques: In practice, solutions to inhomogeneity problems often involve perturbation techniques, where the problem is treated as a small deviation from the homogeneous case, or numerical methods, where the equations are solved using computational techniques.

1.5. Applications

  • Material Science: Eshelby’s inclusion theory is widely used in materials science to predict how inclusions and inhomogeneities affect the mechanical properties of composites, polycrystals, and other heterogeneous materials.
  • Micromechanics: The theory forms the basis for many micromechanical models that predict the behavior of materials with microstructural features such as voids, fibers, or precipitates.
  • Fracture Mechanics: In fracture mechanics, Eshelby’s theory is used to understand how cracks and other defects influence the stress distribution in materials, which is crucial for predicting failure.

1.6. Extensions and Generalizations

  • Non-Ellipsoidal Inclusions: While Eshelby’s original work focused on ellipsoidal inclusions, subsequent research has extended the theory to non-ellipsoidal shapes, though these cases generally require numerical solutions or approximations.
  • Anisotropic Materials: The theory has also been extended to anisotropic materials, where the material properties differ in different directions, adding another layer of complexity to the problem.

1.7. Limitations and Challenges

  • Finite Boundaries: Eshelby’s solution assumes an infinite medium, which is an idealization. In real-world applications, the finite size of the material can influence the stress and strain fields, requiring corrections or alternative methods.
  • Nonlinear Materials: The theory is based on linear elasticity, and its application to nonlinear materials is limited. In such cases, more advanced models are needed.
In summary, Eshelby’s Inclusion and Inhomogeneity problem provides a powerful framework for understanding how embedded regions within a material interact with their surroundings and affect the material’s overall properties. The theory’s simplicity and analytical nature make it a cornerstone of material science, despite the challenges in extending it to more complex scenarios.
We start with the derivation of the 2 double derivatives d x j ( d x i ( r ) ) and d x i d x j d x k d x l r 3 , where r is the radius in Cartesian coordinates r = x 1 2 + x 2 2 + x 3 2 and their integration over the surface S of the inclusion. These double derivatives and their integration over the surface S of the inclusion will be used to find the displacement impressed on the material in stage III (using the Love 1927 equation (shown in Section 6) of Displacement at r due to point force F i at r ) due to the application of the force distribution F j = p j k T n k over S to make the body free of external force (but in a state of self-stress because of the transformation of the inclusion).

2. Derivation of Double Derivative and Their Integration over the Surface §

In this section, we compute the double derivative d x j ( d x i ( r ) ) where r is the radius in Cartesian coordinates r = x 1 2 + x 2 2 + x 3 2 and their integration over the surface S. We start with computing d x j ( d x i ( r ) ) :
Step 1: The first derivative of the radius with respect to x i is:
d x i ( r ) = r x i = x i r
Step 2: The second derivative is:
d x j ( d x i ( r ) ) = x j x i r
Applying the product rule:
d x j ( d x i ( r ) ) = δ i j r x i x j r 3
where δ i j is the Kronecker delta. Let’s now integrate the above double derivative multiplied with the vector n j over the surface. Before doing note the following important identity for the solid angle ω
d ω = n i r i d S r 3 = n i l i d S r 2
n i l i d S = r 2 d ω
Therefore the integration of the above double derivative d x j ( d x i ( r ) ) multiplied with the vector n j over the surface S shall be
S d x j ( d x i ( r ) ) n j d S = S ( δ i j r x i x j r 3 ) n j d S
S d x j ( d x i ( r ) ) n j d S = S ( δ i j r ) n j d S S ( l i l j r ) n j d S
Substituting the relations n j l j d S = r 2 d ω in the above boxed equation
S d x j ( d x i ( r ) ) n j d S = S ( δ i j n j r ) d S 0 4 π ( l i r ) r 2 d ω
S d x j ( d x i ( r ) ) n j d S = S ( δ i j n j r ) d S 0 4 π l i r d ω
S d x j ( d x i ( r ) ) n j d S = S ( δ i j n j r ) d S 0 4 π r ( l ) l i d ω ( l )
Let’s define the tensor F in indicial notation as follows
F i j = δ i j r
Using the Gauss Divergence theorem we can therefore say that
S ( δ i j n j r ) d S = S ( F . n ) d S = V ( . F ) d V = V F i j x j d V
Now note that
V F i j x j d V = V x j ( δ i j r ) d V = V δ i j x j ( 1 r ) d V = V δ i j ( x j r 3 ) d V = V ( x i r 3 ) d V
S ( δ i j n j r ) d S = V ( x i r 3 ) d V
Let us compute the volume integral by integrating over an elementary cone d Ω centred on the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) with its vertex at x. The volume d V of this elementary cone is
d V = r 2 d r d ω
Therefore the integral V ( x i r 3 ) d V can be therefore written as
V ( x i r 3 ) d V = 0 4 π 0 r x i r d r d ω = 0 4 π 0 r l i d r d ω = 0 4 π l i 0 r d r d ω = 0 4 π r ( l ) l i d ω ( l )
S ( δ i j n j r ) d S = 0 4 π r ( l ) l i d ω ( l )
Therefore we have
S d x j ( d x i ( r ) ) n j d S = S ( δ i j n j r ) d S 0 4 π l i r d ω
S d x j ( d x i ( r ) ) n j d S = 0 4 π r ( l ) l i d ω ( l ) 0 4 π r ( l ) l i d ω ( l )
S d x j ( d x i ( r ) ) n j d S = 2 0 4 π r ( l ) l i d ω ( l )

3. Derivation of Quadruple Derivative and Their Integration over the Surface S

In this section, we do the derivation of the double derivative d x i d x j d x k d x l r 3 , where r is the radius in Cartesian coordinates r = x 1 2 + x 2 2 + x 3 2 and their integration over the surface S using 3 different methods.

3.1. Method 1: With Using the Exchanging of Integration and Differentiation Operator

We now derive d x i d x j d x k d x l r 3 , where r = x 1 2 + x 2 2 + x 3 2 . Therefore we have
r 3 = ( x 1 2 + x 2 2 + x 3 2 ) 3 / 2
Step 1: The first derivative of r 3 with respect to x l is:
d x l r 3 = r 3 x l = 3 r 2 x l r = 3 r x l
Step 2: The second derivative is:
d x k 3 r x l = 3 r x k x l + r x l x k
= 3 x k r x l + r δ l k
= 3 x k x l r + r δ l k
Therefore the integration of the second derivative shall be
S d x k ( d x l ( r 3 ) ) n l d S = S 3 x k x l r + r δ l k n l d S = S 3 x k l l + r δ l k n l d S
S d x k ( d x l ( r 3 ) ) n l d S = 3 S ( l l x k ) n l d S + 3 S ( r δ l k ) n l d S
S d x k ( d x l ( r 3 ) ) n l d S = 3 V x l ( l l x k ) d V + 3 V x l ( r δ l k ) d V = 3 V x l ( l l x k ) d V + 3 V x k ( r ) d V
S d x k ( d x l ( r 3 ) ) n l d S = 3 V x l ( l l x k ) d V + 3 V l k d V
Now note that we have
x l ( l l x k ) = l l x l ( x k ) + x k x l ( l l ) = l l x l ( x k ) + x k x l ( x l r )
x l ( l l x k ) = l l x l ( x k ) + x k ( r x l ( x l ) x l x l r ) r 2 = l l x l ( x k ) + x k ( 3 r x l x l r ) r 2 = l l δ k l + x k ( 3 r r ) r 2
x l ( l l x k ) = 3 l k
Therefore we have
S d x k ( d x l ( r 3 ) ) n l d S = 12 V l k d V
Differentiating both sides with respect to x i and x j we get
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 2 x i x j V l k d V
Integration ∫ and differentiation 2 x i x j can be exchanged since they are linear operators
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 V 2 x i x j ( l k ) d V
The differentiation 2 x i x j ( l k ) can be written as
2 x i x j ( l k ) = 2 x i x j ( x k r ) = x i ( x j ( x k r ) ) = x i ( ( r x j x k x k x j r ) r 2 )
2 x i x j ( l k ) = x i ( ( r δ j k x k x j r ) r 2 ) = x i ( δ j k r x k x j r 3 )
2 x i x j ( l k ) = δ j k x i ( 1 r ) x i ( x k x j r 3 ) = δ j k x i r 3 1 r 3 ( x k x i x j + x j x i x k ) x k x j x i ( 1 r 3 )
2 x i x j ( l k ) = δ j k x i r 3 1 r 3 ( x k δ j i + x j δ k i ) x k x j x i ( 1 r 3 )
2 x i x j ( l k ) = δ j k x i r 3 1 r 3 ( x k δ j i + x j δ k i ) + 3 x k x j x i r 5
2 x i x j ( l k ) = 1 r 3 ( x i δ j k + x k δ j i + x j δ k i ) + 3 x k x j x i r 5
2 x i x j ( l k ) = 1 r 2 ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i )
Therefore we have
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 V 1 r 2 ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d V
Let us compute the volume integral by integrating over an elementary cone d Ω centred on the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) with its vertex at x. The volume d V of this elementary cone is
d V = r 2 d r d ω
Therefore we have
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π 0 r 1 r 2 ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) r 2 d r d ω
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) 0 r d r d ω
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π r ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω

3.2. Method 2: With Using Exchanging of Integration and Differentiation Operator

We now derive d x i d x j d x k d x l r 3 , where r = x 1 2 + x 2 2 + x 3 2 . Therefore we have
r 3 = ( x 1 2 + x 2 2 + x 3 2 ) 3 / 2
Step 1: The first derivative of r 3 with respect to x l is:
d x l r 3 = r 3 x l = 3 r 2 x l r = 3 r x l
Step 2: The second derivative is:
d x k 3 r x l = 3 r x k x l + r x l x k
= 3 x k r x l + r δ l k
= 3 x k x l r + r δ l k
Therefore the integration of the second derivative shall be
S d x k ( d x l ( r 3 ) ) n l d S = S 3 x k x l r + r δ l k n l d S = S 3 x k l l + r δ l k n l d S
S d x k ( d x l ( r 3 ) ) n l d S = 3 S ( l l n l ) x k d S + 3 S ( r δ l k ) n l d S
We know that
l l n l d S = r ( l ) 2 d ω
Here r ( l ) is the distance between the point inside the volume V to the surface d S in the direction of l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) . Therefore we have
S d x k ( d x l ( r 3 ) ) n l d S = 3 0 4 π r ( l ) 2 x k d ω + 3 S ( r ( l ) n k ) d ω
Differentiating both sides with respect to x i and x j we get
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 2 x i x j ( 0 4 π r ( l ) 2 x k d ω ) + 3 2 x i x j ( S r ( l ) n k d S )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 2 x i x j ( 0 4 π r 2 x k d ω ) + 3 2 x i x j ( V ( x k r ) d V )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 2 x i x j ( 0 4 π r 2 x k d ω ) + 3 2 x i x j ( V ( x k r ) d V )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 2 x i x j ( 0 4 π r 2 x k d ω ) + 3 2 x i x j ( V l k d V )
Let us compute the volume integral by integrating over an elementary cone d Ω centred on the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) with its vertex at x. The volume d V of this elementary cone is
d V = r 2 d r d ω
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 2 x i x j ( 0 4 π r 2 x k d ω ) + 3 2 x i x j ( 0 4 π 0 r ( l k r 2 ) d r d ω )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 2 x i x j ( 0 4 π r 2 x k d ω ) + 3 2 x i x j ( 0 4 π l k 0 r r 2 d r d ω )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 2 x i x j ( 0 4 π r 2 x k d ω ) + 2 x i x j ( 0 4 π l k r 3 d r d ω )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π r 2 x k d ω )
We shall now show that the Integration operator ∫ and Differentiation operator 2 x i x j can be exchanged if and only if when the measure of the integration variable is independent of the Differentiation operator variables x i , x j . We shall explain this concept using 6 different cases which are as follows
  • Assuming that the incremental solid angle measure d ω is independent of the Differentiation operator variables x i , x j .
  • Assuming that the incremental volume measure d V = r ( l ) 3 3 d ω (subtended by a cone emanating from a point inside the volume V to the surface d S which is at a distance r ( l ) in the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) from the point) is independent of the Differentiation operator variables x i , x j .
  • Assuming that the incremental surface area measure d S = r ( l ) 2 d ω (subtended by a cone emanating from a point inside the volume V to the surface d S which is at a distance r ( l ) in the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) from the point) is independent of the Differentiation operator variables x i , x j .
  • Assuming that the incremental volume measure d V = r 2 d r d ω is independent of the Differentiation operator variables x i , x j .
  • Assuming that the incremental angular direction measure d ω and incremental volume measure d V is independent of the Differentiation operator variables x i , x j in the first integration and second integration respectively.
  • Assuming that the incremental angular direction measure l k d ω is independent of the Differentiation operator variables x i , x j .
We shall prove that all the above-mentioned 6 cases except the Case 4: "Assuming that the incremental volume measure d V = r 2 d r d ω
is independent of the Differentiation operator variables x i , x j " will lead to a wrong answer. Only Case 4: "Assuming that the incremental
volume measure d V = r 2 d r d ω is independent of the Differentiation operator variables x i , x j " will lead to the Correct answer as derived in Method 1.

3.2.1. Assuming that the Incremental Solid Angle Measure d ω is Constant with Respect to the Differentiation Operator Variables x i , x j

In this subsection, we shall show that Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged when we assume that the incremental solid angle measure d ω is independent of the Differentiation operator variables x i , x j . Note that we have earlier derived
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π r 2 x k d ω )
Let’s see what happens when we interchange Integration operator ∫ and Differentiation operator 2 x i x j in the RHS of the above equation
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π r 2 x k d ω ) = 4 ( 0 4 π 2 x i x j ( r 2 x k ) d ω )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 ( 0 4 π 2 x i x j ( r 2 x k ) d ω )
The first differentiation 2 x i x j ( r 2 x k ) can be written as
2 x i x j ( r 2 x k ) = x i ( x j ( r 2 x k ) ) = x i ( x k x j ( r 2 ) + r 2 x j ( x k ) ) = x i ( 2 x j x k + r 2 δ j k )
2 x i x j ( r 2 x k ) = x i ( 2 x j x k + r 2 δ j k ) = 2 x j x i x k + 2 x k x i x j + x i ( r 2 δ j k )
2 x i x j ( r 2 x k ) = 2 x j δ k i + 2 x k δ j i + δ j k x i ( r 2 ) = 2 x j δ k i + 2 x k δ j i + 2 x i δ j k
2 x i x j ( r 2 x k ) = 2 x j δ k i + 2 x k δ j i + 2 x i δ j k = 2 r ( l j δ k i + l k δ j i + l i δ j k )
Therefore we can write
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 ( 0 4 π 2 x i x j ( r 2 x k ) d ω ) = 8 S r ( l j δ k i + l k δ j i + l i δ j k ) d ω
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 8 S r ( l j δ k i + l k δ j i + l i δ j k ) d ω
As we discussed earlier in Section 3.1, the correct answer for the 2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) shall be
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( ( x k δ i j + x j δ i k + x i δ k j ) + 3 x k x j x i r 2 ) d ω
Therefore exchanging the Integration operator ∫ and Differentiation operator 2 x i x j under the assumption that the incremental solid angle measure d ω is independent of the Differentiation operator variables x i , x j leads to wrong answer. Hence Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged under the assumption that the incremental solid angle measure d ω is independent of the Differentiation operator variables x i , x j .

3.2.2. Assuming that the Incremental Volume Measure d V = r 3 3 d ω is Independent of the Differentiation Operator Variables x i , x j

In this subsection, we shall show that Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged when we assume that the incremental volume measure d V = r 3 3 d ω (subtended by a cone emanating from a point inside the volume V to the surface d S which is at a distance r from the point) is independent of the Differentiation operator variables x i , x j . Note that we have earlier derived
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π r 2 x k d ω )
Let’s see what happens when we interchange Integration operator ∫ and Differentiation operator 2 x i x j in the RHS of the above equation
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π r 2 x k d ω ) = 4 0 4 π 2 x i x j ( r 2 x k d ω )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 0 4 π 2 x i x j ( r 2 x k d ω )
Since we have assumed that d V = r 3 3 d ω is independent of the Differentiation operator variables x i , x j , the differentiation with respect to x i shall be zero i.e. x i d V = 0 , therefore we can write
x i d V = x i ( r 3 3 d ω ) = x i ( r 3 3 ) d ω + r 3 3 x i d ω
x i d V = r ( r 3 3 ) x i ( r ) d ω + r 3 3 x i d ω
x i d V = r 2 ( x i r ) d ω + r 3 3 x i ( d ω )
x i d V = r x i d ω + r 3 3 x i ( d ω )
The incremental volume d V remains constant irrespective change of x i , i.e x i d V = 0 therefore we have
0 = r x i d ω + r 3 3 x i ( d ω )
( 3 r x i r 3 ) d ω = x i ( d ω )
x i ( d ω ) = ( 3 x i r 2 ) d ω
This shows that the solid angle integration measure d ω is not independent of Differentiation operator variables x i , x j when we assume that the incremental volume measure d V = r 3 3 d ω is independent of the Differentiation operator variables x i , x j . The differentiation 2 x i x j ( r 2 x k d ω ) within the integral 0 4 π 2 x i x j ( r 2 x k d ω ) can be written as
2 x i x j ( r 2 x k d ω ) = x i ( x j ( r 2 x k d ω ) ) = x i ( x k x j ( r 2 d ω ) + r 2 d ω x j ( x k ) )
2 x i x j ( r 2 x k d ω ) = x i ( x k ( r 2 x j ( d ω ) + d ω x j ( r 2 ) ) + r 2 d ω x j ( x k ) )
2 x i x j ( r 2 x k d ω ) = x i ( x k r 2 x j ( d ω ) + x k d ω x j ( r 2 ) + r 2 d ω x j ( x k ) )
2 x i x j ( r 2 x k d ω ) = x i ( x k r 2 x j ( d ω ) + 2 x k x j d ω + r 2 d ω δ k j )
Substituting the equation x j ( d ω ) = 3 x j r 2 d ω which we derived earlier into the above boxed equation we get
2 x i x j ( r 2 x k d ω ) = x i ( 3 x j x k d ω + 2 x k x j d ω + r 2 d ω δ k j )
2 x i x j ( r 2 x k d ω ) = x i ( x j x k d ω + r 2 d ω δ k j ) = x i ( x j x k d ω ) + δ k j x i ( r 2 d ω )
2 x i x j ( r 2 x k d ω ) = ( x k d ω x i x j x j d ω x i x k x k x j x i d ω ) + δ k j ( d ω x i ( r 2 ) + r 2 x i ( d ω ) )
2 x i x j ( r 2 x k d ω ) = ( x k d ω x i x j x j d ω x i x k x k x j x i d ω ) + δ k j ( 2 x i d ω + r 2 x i ( d ω ) )
Substituting the equation x i ( d ω ) = 3 x i r 2 d ω which we derived earlier into the above boxed equation we get
2 x i x j ( r 2 x k d ω ) = ( x k δ i j d ω x j δ i k d ω + 3 x k x j x i r 2 d ω ) + δ k j ( 2 x i d ω 3 x i d ω )
2 x i x j ( r 2 x k d ω ) = ( x k δ i j d ω x j δ i k d ω + 3 x k x j x i r 2 d ω ) x i δ k j d ω
2 x i x j ( r 2 x k d ω ) = ( x k δ i j + x j δ i k + x i δ k j ) d ω + 3 x k x j x i r 2 d ω
Therefore we can write
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 ( 0 4 π 2 x i x j ( r 2 x k d ω ) ) = 4 0 4 π ( ( x k δ i j + x j δ i k + x i δ k j ) + 3 x k x j x i r 2 ) d ω
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 0 4 π ( ( x k δ i j + x j δ i k + x i δ k j ) + 3 x k x j x i r 2 ) d ω
As we discussed earlier in Section 3.1, the correct answer for the 2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) shall be
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( ( x k δ i j + x j δ i k + x i δ k j ) + 3 x k x j x i r 2 ) d ω
Therefore exchanging the Integration operator ∫ and Differentiation operator 2 x i x j under the assumption that the incremental volume measure d V = r 3 3 d ω (subtended by a cone emanating from a point inside the volume V to the surface d S which is at a distance r from the point) is independent of the Differentiation operator variables x i , x j leads to wrong answer. Hence Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged under the assumption that the incremental volume measure d V = r 3 3 d ω .

3.2.3. Assuming that the Incremental Volume Measure d V = r 2 d r d ω is Independent of the Differentiation Operator Variables x i , x j

In this subsection, we shall show that Integration operator ∫ and Differentiation operator 2 x i x j can be exchanged when we assume that the incremental volume measure d V = r 2 d r d ω is independent of the Differentiation operator variables x i , x j . Note that here r is not the distance between the point inside the volume V and surface d S . We have earlier derived that
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π r ( l ) 2 x k d ω )
This can be alternatively written as
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π l k r ( l ) 3 d ω ) = 4 2 x i x j ( 0 4 π 3 l k 0 r ( l ) r 2 d r d ω ) = 12 2 x i x j ( 0 4 π l k 0 r ( l ) r 2 d r d ω )
Since the direction cosine l k remains constant as we integrate from 0 to r ( l ) , we can take the direction cosine l k inside the second integral 0 r ( l ) . Therefore we can write the above equation as
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 2 x i x j ( 0 4 π 0 r ( l ) l k r 2 d r d ω )
Let’s see what happens when we interchange Integration operator ∫ and Differentiation operator 2 x i x j in the RHS of the above equation
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 2 x i x j ( 0 4 π 0 r ( l ) l k r 2 d r d ω ) = 12 0 4 π 0 r ( l ) 2 x i x j ( l k r 2 d r d ω )
Since r 2 d r d ω is independent of the Differentiation operator variables x i , x j , we have 2 x i x j ( r 2 d r d ω ) = 0 . We can therefore write
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π 0 r ( l ) ( 2 x i x j l k ) r 2 d r d ω
The differentiation 2 x i x j ( l k ) can be written as
2 x i x j ( l k ) = 2 x i x j ( x k r ) = x i ( x j ( x k r ) ) = x i ( ( r x j x k x k x j r ) r 2 )
2 x i x j ( l k ) = x i ( ( r δ j k x k x j r ) r 2 ) = x i ( δ j k r x k x j r 3 )
2 x i x j ( l k ) = δ j k x i ( 1 r ) x i ( x k x j r 3 ) = δ j k x i r 3 1 r 3 ( x k x i x j + x j x i x k ) x k x j x i ( 1 r 3 )
2 x i x j ( l k ) = δ j k x i r 3 1 r 3 ( x k δ j i + x j δ k i ) x k x j x i ( 1 r 3 )
2 x i x j ( l k ) = δ j k x i r 3 1 r 3 ( x k δ j i + x j δ k i ) + 3 x k x j x i r 5
2 x i x j ( l k ) = 1 r 3 ( x i δ j k + x k δ j i + x j δ k i ) + 3 x k x j x i r 5
2 x i x j ( l k ) = 1 r 2 ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i )
Therefore we have
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 V 1 r 2 ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d V
Let us compute the volume integral by integrating over an elementary cone d Ω centred on the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) with its vertex at x. The volume d V of this elementary cone is
d V = r 2 d r d ω
Therefore we have
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π 0 r 1 r 2 ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) r 2 d r d ω
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) 0 r d r d ω
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π r ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω
As we discussed earlier in Section 3.1, the correct answer for the 2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) shall be
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( ( x k δ i j + x j δ i k + x i δ k j ) + 3 x k x j x i r 2 ) d ω
Therefore exchanging the Integration operator ∫ and Differentiation operator 2 x i x j under the assumption that the incremental volume measure d V = r 2 d r d ω is independent of the Differentiation operator variables x i , x j leads to correct answer. Hence Integration operator ∫ and Differentiation operator 2 x i x j can be exchanged under the assumption that the incremental volume measure d V = r 2 d r d ω is independent of the Differentiation operator variables x i , x j .
Reason why we are getting different answers
The reason why we get a different answer than when we consider the incremental volume measure d V = r 3 3 d ω to be independent of the Differentiation operator variables x i , x j is due to the following reason: If we consider incremental volume measure d V = r 2 d r d ω to be independent of the Differentiation operator variables x i , x j , i.e x i ( r 2 d r d ω ) shall be zero. Then we get
x i d V = x i ( r 2 d r d ω ) = r 2 d r x i ( d ω ) + r 2 d ω x i ( d r ) + d ω d r x i ( r 2 )
x i d V = r 2 d r x i ( d ω ) + r 2 d ω x i ( d r ) + d ω d r x i ( r 2 )
x i d V = r 2 d r x i ( d ω ) + 2 x i d ω d r
The incremental volume d V remains constant irrespective change of x i , i.e x i d V = 0 therefore we have
0 = r 2 d r x i ( d ω ) + 2 x i d ω d r
2 x i d ω = r 2 x i ( d ω )
x i ( d ω ) = 2 x i r 2 d ω
But if we consider the incremental volume measure d V = r 3 3 d ω to be independent of the Differentiation operator variables x i , x j , i.e x i ( r 3 3 d ω ) shall be zero. Then we get
x i d V = x i ( r 3 3 d ω ) = x i ( r 3 3 ) d ω + r 3 3 x i d ω
x i d V = r ( r 3 3 ) x i ( r ) d ω + r 3 3 x i d ω
x i d V = r 2 ( x i r ) d ω + r 3 3 x i ( d ω )
x i d V = r x i d ω + r 3 3 x i ( d ω )
The incremental volume d V remains constant irrespective change of x i , i.e x i d V = 0 therefore we have
0 = r x i d ω + r 3 3 x i ( d ω )
( 3 r x i r 3 ) d ω = x i ( d ω )
x i ( d ω ) = ( 3 x i r 2 ) d ω
The x i ( d ω ) is different for different assumptions of integration measure independence with respect to the Differentiation operator variables x i , x j . This leads to different values of 2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) .

3.2.4. Assuming that the Incremental Angular Direction Measure l k d ω is Independent of the Differentiation Operator Variables x i , x j

In this subsection, we shall show that Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged when we assume that the incremental angular direction measure l k d ω is independent of the Differentiation operator variables x i , x j . Note that here r is not the distance between the point inside the volume V and surface d S . We have earlier derived that
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π r ( l ) 2 x k d ω )
Let’s see what happens when we interchange Integration operator ∫ and Differentiation operator 2 x i x j in the RHS of the above equation
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 2 x i x j ( 0 4 π r ( l ) 2 x k d ω ) = 4 0 4 π 2 x i x j ( r ( l ) 2 x k d ω )
Since we have assumed that the incremental angular direction measure l k d ω is independent of the Differentiation operator variables x i , x j , the differentiation with respect to x i shall be zero i.e. x i ( l k d ω ) = 0 , therefore we can write
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 ( 0 4 π 2 x i x j ( r ( l ) 2 x k ) d ω ) = 4 ( 0 4 π 2 x i x j ( r ( l ) 3 l k ) d ω ) = 4 ( 0 4 π l k 2 x i x j ( r ( l ) 3 ) d ω )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 4 ( 0 4 π l k 2 x i x j ( r ( l ) 3 ) d ω )
We know that
2 x i x j ( r 3 ) = 3 x i x j r + r δ i j = 3 r l i l j + r δ i j
2 x i x j ( r 3 ) = 3 r ( l i l j + δ i j )
Therefore can write
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 ( 0 4 π r ( l i l j l k + δ i j l k ) d ω )
As we discussed earlier in Section 3.1, the correct answer for the 2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) shall be
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( ( x k δ i j + x j δ i k + x i δ k j ) + 3 x k x j x i r 2 ) d ω
Therefore exchanging the Integration operator ∫ and Differentiation operator 2 x i x j under the assumption that the incremental angular direction measure l k d ω is independent of the Differentiation operator variables x i , x j leads to wrong answer. Hence Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged under the assumption that the incremental angular direction measure l k d ω is independent of the Differentiation operator variables x i , x j .

3.2.5. Assuming that the Incremental Angular Direction Measure d ω and Incremental Volume Measure d V is Independent of the Differentiation Operator Variables x i , x j in the First Integration and Second Integration Respectively

In this subsection, we shall show that Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged when we assume that the incremental angular direction measure d ω and incremental volume measure d V is independent of the Differentiation operator variables x i , x j in the first integration and second integration respectively. Note that here r is not the distance between the point inside the volume V and surface d S . We have earlier derived that
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 2 x i x j ( 0 4 π r 2 x k d ω ) + 3 2 x i x j ( V l k d V )
Let’s see what happens when we interchange Integration operator ∫ and Differentiation operator 2 x i x j in the RHS of the above equation
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 0 4 π 2 x i x j ( r 2 x k d ω ) + 3 V 2 x i x j ( l k d V )
Since the incremental angular direction measure d ω and incremental volume measure d V is independent of the Differentiation operator variables x i , x j in the first integration and second integration respectively, we can write the above equation as
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 3 0 4 π 2 x i x j ( r 2 x k ) d ω + 3 V 2 x i x j ( l k ) d V
The first differentiation 2 x i x j ( r 2 x k ) can be written as
2 x i x j ( r 2 x k ) = x i ( x k x j ( r 2 ) + r 2 x j ( x k ) ) = x i ( 2 x k x j + r 2 δ k j )
2 x i x j ( r 2 x k ) = x i ( 2 x k x j + r 2 δ k j ) = 2 ( x k x i x j + x j x i x k ) + δ k j x i r 2
2 x i x j ( r 2 x k ) = x i ( 2 x k x j + r 2 δ k j ) = 2 ( x k δ i j + x j δ i k ) + 2 x i δ k j
2 x i x j ( r 2 x k ) = 2 ( x k δ i j + x j δ i k + x i δ k j )
2 x i x j ( r 2 x k ) = 2 r ( l k δ i j + l j δ i k + l i δ k j )
The second differentiation 2 x i x j ( x k r ) can be written as
2 x i x j ( x k r ) = x i ( x j ( x k r ) ) = x i ( ( r x j x k x k x j r ) r 2 )
2 x i x j ( x k r ) = x i ( ( r δ j k x k x j r ) r 2 ) = x i ( δ j k r x k x j r 3 )
2 x i x j ( x k r ) = δ j k x i ( 1 r ) x i ( x k x j r 3 ) = δ j k x i r 3 1 r 3 ( x k x i x j + x j x i x k ) x k x j x i ( 1 r 3 )
2 x i x j ( x k r ) = δ j k x i r 3 1 r 3 ( x k δ j i + x j δ k i ) x k x j x i ( 1 r 3 )
2 x i x j ( x k r ) = δ j k x i r 3 1 r 3 ( x k δ j i + x j δ k i ) + 3 x k x j x i r 5
2 x i x j ( x k r ) = 1 r 3 ( x i δ j k + x k δ j i + x j δ k i ) + 3 x k x j x i r 5
2 x i x j ( x k r ) = 1 r 2 ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i )
Note that we have d V = r 3 3 d ω . Therefore can write
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 6 0 4 π r ( l k δ i j + l j δ i k + l i δ k j ) d ω + 3 V ( 1 r 2 ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) ) d V
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 6 0 4 π r ( l k δ i j + l j δ i k + l i δ k j ) d ω + V ( r ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) ) d ω
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 0 4 π 5 r ( l k δ i j + l j δ i k + l i δ k j ) d ω + V ( 3 l k l j l i ) d ω
As we discussed earlier in Section 3.1, the correct answer for the 2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) shall be
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( ( x k δ i j + x j δ i k + x i δ k j ) + 3 x k x j x i r 2 ) d ω
Therefore exchanging the Integration operator ∫ and Differentiation operator 2 x i x j under the assumption that the the incremental angular direction measure d ω and incremental volume measure d V is independent of the Differentiation operator variables x i , x j in the first integration and second integration respectively leads to wrong answer. Hence Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged under the assumption that the incremental angular direction measure d ω and incremental volume measure d V is independent of the Differentiation operator variables x i , x j in the first integration and second integration respectively.

3.2.6. Assuming that the Incremental Surface Area Measure d S = r ( l ) 2 d ω is Independent of the Differentiation Operator Variables x i , x j

In this subsection, we shall show that Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged when we assume that the incremental surface area measure d S = r ( l ) 2 d ω (subtended by a cone emanating from a point inside the volume V to the surface d S which is at a distance r ( l ) in the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) from the point) is independent of the Differentiation operator variables x i , x j . Note that we have earlier derived
S d x k ( d x l ( r ( l ) 3 ) ) n l d S = 3 S ( l l n l ) x k d S + 3 S ( r ( l ) δ l k ) n l d S
S d x k ( d x l ( r ( l ) 3 ) ) n l d S = 3 S ( l l x k ) n l d S + 3 S ( r ( l ) δ l k ) n l d S
Using Gauss Divergence theorem on the first integration, we can write the above equation as
S d x k ( d x l ( r ( l ) 3 ) ) n l d S = 3 V x l ( l l x k ) d V + 3 S ( r ( l ) δ l k ) n l d S
Now note that we have
x l ( l l x k ) = l l x l ( x k ) + x k x l ( l l ) = l l x l ( x k ) + x k x l ( x l r )
x l ( l l x k ) = l l x l ( x k ) + x k ( r x l ( x l ) x l x l r ) r 2 = l l x l ( x k ) + x k ( 3 r x l x l r ) r 2 = l l δ k l + x k ( 3 r r ) r 2
x l ( l l x k ) = 3 l k
Therefore we have
S d x k ( d x l ( r ( l ) 3 ) ) n l d S = 9 V l k d V + 3 S ( r ( l ) δ l k ) n l d S
S d x k ( d x l ( r ( l ) 3 ) ) n l d S = 9 V l k d V + 3 S r ( l ) n k d S
Now not that we have x k ( r ) = l k . Therefore using the Gauss Divergence theorem, we can write
V l k d V = V x k ( r ) d V = S r ( l ) n k d S
V l k d V = S r ( l ) n k d S
Therefore we can write
S d x k ( d x l ( r 3 ) ) n l d S = 12 S r ( l ) n k d S
Let’s see what happens when we interchange Integration operator ∫ and Differentiation operator 2 x i x j in the RHS of the above equation
2 x i x j ( S d x k ( d x l ( r ( l ) 3 ) ) n l d S ) = 12 2 x i x j ( S r ( l ) n k d S ) = 12 0 4 π 2 x i x j ( r ( l ) n k d S )
2 x i x j ( S d x k ( d x l ( r ( l ) 3 ) ) n l d S ) = 12 0 4 π 2 x i x j ( r ( l ) n k d S )
Since we assumed that the incremental surface area measure d S = r ( l ) 2 d ω (subtended by a cone emanating from a point inside the volume V to the surface d S which is at a distance r ( l ) in the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) from the point) is independent of the Differentiation operator variables x i , x j , we have x i x j ( d S ) = 0 . Therefore the above-mentioned boxed equation becomes
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π 2 x i x j ( r ( l ) n k ) d S
Now we know that
l k n k = r ( l ) 2 d ω = d S
Differentiating both sides of the equation with respect to x j we get
x j ( l k n k ) = x j ( d S )
Since we assumed that the incremental surface area measure d S = r ( l ) 2 d ω is independent of the Differentiation operator variables x i , x j , we have x j ( d S ) = 0
l k x j n k + n k x j l k = 0
l k x j n k + n k x j ( x k r ) = 0
l k x j n k + n k ( r x j x k x k x j r ) r 2 = 0
l k x j n k + n k ( r δ j k x k x j r ) r 2 = 0
x j n k = ( n j r n k x k x j r 3 ) l k
x j n k = ( n j n k x k x j r 2 ) x k
x j n k = n j x k n k x j r 2
Similarly, we also have l j n j = r ( l ) 2 d ω = d S and x i ( d S ) = 0 , therefore we can similarly write
l j x i n j + n j x i l j = 0
l j x i n j + n j x i ( x j r ) = 0
l j x i n j + n j ( r x i x j x j x i r ) r 2 = 0
l j x i n j + n j ( r δ i j x j x i r ) r 2 = 0
x i n j = ( n i r n j x j x i r 3 ) l k
x i n j = ( n i n j x j x i r 2 ) x j
x i n j = n i x j n j x i r 2
To compute the double derivative 2 x i x j ( r n k ) = x i ( x j ( r n k ) ) . We need to first compute the 2 derivatives x j ( r n k ) and x i ( r n j ) . Note that we can write using the product rule
x j ( r n k ) = r x j n k + n k x j r
x i ( r n j ) = r x i n j + n j x i r
Substituting the expression for x j n k that we got earlier in the above equation x j ( r n k ) = r x j n k + n k x j r we get
x j ( r n k ) = r ( n j x k n k x j r 2 ) + n k x j r = r n j x k n k x j r + n k x j r
x j ( r n k ) = r n j x k
Substituting the expression for x i n j that we got earlier in the equation x i ( r n j ) = r x i n j + n j x i r we get
x i ( r n j ) = r ( n i x j n j x i r 2 ) + n j x i r = r n i x j n j x i r + n j x i r
x i ( r n j ) = r n i x j
Using the expression for x j ( r n k ) that we derived earlier, we can write the double derivative 2 x i x j ( r n k ) = x i ( x j ( r n k ) ) as
2 x i x j ( r n k ) = x i ( x j ( r n k ) ) = x i ( r n j x k )
2 x i x j ( r n k ) = ( x k x i ( r n j ) r n j x i ( x k ) ) x k 2
Using the expression for x i ( r n j ) that we derived earlier,
2 x i x j ( r n k ) = ( r n i x k x j r n j δ k i ) x k 2
2 x i x j ( r n k ) = r n i x j x k r n j δ k i x k 2
2 x i x j ( r n k ) = r n i x j x k r n j x i 2
Now we have derived earlier that
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π 2 x i x j ( r ( l ) n k ) d S
Substituting the above-mentioned boxed equation into the previous equation we get
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( r n i x j x k r n j x i 2 ) d S
As we discussed earlier in Section 3.1, the correct answer for the 2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) shall be
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π ( ( x k δ i j + x j δ i k + x i δ k j ) + 3 x k x j x i r 2 ) d ω
Therefore exchanging the Integration operator ∫ and Differentiation operator 2 x i x j under the assumption that the incremental surface area measure d S = r ( l ) 2 d ω (subtended by a cone emanating from a point inside the volume V to the surface d S which is at a distance r ( l ) in the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) from the point) is independent of the Differentiation operator variables x i , x j leads to wrong answer. Hence Integration operator ∫ and Differentiation operator 2 x i x j cannot be exchanged under the assumption that the incremental surface area measure d S = r ( l ) 2 d ω is independent of the Differentiation operator variables x i , x j .

3.3. Method 3: Without Using Exchanging of Integration and Differentiation Operator

We now derive d x i d x j d x k d x l r 3 , where r = x 1 2 + x 2 2 + x 3 2 . Therefore we have
r 3 = ( x 1 2 + x 2 2 + x 3 2 ) 3 / 2
Step 1: The first derivative of r 3 with respect to x l is:
d x l r 3 = r 3 x l = 3 r 2 x l r = 3 r x l
Step 2: The second derivative is:
d x k 3 r x l = 3 r x k x l + r x l x k
= 3 x k r x l + r δ l k
= 3 x k x l r + r δ l k
Step 3: The third derivative is:
d x j 3 x k x l r + r δ l k = 3 x j x k x l r + x j r δ l k
x j x k x l r = δ j k x l + δ j l x k r x k x l x j r 3
x j r δ l k = x j r δ l k
d x j 3 x k x l r + r δ l k = 3 δ j k x l + δ j l x k r x k x l x j r 3 + x j r δ l k
Step 4: The fourth derivative is:
d x i 3 δ j k x l + δ j l x k r x k x l x j r 3 + x j r δ l k
For the first term:
d x i δ j k x l + δ j l x k r = δ j k δ i l + δ j l δ i k r δ j k x l + δ j l x k x i r 3
For the second term:
d x i x k x l x j r 3 = 3 x k x l x j x i r 5 + δ i k x l x j + δ i l x k x j + δ i j x k x l r 3
For the third term:
d x i x j δ l k r = δ i j δ l k r x j δ l k x i r 3
Combining all the terms, we get:
d x i d x j d x k d x l r 3 = 3 δ j k δ i l + δ j l δ i k + δ i j δ l k r δ j k x l x i + δ j l x k x i + δ i l x k x j + δ i k x l x j + δ l k x j x i + δ i j x k x l r 3 + 3 x k x l x j x i r 5
Therefore the integral S d x i d x j d x k d x l r 3 n l d S shall be
S d x i d x j d x k d x l r 3 n l d S = S 3 δ j k δ i l + δ j l δ i k + δ i j δ l k r δ j k x l x i + δ j l x k x i + δ i l x k x j + δ i k x l x j + δ l k x j x i + δ i j x k x l r 3 + 3 x k x l x j x i r 5 n l d S
S d x i d x j d x k d x l r 3 n l d S = S 3 δ j k δ i l + δ j l δ i k + δ i j δ l k r δ j k l l l i + δ j l l k l i + δ i l l k l j + δ i k l l l j + δ l k l j l i + δ i j l k l l r + 3 l k l l l j l i r n l d S
Note that we can write
n i l i d S = r 2 d ω
Therefore we can write
S d x i d x j d x k d x l r 3 n l d S = S 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k r ) n l d S + 0 4 π r ( 9 l k l j l i 3 ( δ j k l i + δ i k l j + δ i j l k ) ) d ω S 3 ( δ j l l k l i + δ i l l k l j + δ l k l j l i r ) n l d S
We shall compute each of the 2 surface integrals on the RHS one by one. Let us start by computing the 1st surface integral S 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k r ) n l d S . For this we define a tensor F in indicial notation as follows:
F i j k l = 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k r )
Using the Gauss Divergence theorem we can write
S ( F . n ) d S = V ( . F ) d V = V ( F i j k l x l ) d V
The volume integral V ( F i j k l x l ) d V can be written as
V ( F i j k l x l ) d V = V x l ( 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k r ) ) d V = 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k ) V x l ( 1 r ) d V = 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k ) V ( x l r 3 ) d V
Let us compute the volume integral by integrating over an elementary cone d ω centred on the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) with its vertex at x. The volume d V of this elementary cone is
d V = r 2 d r d ω
Therefore the integral V ( x l r 3 ) d V can be therefore written as
V ( x l r 3 ) d V = 0 4 π 0 r ( l ) x l r d r d ω = 0 4 π 0 r ( l ) l l d r d Ω = 0 4 π l l 0 r ( l ) d r d ω = 0 4 π l l r ( l ) d ω
V ( x l r 3 ) d V = 0 4 π r ( l ) l l d ω
Therefore the volume integral shall be
V ( F i j k l x l ) d V = 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k ) V ( x l r 3 ) d V = 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k ) 0 4 π r ( l ) l l d ω ( l )
V ( F i j k l x l ) d V = 3 0 4 π r ( l ) ( δ j k l i + δ i k l j + δ i j l k ) d ω ( l )
Therefore the surface integral shall be
S 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k r ) n l d S = 3 0 4 π r ( l ) ( δ j k l i + δ i k l j + δ i j l k ) d ω ( l )
Let us now compute the 2nd surface integral S ( δ j l l k l i + δ i l l k l j + δ l k l j l i r ) n l d S . For this define a tensor G in indicial notation as follows:
G i j k l = δ j l l k l i + δ i l l k l j + δ l k l j l i r
Using the Gauss Divergence theorem we can write
S ( G . n ) d S = V ( . G ) d V = V ( G i j k l x l ) d V
The volume integral V ( G i j k l x l ) d V can be written as
V ( G i j k l x l ) d V = V x l ( δ j l l k l i + δ i l l k l j + δ l k l j l i r ) d V
We need to do differentiate ( δ j l l k l i + δ i l l k l j + δ k l l j l i r ) before doing the volume integral. Given the function:
u = δ j l l k l i + δ i l l k l j + δ k l l j l i r
Substituting l i = x i r , l j = x j r , and l k = x k r into the above expression for u, we get:
u = 1 r ( δ j l x k r x i r + δ i l x k r x j r + δ k l x j r x i r )
u = 1 r 3 δ j l x k x i + δ i l x k x j + δ k l x j x i
Now, let’s compute the differentiation of u with respect to x l :
u x l = x l 1 r 3 δ j l x k x i + δ i l x k x j + δ k l x j x i + 1 r 3 x l δ j l x k x i + δ i l x k x j + δ k l x j x i
u x l = r 1 r 3 x l ( r ) δ j l x k x i + δ i l x k x j + δ k l x j x i + 1 r 3 x l δ j l x k x i + δ i l x k x j + δ k l x j x i
u x l = 3 x l r 5 δ j l x k x i + δ i l x k x j + δ k l x j x i + 1 r 3 x l δ j l x k x i + δ i l x k x j + δ k l x j x i
Now, let’s compute the differentiation x l δ j l x k x i + δ i l x k x j + δ k l x j x i term by term. For the first term we have:
x l δ j l x k x i = δ j l x l x k x i = δ j l x i x l x k + x k x l x i = δ j l δ k l x i + δ i l x k
x l δ j l x k x i = δ j k x i + δ j i x k
For the second term we have:
x l δ i l x k x j = δ i l x l x k x j = δ i l ( x l x k x j + x l x j x k ) = δ i l ( δ k l x j + δ j l x k )
x l δ i l x k x j = δ i k x j + δ i j x k
For the third term we have:
x l δ k l x j x i = δ k l x l x j x i = δ k l ( x i x l ( x j ) + x j x l ( x i ) ) = δ k l ( δ j l x i + δ i l x j )
x l δ k l x j x i = δ k j x i + δ k i x j
Therefore we have
x l δ j l x k x i + δ i l x k x j + δ k l x j x i = 2 ( δ k j x i + δ k i x j + δ i j x k )
Therefore the differentiation x l δ j l x k x i + δ i l x k x j + δ k l x j x i shall be
u x l = 3 x l r 5 δ j l x k x i + δ i l x k x j + δ k l x j x i + 2 r 3 ( δ k j x i + δ k i x j + δ i j x k )
u x l = 3 r 5 x j x k x i + x i x k x j + x k x j x i + 2 r 3 ( δ k j x i + δ k i x j + δ i j x k )
u x l = 9 r 5 x i x k x j + 2 r 3 ( δ k j x i + δ k i x j + δ i j x k )
u x l = 9 r 2 l i l k l j + 2 r 2 ( δ k j l i + δ k i l j + δ i j l k )
d d x l δ j l l k l i + δ i l l k l j + δ k l l j l i r = 1 r 2 ( 9 l i l k l j + 2 δ k j l i + 2 δ k i l j + 2 δ i j l k )
This is the fully simplified expression for the derivative. Therefore the volume integral shall be
V ( d d x l δ j l l k l i + δ i l l k l j + δ k l l j l i r ) d V = V 1 r 2 ( 9 l i l k l j + 2 δ k j l i + 2 δ k i l j + 2 δ i j l k ) d V
Let us compute the volume integral by integrating over an elementary cone d ω centred on the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) with its vertex at x. The volume d V of this elementary cone is
d V = r 2 d r d ω
Therefore the volume integral shall be
V ( d d x l δ j l l k l i + δ i l l k l j + δ k l l j l i r ) d V = 0 4 π 0 r ( l ) 1 r 2 ( 9 l i l k l j + 2 δ k j l i + 2 δ k i l j + 2 δ i j l k ) r 2 d r d ω
V ( d d x l δ j l l k l i + δ i l l k l j + δ k l l j l i r ) d V = 0 4 π ( 9 l i l k l j + 2 δ k j l i + 2 δ k i l j + 2 δ i j l k ) 0 r ( l ) d r d ω
V ( d d x l δ j l l k l i + δ i l l k l j + δ k l l j l i r ) d V = 0 4 π r ( l ) ( 9 l i l k l j + 2 δ k j l i + 2 δ k i l j + 2 δ i j l k ) d ω
Therefore the surface integral shall be
S ( δ j l l k l i + δ i l l k l j + δ l k l j l i r ) n l d S = 0 4 π r ( l ) ( 9 l i l k l j + 2 δ k j l i + 2 δ k i l j + 2 δ i j l k ) d ω
Now we earlier derived that
S d x i d x j d x k d x l r 3 n l d S = S 3 ( δ j k δ i l + δ j l δ i k + δ i j δ l k r ) n l d S + 0 4 π r ( 9 l k l j l i 3 ( δ j k l i + δ i k l j + δ i j l k ) ) d ω S 3 ( δ j l l k l i + δ i l l k l j + δ l k l j l i r ) n l d S
Substituting the expressions of the 2 surface integrals that we derived earlier into the above equation we get
S d x i d x j d x k d x l r 3 n l d S = 3 0 4 π r ( l ) ( δ j k l i + δ i k l j + δ i j l k ) d ω ( l ) + 0 4 π r ( l ) ( 9 l k l j l i 3 ( δ j k l i + δ i k l j + δ i j l k ) ) d ω ( l )
3 0 4 π r ( l ) ( 9 l i l k l j + 2 δ k j l i + 2 δ k i l j + 2 δ i j l k ) d ω ( l )
S d x i d x j d x k d x l r 3 n l d S = 0 4 π r ( l ) ( 36 l k l j l i 12 δ j k l i 12 δ i k l j 12 δ i j l k ) d ω ( l )
S d x i d x j d x k d x l r 3 n l d S = 12 0 4 π r ( l ) ( 3 l k l j l i δ j k l i δ i k l j δ i j l k ) d ω ( l )

4. Betti’s Theorem and Reciprocity

Betti’s Theorem is a fundamental result in linear elasticity that provides a relationship between two different states of stress and strain within an elastic body. The theorem states that if a body is subjected to two different sets of equilibrating forces, the work done by one set of forces during the displacements caused by the other set is the same as the work done by the second set during the displacements caused by the first set.
Let:
  • u ( 1 ) be the displacement field due to traction force t ( 1 ) and body force b ( 1 ) .
  • u ( 2 ) be the displacement field due to traction force t ( 2 ) and body force b ( 2 ) .
Then Betti’s Theorem is expressed as:
S t i ( 1 ) u i ( 2 ) d S + V b i ( 1 ) u i ( 2 ) d V = S t i ( 2 ) u i ( 1 ) d S + V b i ( 2 ) u i ( 1 ) d V
Or, in component form:
S t i ( 1 ) u i ( 2 ) d S + V b i ( 1 ) u i ( 2 ) d V = S t i ( 2 ) u i ( 1 ) d S + V b i ( 2 ) u i ( 1 ) d V ( 1.82 )
To prove Betti’s Theorem, we begin by noting the stress-strain relationships in each state:
σ i j ( 1 ) = C i j k l e k l ( 1 ) and σ i j ( 2 ) = C i j k l e k l ( 2 )
The strain energy densities for the two states are given by:
σ i j ( 1 ) e i j ( 2 ) = C i j k l e k l ( 1 ) e i j ( 2 ) and σ i j ( 2 ) e i j ( 1 ) = C i j k l e k l ( 2 ) e i j ( 1 )
Since C i j k l is symmetric under interchange of the first and second pairs of indices (i.e., C i j k l = C k l i j ), we have:
σ i j ( 1 ) e i j ( 2 ) = σ i j ( 2 ) e i j ( 1 )
Integrating this over the volume V, we obtain:
V σ i j ( 1 ) e i j ( 2 ) d V = V σ i j ( 2 ) e i j ( 1 ) d V ( 1.83 )
Expanding the left-hand side:
V σ i j ( 1 ) e i j ( 2 ) d V = V σ i j ( 1 ) u j , i ( 2 ) d V
Applying the divergence theorem:
V σ i j ( 1 ) u j , i ( 2 ) d V = V ( σ i j ( 1 ) u j ( 2 ) ) , i d V V σ i j , i ( 1 ) u j ( 2 ) d V
Using the equilibrium condition σ i j , i ( 1 ) + b i ( 1 ) = 0 :
V σ i j ( 1 ) e i j ( 2 ) d V = S t i ( 1 ) u i ( 2 ) d S + V b i ( 1 ) u i ( 2 ) d V ( 1.85 )
Similarly, for the second state:
V σ i j ( 2 ) e i j ( 1 ) d V = S t i ( 2 ) u i ( 1 ) d S + V b i ( 2 ) u i ( 1 ) d V
Equating these expressions gives Betti’s Theorem:
S t i ( 1 ) u i ( 2 ) d S + V b i ( 1 ) u i ( 2 ) d V = S t i ( 2 ) u i ( 1 ) d S + V b i ( 2 ) u i ( 1 ) d V
Betti’s Theorem can be used to derive the reciprocity relation for the Green’s function in elasticity. The Green’s function G i j ( x , x ) represents the displacement in the i-th direction at point x due to a unit point force in the j-th direction at point x . The reciprocity relation states that:
G i j ( x , x ) = G j i ( x , x ) ( 1.86 )
Consider two specific states:
  • b i ( 1 ) = F i δ ( x x ( 1 ) ) , i.e., a point force F i at x ( 1 ) .
  • b i ( 2 ) = H i δ ( x x ( 2 ) ) , i.e., a point force H i at x ( 2 ) .
The corresponding displacement fields are:
u i ( 1 ) ( x ) = G i j ( x , x ( 1 ) ) F j and u i ( 2 ) ( x ) = G i j ( x , x ( 2 ) ) H j
Substituting these into Betti’s Theorem:
V F i δ ( x x ( 1 ) ) G i j ( x , x ( 2 ) ) H j d V = V H j δ ( x x ( 2 ) ) G j i ( x , x ( 1 ) ) F i d V
Simplifying using the properties of the delta function:
F i H j G i j ( x ( 1 ) , x ( 2 ) ) = F i H j G j i ( x ( 2 ) , x ( 1 ) )
Since F i and H j are arbitrary, this implies the reciprocity of the Green’s function:
G i j ( x ( 1 ) , x ( 2 ) ) = G j i ( x ( 2 ) , x ( 1 ) ) ( 1.90 )
The derivation provided above rigorously covers the theory and equations presented in Section 1.7 of the document. Betti’s Theorem provides a powerful tool for understanding the relationship between different states of stress and strain in an elastic body, and it leads directly to the important reciprocity relation for Green’s functions, which is essential in solving a wide range of problems in elasticity.

5. Derivation of the Green’s Function for Isotropic Medium

To derive the Love 1927 Solution of Displacement we need to first derive the Navier-Cauchy equations. The Navier-Cauchy equations describe the equilibrium state of an isotropic elastic medium under applied forces. The equilibrium equation in the absence of body forces is:
· σ + f = 0
where σ is the stress tensor, and f represents the body forces. For an isotropic material, the stress tensor σ is related to the strain tensor ϵ via Hooke’s law:
σ i j = λ δ i j ϵ k k + 2 μ ϵ i j
where λ and μ are the Lamé constants, and the strain tensor ϵ i j is:
ϵ i j = 1 2 u i x j + u j x i
Substituting the strain tensor into the stress tensor and then into the equilibrium equation gives:
x j λ δ i j u k x k + 2 μ ϵ i j + f i = 0
Expanding this, we have:
λ x i u k x k + μ 2 u i x j 2 + μ 2 u j x i x j + f i = 0
Recognizing that 2 u j x i x j = 2 u i x j 2 , we simplify to:
( λ + μ ) x i u k x k + μ 2 u i + f i = 0
In vector form, this is:
μ 2 u + ( λ + μ ) ( · u ) + f = 0
This is the Navier-Cauchy equation in its general form. To find the Green’s function G i j ( x , x 0 ) , we consider the Navier-Cauchy equation with a point force applied at x 0 :
μ 2 G i j ( x , x 0 ) + ( λ + μ ) x i G k j ( x , x 0 ) x k = δ i j δ ( x x 0 )
Here, G i j ( x , x 0 ) represents the ith component of displacement at x due to a unit point force in the jth direction at x 0 . We take the Fourier transform F of the Navier-Cauchy equation. The Fourier transform of a function f ( r ) in three dimensions is defined by:
F { f ( r ) } ( k ) = f ^ ( k ) = R 3 e i k · r f ( r ) d r
We will now compute the Fourier transform of each term in the Navier-Cauchy equation separately. The term 2 u is the Laplacian of the displacement field. In three dimensions, the Laplacian of a vector field u ( r ) is:
2 u = 2 u i x j 2
Taking the Fourier transform:
F { 2 u i ( r ) } ( k ) = F 2 u i ( r ) x j 2 ( k )
Using the property of Fourier transforms:
F 2 u i ( r ) x j 2 ( k ) = k j 2 u ^ i ( k )
we have:
F { 2 u i ( r ) } ( k ) = k 2 u ^ i ( k )
where k 2 = k j k j . The term ( · u ) involves the gradient of the divergence of the displacement field. First, compute the Fourier transform of · u :
F { · u ( r ) } ( k ) = F u i x i ( k ) = i k i u ^ i ( k )
Now, take the gradient:
F { ( · u ) } ( k ) = F x j u i x i ( k ) = k j k i u ^ i ( k )
For completeness, we also compute the Fourier transform of the body force term f :
F { f ( r ) } ( k ) = f ^ ( k )
In the case of a point force applied at r = r 0 , this becomes:
F { δ ( r r 0 ) } ( k ) = e i k · r 0
Substituting these results into the Navier-Cauchy equation:
μ ( k 2 u ^ i ( k ) ) + ( λ + μ ) ( k i k j u ^ j ( k ) ) = f ^ i ( k )
μ k 2 u ^ i ( k ) ( λ + μ ) k i k j u ^ j ( k ) = f ^ i ( k )
μ k 2 u ^ i ( k ) + ( λ + μ ) k i k j u ^ j ( k ) = f ^ i ( k )
In matrix form, this can be written as:
μ k 2 δ i j + ( λ + μ ) k i k j u ^ j ( k ) = f ^ i ( k )
To find u ^ j ( k ) , we need to invert the matrix A i j = μ k 2 δ i j + ( λ + μ ) k i k j . The matrix A i j can be decomposed into:
A i j = μ k 2 δ i j + ( λ + μ ) k i k j
We need to find A i j 1 such that A i j A j k 1 = δ i k . We propose the inverse to be of the form:
A i j 1 = 1 μ k 2 α δ i j + β k i k j k 2
Multiplying A i j with the proposed inverse A i j 1 :
A i j A i j 1 = μ k 2 δ i j + ( λ + μ ) k i k j α δ j k + β k j k k k 2 μ k 2
A i j A i j 1 = α δ i k + β + α ( λ + μ ) μ + β ( λ + μ ) μ k i k k k 2
For this to equal δ i k , we need:
α = 1 , β = λ + μ λ + 2 μ
Therefore we have:
A i j 1 = 1 μ k 2 δ i j λ + μ λ + 2 μ k i k j k 2
The displacement field u ^ j ( k ) is then:
u ^ j ( k ) = A i j 1 f ^ i ( k )
Substituting the inverse:
u ^ j ( k ) = e i k · r 0 μ k 2 δ i j λ + μ λ + 2 μ k i k j k 2
The Green’s function G i j ( r , r 0 ) in real space is given by the inverse Fourier transform:
G i j ( r , r 0 ) = F 1 G ^ i j ( k ) ( r ) = 1 ( 2 π ) 3 R 3 G ^ i j ( k ) e i k · ( r r 0 ) d k
Substituting G ^ i j ( k ) :
G i j ( r , r 0 ) = 1 ( 2 π ) 3 R 3 1 μ k 2 δ i j λ + μ λ + 2 μ k i k j k 2 e i k · ( r r 0 ) d k
The first integral involves R 3 1 k 2 e i k · r d k . This is known to yield:
F 1 1 | k | 2 ( r ) = 1 4 π r
Let’s rigorously derive the Fourier transform of 1 | k | 2 in three dimensions. The inverse Fourier transform of f ^ ( k ) = 1 | k | 2 is:
f ( r ) = F 1 1 | k | 2 ( r ) = 1 ( 2 π ) 3 R 3 e i k · r | k | 2 d k
Here, k and r are three-dimensional vectors. The function 1 | k | 2 is spherically symmetric in k -space. Therefore, we can simplify the problem by transforming to spherical coordinates. In three dimensions, the spherical coordinates are:
k = ( k x , k y , k z ) = ( k sin θ cos ϕ , k sin θ sin ϕ , k cos θ )
where: k = | k | , θ is the polar angle (angle with the z-axis), ϕ is the azimuthal angle (angle in the x y -plane). The volume element in spherical coordinates is d k = k 2 sin θ d k d θ d ϕ . The integral becomes:
f ( r ) = 1 ( 2 π ) 3 0 0 π 0 2 π e i k r cos θ k 2 k 2 sin θ d ϕ d θ d k
Here, r = | r | and without loss of generality, we have aligned r along the z-axis, so k · r = k r cos θ . Simplifying the integrand:
f ( r ) = 1 ( 2 π ) 3 0 0 π 0 2 π e i k r cos θ sin θ d ϕ d θ d k
The integral over ϕ is straightforward 0 2 π d ϕ = 2 π . This simplifies the expression to:
f ( r ) = 1 ( 2 π ) 2 0 0 π e i k r cos θ sin θ d θ d k
Next, consider the integral over the polar angle θ :
0 π e i k r cos θ sin θ d θ
Let u = cos θ , hence d u = sin θ d θ , and the limits change from θ = 0 to θ = π correspond to u = 1 to u = 1 :
0 π e i k r cos θ sin θ d θ = 1 1 e i k r u ( d u ) = 1 1 e i k r u d u
This is a standard integral 1 1 e i k r u d u = 2 sin ( k r ) k r . Thus, the expression for f ( r ) becomes:
f ( r ) = 1 ( 2 π ) 2 r 0 sin ( k r ) d k
Finally, we evaluate the integral over k:
0 sin ( k r ) k d k = π 2
Thus:
f ( r ) = 1 ( 2 π ) 2 r · π 2 = 1 4 π r
The inverse Fourier transform of 1 k 2 is:
F 1 1 | k | 2 ( r ) = 1 4 π r
So, the first integral contributes:
1 μ δ i j 4 π r
The second integral involves R 3 k i k j k 4 e i k · r d k . The Fourier transform of k i k j / k 4 gives:
R 3 k i k j k 4 e i k · r d k = 2 r i r j 4 π | r | = 3 r i r j r 2 δ i j r 5
Combining the results, the Green’s function for an isotropic elastic medium is:
G i j ( r , r 0 ) = 1 8 π μ r ( 3 4 ν ) δ i j + r i r j r 2
where r = | r r 0 | and ν is Poisson’s ratio. This derivation rigorously follows from the Navier-Cauchy equations, through the application of Fourier transforms, detailed inversion of the resulting matrix in Fourier space, and evaluation of inverse Fourier transforms. The final Green’s function describes the displacement field due to a point force in an infinite isotropic elastic medium.

6. Displacement due to Spontaneous Change of Form of Inclusion

The Displacement at r due to point force F i at r is (Love 1927)
U i ( r r ) = 1 4 π μ F i | r r | 1 16 π μ ( 1 σ ) F j 2 x j x i | r r |
which is Equation (2.5) of Eshelby’s classic paper 1957. Note that we have earlier derived
d x j ( d x i ( | r r | ) ) = δ i j | r r | x i x j | r r | 3
The Displacement (Love 1927) at r due to point force F i at r can be therefore written as
U i ( r r ) = 1 4 π μ F i | r r | 1 16 π μ ( 1 σ ) F j ( δ i j | r r | x i x j | r r | 3 )
U i ( r r ) = 1 4 π μ F j δ j i | r r | 1 16 π μ ( 1 σ ) F j ( δ j i | r r | x i x j | r r | 3 )
U i ( r r ) = F j 16 π μ ( 1 σ ) 4 ( 1 σ ) δ j i | r r | F j 16 π μ ( 1 σ ) δ j i | r r | + 1 16 π μ ( 1 σ ) F j x i x j | r r | 3
U i ( r r ) = F j 16 π μ ( 1 σ ) ( 3 4 σ ) δ j i | r r | + 1 16 π μ ( 1 σ ) F j x i x j | r r | 3
U i ( r r ) = F j 16 π μ ( 1 σ ) ( ( 3 4 σ ) | r r | δ j i + x i x j | r r | 3 )
which is Equation (2.14) of Eshelby’s classic paper 1957. There are multiple ways to do the surface integral over S of the above-boxed quantity. We shall mention 2 different methods in this article

6.1. Method 1: Simple Method using only Gauss Divergence Theorem

In Stage III, we apply a force distribution F j = p j k T n k over S to make the body free of external force (but in a state of self-stress because of the transformation of the inclusion). The displacement impressed on the material in stage III is
u i C = S U i ( r r ) d S
u i C = S ( F j 16 π μ ( 1 σ ) ( ( 3 4 σ ) | r r | δ j i + x i x j | r r | 3 ) ) d S
u i C = 1 16 π μ ( 1 σ ) S F j ( ( 3 4 σ ) | r r | δ j i + x i x j | r r | 3 ) d S
u i C = 1 16 π μ ( 1 σ ) S p j k T ( ( 3 4 σ ) | r r | δ j i + x i x j | r r | 3 ) n k d S
Let’s proceed with a rigorous application of the Gauss Divergence Theorem to compute the given surface integral. We need to evaluate the surface integral:
u i C = 1 16 π μ ( 1 σ ) S p j k T 3 4 σ r δ i j + x j x i r 3 n k d S
where:
  • p j k T is the stress tensor, which is constant throughout the volume.
  • n k is the k-th component of the unit normal vector to the surface S.
  • δ i j is the Kronecker delta.
  • r = x 1 2 + x 2 2 + x 3 2 is the radial distance from the origin.
Let us now apply the Gauss Divergence Theorem. The Gauss Divergence Theorem states:
S F k n k d S = V F k x k d V
where F k is a vector field. However, in this problem, we have a tensor field F i k defined as:
F i k = p j k T 3 4 σ r δ i j + x j x i r 3
The surface integral becomes:
I = S F i k n k d S
Using the Gauss Divergence Theorem for tensor fields, this surface integral can be converted into a volume integral:
I = V F i k x k d V
We need to compute the divergence of the Tensor Field F i k x k :
F i k x k = x k p j k 3 4 σ r δ i j + x j x i r 3
Since p j k T is constant, it can be factored out:
F i k x k = p j k T x k 3 4 σ r δ i j + x j x i r 3
Let’s differentiate each term of the above-mentioned boxed equation. The differentiation of the first term shall be
x k 3 4 σ r δ i j = δ i j x k 3 4 σ r
Now note that we have
x k 1 r = x k r 3
Therefore we can write
δ i j x k 3 4 σ r = δ i j x k ( 3 4 σ ) r 3
The differentiation of the second term shall be
x k x j x i r 3
This expands using the product rule:
x k x j x i r 3 = x j x k · x i r 3 + x j x k x i r 3
Since x j x k = δ j k , the above equation can be written as
x k x j x i r 3 = δ j k x i r 3 + x j x k x i r 3
Note that for the second term of the above-equation we can write
x k x i r 3 = δ i k r 3 3 x i x k r 5
Thus, the divergence of the tensor field F i k becomes:
F i k x k = p j k T δ i j x k ( 3 4 σ ) r 3 + δ j k x i r 3 + x j δ i k r 3 3 x i x k r 5
The surface integral now transforms into a volume integral:
u i C = 1 16 π μ ( 1 σ ) V p j k T δ i j x k ( 3 4 σ ) r 3 + δ j k x i r 3 + x j δ i k r 3 3 x i x k r 5 d V
Now note that the above volume integral can be written as
u i C = 1 16 π μ ( 1 σ ) V p j k T δ i j x k ( 4 4 σ ) r 3 + δ i j x k r 3 + δ j k x i r 3 + x j δ i k r 3 3 x i x k r 5 d V
u i C = 1 16 π μ ( 1 σ ) V p j k T δ i j x k ( 4 4 σ ) r 3 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j x k r 3 + δ j k x i r 3 + δ i k x j r 3 3 x i x j x k r 5 ) d V
u i C = 1 16 π μ ( 1 σ ) V p i k T x k ( 4 4 σ ) r 3 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j x k r 3 + δ j k x i r 3 + δ i k x j r 3 3 x i x j x k r 5 ) d V
u i C = 1 16 π μ ( 1 σ ) V p i k T l k ( 4 4 σ ) r 2 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j l k r 2 + δ j k l i r 2 + δ i k l j r 2 3 l i l j l k r 2 ) d V
u i C = 1 16 π μ ( 1 σ ) V p i k T l k ( 4 4 σ ) r 2 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
u i C = 1 16 π μ ( 1 σ ) V ( p i k l k + p i j T l j ) ( 2 2 σ ) r 2 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
u i C = 1 16 π μ ( 1 σ ) V ( p i k T l k + p i j T l j ) 2 ( 1 σ ) r 2 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
Due to balance of the angular momentum, the transformation stress matrix p T is symmetric, i.e
p i j T = p j i T
Therefore we can write
u i C = 1 16 π μ ( 1 σ ) V ( p i k T l k + p j i T l j ) 2 ( 1 σ ) r 2 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
Now we can write
p i k T l k = p j k T δ i j l k
p j i T l j = p j k T δ i k l j
Therefore we can write
u i C = 1 16 π μ ( 1 σ ) V ( p j k T δ i j l k + p j k T δ i k l j ) 2 ( 1 σ ) r 2 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
u i C = 1 16 π μ ( 1 σ ) V p j k T ( δ i j l k + δ i k l j ) 2 ( 1 σ ) r 2 d V + 1 16 π μ ( 1 σ ) V p j k T ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
Since p j k T is constant inside the volume V, we can take p j k T outside the volume integral
u i C = p j k T 16 π μ ( 1 σ ) V ( δ i j l k + δ i k l j ) 2 ( 1 σ ) r 2 d V + p j k T 16 π μ ( 1 σ ) V ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
u i C = p j k T 16 π μ ( 1 σ ) V 2 ( δ i j l k + δ i k l j ) ( 1 σ ) + ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
u i C = p j k T 16 π μ ( 1 σ ) V ( δ i j l k + δ i k l j ) ( 1 2 σ ) ( δ i j l k + δ i k l j ) + ( δ i j l k + δ j k l i + δ i k l j 3 l i l j l k ) r 2 d V
u i C = p j k T 16 π μ ( 1 σ ) V ( δ i j l k + δ i k l j ) ( 1 2 σ ) δ i j l k δ i k l j + δ i j l k + δ j k l i + δ i k l j 3 l i l j l k r 2 d V
u i C = p j k T 16 π μ ( 1 σ ) V ( δ i j l k + δ i k l j ) ( 1 2 σ ) + δ j k l i 3 l i l j l k r 2 d V
u i C = p j k T 16 π μ ( 1 σ ) V ( 1 2 σ ) ( δ i j l k + δ i k l j ) δ j k l i + 3 l i l j l k r 2 d V
which is the first part of the Equation (2.15) of Eshelby’s Classic Paper.

6.2. Method 2: Eshelby’s Method of using the Gauss Divergence Theorem and a special variation of Stokes’s theorem

Note the following 2 important identities in spherical co-ordiantes ( r , θ , ϕ )
2 r = 1 r 2 r ( r 2 r ( r ) ) + 1 r 2 sin θ θ ( sin θ θ ( r ) ) + 1 r 2 sin 2 θ 2 θ 2 ( r ) = 2 r
2 r 3 = 1 r 2 r ( r 2 r ( r 3 ) ) + 1 r 2 sin θ θ ( sin θ θ ( r 3 ) ) + 1 r 2 sin 2 θ 2 θ 2 ( r ) = 12 r
The above-mentioned 2 identities can be written in cartesian co-ordinates as
2 x l x l ( r ) = 2 r
2 x l x l ( r 3 ) = 12 r
The Displacement (Love 1927) at r due to point force F i at r can be therefore written as
U i ( r r ) = 1 8 π μ F i 2 x l x l ( r ) 1 192 π μ ( 1 σ ) F j 2 x j x i ( 2 x l x l ( r 3 ) )
U i ( r r ) = 1 8 π μ F i 2 x l x l ( r ) 1 192 π μ ( 1 σ ) F j 4 x j x i x l x l ( r 3 ) )
U i ( r r ) = 1 8 π μ F j δ j i 2 x l x l ( r ) 1 192 π μ ( 1 σ ) F j 4 x j x i x l x l ( r 3 )
U i ( r r ) = 1 192 π μ ( 1 σ ) F j ( 24 ( 1 σ ) 2 x l x l ( r ) δ j i 4 x j x i x l x l ( r 3 ) )
Now note that we have F j = p j k T n k . Therefore the above-boxed identity can be written as
U i ( r r ) = 1 192 π μ ( 1 σ ) p j k T ( 24 ( 1 σ ) 2 x l x l ( r ) δ j i n k 4 x j x i x l x l ( r 3 ) n k )
In Stage III, we apply a force distribution F j = p j k T n k over S to make the body free of external force (but in a state of self-stress because of the transformation of the inclusion). The displacement impressed on the material in stage III is
u i C = S U i ( r r ) d S
Substituting the equation of displacement U i ( r r ) that we got earlier in the above equation we get
u i C = S ( 1 192 π μ ( 1 σ ) p j k T ( 24 ( 1 σ ) 2 x l x l ( r ) δ j i n k 4 x j x i x l x l ( r 3 ) n k ) ) d S
u i C = p j k T δ j i 8 π μ S ( 2 x l x l ( r ) n k ) d S p j k T 192 π μ ( 1 σ ) S ( 4 x j x i x l x l ( r 3 ) n k ) d S
We need to simplify the 2 integrals S ( 2 x l x l ( r ) n k ) d S and S ( 4 x j x i x l x l ( r 3 ) n k ) d S . For this, we use the Gauss Divergence Theorem twice once to convert the surface integral to volume integral and then convert the volume integral back again to surface integral but with a different surface vector. The integral S ( 2 x l x l ( r ) n k ) d S can be written as
S ( 2 x l x l ( r ) n k ) d S = V x k ( 2 x l x l ( r ) ) d V = V x l ( 2 x k x l ( r ) ) d V = S ( 2 x k x l ( r ) n l ) d S
S ( 2 x l x l ( r ) n k ) d S = S ( 2 x k x l ( r ) n l ) d S
The integral S ( 4 x j x i x l x l ( r 3 ) n k ) d S can be written as
S ( 4 x j x i x l x l ( r 3 ) n k ) d S = V x k ( 4 x j x i x l x l ( r 3 ) ) d V = V x l ( 4 x j x i x k x l ( r 3 ) ) d V = S ( 4 x j x i x k x l ( r 3 ) n l ) d S
S ( 4 x j x i x l x l ( r 3 ) n k ) d S = S ( 4 x j x i x k x l ( r 3 ) n l ) d S
From the above 2 boxed equations, we can write the displacement u i C as
u i C = p j k T δ j i 8 π μ S ( 2 x l x k ( r ) n l ) d S p j k T 192 π μ ( 1 σ ) S ( 4 x j x i x l x k ( r 3 ) n l ) d S
Note that we earlier derived the following 2 identities
S d x l ( d x k ( r ) ) n l d S = 2 0 4 π r ( l ) l k d ω ( l )
2 x i x j ( S d x k ( d x l ( r 3 ) ) n l d S ) = 12 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l )
Substituting the above 2 identities in the equation for the displacement u i C , we get
u i C = p j k T δ j i 8 π μ ( 2 0 4 π r ( l ) l k d ω ( l ) ) p j k T 192 π μ ( 1 σ ) ( 12 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l ) )
u i C = p j k T δ j i 4 π μ ( 0 4 π r ( l ) l k d ω ( l ) ) p j k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l ) )
u i C = p i k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) 4 ( 1 σ ) l k d ω ( l ) ) p j k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l ) )
Since p i k T is constant throughout the volume of the inclusion, we can take p i k T inside of the integral. Therefore we have
u i C = 1 16 π μ ( 1 σ ) ( 0 4 π r ( l ) 4 ( 1 σ ) p i k T l k d ω ( l ) ) p j k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l ) )
u i C = 1 16 π μ ( 1 σ ) ( 0 4 π r ( l ) 2 ( 1 σ ) ( p i k T l k + p i j T l j ) d ω ( l ) ) p j k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l ) )
Due to balance of the angular momentum, the transformation stress matrix p T is symmetric, i.e
p i j T = p j i T
Therefore we can write
u i C = 1 16 π μ ( 1 σ ) ( 0 4 π r ( l ) 2 ( 1 σ ) ( p i k T l k + p j i T l j ) d ω ( l ) ) p j k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l ) )
Now we can write
p i k T l k = p j k T δ i j l k
p j i T l j = p j k T δ i k l j
Therefore we can write
u i C = 1 16 π μ ( 1 σ ) ( 0 4 π r ( l ) 2 ( 1 σ ) ( p j k T δ i j l k + p j k T δ i k l j ) d ω ( l ) ) p j k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l ) )
u i C = p j k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) 2 ( 1 σ ) ( δ i j l k + δ i k l j ) d ω ( l ) ) p j k T 16 π μ ( 1 σ ) ( 0 4 π r ( l ) ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l ) )
u i C = p j k T 16 π μ ( 1 σ ) 0 4 π r ( l ) ( 2 ( 1 σ ) ( δ i j l k + δ i k l j ) + ( ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) ) d ω ( l )
u i C = p j k T 16 π μ ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( δ i j l k + δ i k l j ) + ( δ i j l k + δ i k l j ) ( l i δ j k + l k δ j i + l j δ k i ) + 3 l k l j l i ) d ω ( l )
u i C = p j k T 16 π μ ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( δ i j l k + δ i k l j ) δ j k l i + 3 l k l j l i ) d ω ( l )
To convert this into volume integral we note that the incremental solid angle is related to the incremental volume by the following relation
r ( l ) d ω ( l ) = d V r ( l ) 2
Therefore the equation for displacement u i C can be written as
u i C = p j k T 16 π μ ( 1 σ ) V ( ( 1 2 σ ) ( δ i j l k + δ i k l j ) δ j k l i + 3 l k l j l i ) r ( l ) 2 d V
which is the first part of the Equation (2.15) of Eshelby’s Classic Paper. To derive Equation (3.1) and also the second part of the Equation (2.15), we need to first state Hooke’s law
p j k T = λ e m m T δ j k + 2 μ e j k T
Substituting the above equation into the equation of displacement u i C that we derived earlier we get
u i C = ( λ e m m T δ j k + 2 μ e j k T ) 16 π μ ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( δ i j l k + δ i k l j ) δ j k l i + 3 l k l j l i ) d ω ( l )
u i C = 1 16 π μ ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( λ e m m T δ j k + 2 μ e j k T ) ( δ i j l k + δ i k l j ) ( λ e m m T δ j k + 2 μ e j k T ) δ j k l i + 3 ( λ e m m T δ j k + 2 μ e j k T ) l k l j l i ) d ω ( l )
We can simplify the first term of the integrand in the above equation as
r ( l ) ( 1 2 σ ) ( λ e m m T δ j k + 2 μ e j k T ) ( δ i j l k + δ i k l j ) = r ( l ) ( 1 2 σ ) ( λ e m m T δ j k δ i j l k + 2 μ e j k T δ i j l k + λ e m m δ j k δ i k l j + 2 μ e j k T δ i k l j )
r ( l ) ( 1 2 σ ) ( λ e m m T δ j k + 2 μ e j k T ) ( δ i j l k + δ i k l j ) = r ( l ) ( 1 2 σ ) ( λ e m m T δ i k l k + 2 μ e i k T l k + λ e m m T δ j i l j + 2 μ ϵ j i l j )
r ( l ) ( 1 2 σ ) ( λ e m m T δ j k + 2 μ e j k T ) ( δ i j l k + δ i k l j ) = r ( l ) ( 1 2 σ ) ( λ e m m T l i + 2 μ e i k T l k + λ e m m T l i + 2 μ e j i T l j )
r ( l ) ( 1 2 σ ) ( λ e m m T δ j k + 2 μ e j k T ) ( δ i j l k + δ i k l j ) = r ( l ) ( 1 2 σ ) ( 2 λ e m m T l i + 2 μ ( e i k T l k + e j i T l j ) )
Similarly the second term of the integrand in the above equation can be simplified as
r ( l ) ( λ e m m T δ j k + 2 μ e j k T ) δ j k l i = r ( l ) ( λ e m m T δ j k δ j k l i + 2 μ e j k T δ j k l i )
r ( l ) ( λ e m m T δ j k + 2 μ e j k T ) δ j k l i = r ( l ) ( λ e m m T δ j j l i + 2 μ e k k T l i )
r ( l ) ( λ e m m T δ j k + 2 μ e j k T ) δ j k l i = r ( l ) ( 3 λ e m m T l i + 2 μ e k k T l i )
Similarly the third term of the integrand in the above equation can be simplified as
r ( l ) ( 3 ( λ e m m T δ j k + 2 μ e j k T ) l k l j l i ) = r ( l ) ( 3 λ e m m T δ j k l k l j l i + 6 μ e j k T l k l j l i )
r ( l ) ( 3 ( λ e m m δ j k + 2 μ ϵ j k ) l k l j l i ) = r ( l ) ( 3 λ e m m l j l j l i + 6 μ ϵ j k l k l j l i )
Since l j l j = x j x j r = r 2 r = r , we can write the above equation as
r ( l ) ( 3 ( λ e m m δ j k + 2 μ ϵ j k ) l k l j l i ) = r ( l ) ( 3 λ e m m l i + 6 μ ϵ j k l k l j l i )
Therefore we can write the equation of displacement u i C that we derived earlier as
u i C = 1 16 π μ ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( 2 λ e m m T l i + 2 μ ( e i k T l k + e j i T l j ) ) ( 3 λ e m m T l i + 2 μ e k k T l i ) + ( 3 λ e m m T l i + 6 μ e j k T l k l j l i ) ) d ω ( l )
u i C = 1 16 π μ ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( 2 λ e m m T l i + 2 μ ( e i k T l k + e j i T l j ) ) 3 λ e m m T l i 2 μ e k k T l i + 3 λ e m m T l i + 6 μ e j k T l k l j l i ) d ω ( l )
u i C = 1 16 π μ ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( 2 λ e m m T l i ) + 2 μ ( 1 2 σ ) ( e i k T l k + e j i T l j ) 2 μ e k k T l i + 6 μ e j k T l k l j l i ) d ω ( l )
Now note that we have e m m T l i = e k k T l i . Therefore we can write the above equation as
u i C = 1 16 π μ ( 1 σ ) 0 4 π r ( l ) ( 2 e m m T l i ( λ ( 1 2 σ ) μ ) + 2 μ ( 1 2 σ ) ( e i k T l k + e j i T l j ) + 6 μ e j k T l k l j l i ) d ω ( l )
Now note that we can express the first lame parameter λ in terms of Poisson’s ratio σ and the second lame parameter μ as
λ = 2 μ σ ( 1 2 σ )
Therefore we can write
λ ( 1 2 σ ) = 2 μ σ
λ ( 1 2 σ ) μ = 2 μ σ μ = μ ( 2 σ 1 )
λ ( 1 2 σ ) μ = μ ( 2 σ 1 )
Therefore we can write the equation for displacement u i C as
u i C = 1 16 π μ ( 1 σ ) 0 4 π r ( l ) ( 2 e m m T l i μ ( 2 σ 1 ) + 2 μ ( 1 2 σ ) ( e i k T l k + e j i T l j ) + 6 μ e j k T l k l j l i ) d ω ( l )
u i C = 1 8 π ( 1 σ ) 0 4 π r ( l ) ( e m m T l i ( 1 2 σ ) + ( 1 2 σ ) ( e i k T l k + e j i T l j ) + 3 e j k T l k l j l i ) d ω ( l )
Now note that we have e m m T l i = e j k T δ j k l i , e i k T l k = e j k T δ i j l k , e j i T l j = e j k T δ i k l j . Therefore we can write the equation for displacement u i C as
u i C = 1 8 π ( 1 σ ) 0 4 π r ( l ) ( e j k T δ j k l i ( 1 2 σ ) + ( 1 2 σ ) ( e j k T δ i j l k + e j k T δ i k l j ) + 3 e j k T l k l j l i ) d ω ( l )
Note that e j k T is constant throughout the inclusion volume. Therefore we take e j k outside of the integral
u i C = e j k T 8 π ( 1 σ ) 0 4 π r ( l ) ( δ j k l i ( 1 2 σ ) + ( 1 2 σ ) ( δ i j l k + ϵ j k δ i k l j ) + 3 l k l j l i ) d ω ( l )
u i C = e j k T 8 π ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( δ i j l k + ϵ j k δ i k l j δ j k l i ) + 3 l k l j l i ) d ω ( l )
which is the Equation (3.1) of Eshelby’s Classic Paper. To convert this into volume integral we note that the incremental solid angle is related to the incremental volume by the following relation
r ( l ) d ω ( l ) = d V r ( l ) 2
Therefore the equation for displacement u i C can be written as
u i C = e j k T 8 π ( 1 σ ) V ( ( 1 2 σ ) ( δ i j l k + ϵ j k δ i k l j δ j k l i ) + 3 l k l j l i ) r ( l ) 2 d V
which is the second part of the Equation (2.15) of Eshelby’s Classic Paper.
d x l ( d x k ( d x j ( d x i ( r ) ) ) ) = 3 δ i j x k x l r 5 5 x i x j x k x l r 7 + δ i l x j x k + δ j l x i x k + δ k l x i x j r 5

7. Strain due to Spontaneous Change of Form of Inclusion

In the previous section, we stated that the Displacement at r due to point force F i at r is (Love 1927)
U i ( r r ) = 1 4 π μ F i | r r | 1 16 π μ ( 1 σ ) F j 2 x j x i | r r |
In Stage III, we apply a force distribution F j = p j k T n k over S to make the body free of external force (but in a state of self-stress because of the transformation of the inclusion). The displacement impressed on the material in stage III is
u i C = S U i ( r r ) d S
u i C = 1 4 π μ S p j k T δ i j | r r | n k d S 1 16 π μ ( 1 σ ) S p j k T 2 x j x i | r r | n k d S
Since p j k T is constant throughout the surface S of the inclusion we can take p j k T outside the integral. Therefore we can write
u i C = p i k T 4 π μ S 1 | r r | n k d S p j k T 16 π μ ( 1 σ ) S 2 x j x i | r r | n k d S
We have earlier proved using the Gauss Divergence theorem that the 2 integrals can be written as
S 1 | r r | n k d S = x k ( V d V | r r | )
S 2 x j x i | r r | n k d S = 3 x j x i x k ( V | r r | d V ) = 3 x i x j x k ( V | r r | d V )
Therefore the displacement impressed on the material in stage III can be written as
u i C = p j k T 16 π μ ( 1 σ ) ψ , i j k p i k T 4 π μ ϕ , k
where
ψ = V d V | r r |
ϕ = V | r r | d V
are the ordinary Newtonian potential and the biharmonic potential of attracting matter of unit density filling the volume V bounded by S. We also proved earlier that
2 ψ = 2 ϕ
4 ψ = 2 2 ϕ = 8 π , inside S 0 , outside S
The strain in the material on stage III can be written as
e i l C = 1 2 ( x l u i C + x i u l C )
e i l C = 1 2 ( x l ( p j k T 16 π μ ( 1 σ ) ψ , i j k p i k T 4 π μ ϕ , k ) + x i ( p j k T 16 π μ ( 1 σ ) ψ , l j k p l k T 4 π μ ϕ , k ) )
e i l C = p j k T 32 π μ ( 1 σ ) ( ψ , i j k l + ψ , l j k i ) p i k T 8 π μ ϕ , k l p l k T 8 π μ ϕ , k i
e i l C = p j k T 16 π μ ( 1 σ ) ψ , i j k l p i k T 8 π μ ϕ , k l p l k T 8 π μ ϕ , k i
Therefore the dilatation in the material shall be
e i i C = p j k T 16 π μ ( 1 σ ) ψ , i j k i p i k T 8 π μ ϕ , k i p i k T 8 π μ ϕ , k i
e i i C = p j k T 16 π μ ( 1 σ ) ψ , i i j k p i k T 4 π μ ϕ , k i
Now we know that
2 ψ = 2 ϕ
2 x i x i ψ = 2 ϕ
2 x j x k ( 2 x i x i ψ ) = 2 2 x j x k ϕ
4 x i x i x j x k ψ = 2 2 x j x k ϕ
ψ i i j k l = 2 ϕ j k
Therefore the dilatation in the material shall be
e i i C = p j k T 16 π μ ( 1 σ ) ψ , i i j k p i k T 4 π μ ϕ , k i = p j k T 8 π μ ( 1 σ ) ϕ , j k p i k T 4 π μ ϕ , k i
Due to the balance of angular momentum, the stress tensor is symmetric i.e. p i k T = p k i T . Therefore we can write the above equation as
e i i C = p j k T 8 π μ ( 1 σ ) ϕ , j k p i k T 4 π μ ϕ , k i = p j k T 8 π μ ( 1 σ ) ϕ , j k p k i T 4 π μ ϕ , k i
Repeated indices i , k on the second term can be given in any other name like k , j . Therefore we can write the above equation as
e i i C = p j k T 8 π μ ( 1 σ ) ϕ , j k p j k T 4 π μ ϕ , j k = p j k T 8 π μ ( 1 σ ) ϕ , j k ( 1 2 ( 1 σ ) )
e i i C = p j k T 8 π μ ( 1 σ ) ϕ , j k p j k T 4 π μ ϕ , j k = p j k T 8 π μ ( 1 σ ) ϕ , j k ( 1 2 σ )
e i i C = p j k T 8 π μ ( 1 σ ) ϕ , j k ( 1 2 σ )
We have earlier derived that the strain in the material in stage III shall be
e i l C = p j k T 16 π μ ( 1 σ ) ψ , i j k l p i k T 8 π μ ϕ , k l p l k T 8 π μ ϕ , k i
Now note that for isotropic materials, we have
p j k T = λ e T δ j k + 2 μ e j k T , p i k T = λ e T δ i k + 2 μ e i k T , p l k T = λ e T δ l k + 2 μ e l k T
Substituting the above relations in the equation for the strain in the material in stage III shall be
e i l C = ( λ e T δ j k + 2 μ e j k T ) 16 π μ ( 1 σ ) ψ , i j k l ( λ e T δ i k + 2 μ e i k T ) 8 π μ ϕ , k l ( λ e T δ l k + 2 μ e l k T ) 8 π μ ϕ , k i
e i l C = ( λ e T δ j k 16 π μ ( 1 σ ) ψ , i j k l + 2 μ e j k T 16 π μ ( 1 σ ) ψ , i j k l ) ( λ e T δ i k 8 π μ ϕ , k l + 2 μ e i k T 8 π μ ϕ , k l ) ( λ e T δ l k 8 π μ ϕ , k i + 2 μ e l k T 8 π μ ϕ , k i )
e i l C = ( λ e T 16 π μ ( 1 σ ) ψ , i j j l + e j k T 8 π ( 1 σ ) ψ , i j k l ) ( λ e T 8 π μ ϕ , i l + e i k T 4 π ϕ , k l ) ( λ e T 8 π μ ϕ , l i + e l k T 4 π ϕ , k i )
e i l C = ( λ e T 16 π μ ( 1 σ ) ψ , j j i l + e j k T 8 π ( 1 σ ) ψ , i j k l ) ( λ e T 8 π μ ϕ , i l + λ e T 8 π μ ϕ , l i ) ( e i k T 4 π ϕ , k l + e l k T 4 π ϕ , k i )
Now we know that ψ , j j i l = 2 ϕ , i l and ϕ , i l = ϕ , l i , we can therefore the equation as
e i l C = ( λ e T 8 π μ ( 1 σ ) ϕ , i l + e j k T 8 π ( 1 σ ) ψ , i j k l ) λ e T 4 π μ ϕ , i l ( e i k T 4 π ϕ , k l + e l k T 4 π ϕ , k i )
e i l C = ( λ e T 8 π μ ( 1 σ ) ϕ , i l λ e T 4 π μ ϕ , i l ) ( e i k T 4 π ϕ , k l + e l k T 4 π ϕ , k i e j k T 8 π ( 1 σ ) ψ , i j k l ) )
e i l C = ( λ e T 8 π μ ( 1 σ ) ϕ , i l ( 1 2 ( 1 σ ) ) ) ( e i k T 4 π ϕ , k l + e l k T 4 π ϕ , k i e j k T 8 π ( 1 σ ) ψ , i j k l ) )
e i l C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i l ) ( e i k T 4 π ϕ , k l + e l k T 4 π ϕ , k i e j k T 8 π ( 1 σ ) ψ , i j k l ) )
e i l C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i l ) ( e i k T 4 π ϕ , k l + e l k T 4 π ϕ , k i e j k T 8 π ( 1 σ ) ψ , i j k l ) )
If e T is a pure dilatation 1 3 e T δ i j , then we have
e i k T = 1 3 e T δ i k , e l k T = 1 3 e T δ l k , e j k T = 1 3 e T δ j k
Therefore the above boxed equation can be written as
e i l C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i l ) ( e T δ i k 12 π ϕ , k l + e T δ l k 12 π ϕ , k i e T δ j k 24 π ( 1 σ ) ψ , i j k l ) )
e i l C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i l ) ( e T 12 π ϕ , i l + e T 12 π ϕ , l i e T 24 π ( 1 σ ) ψ , i j j l ) )
Since ϕ , l i = ϕ , i l , we can write the above equation as
e i l C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i l ) ( e T 6 π ϕ , i l e T 24 π ( 1 σ ) ψ , j j i l ) )
Now we know that ψ , j j i l = 2 ϕ , i l , we can therefore the above equation as
e i l C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i l ) ( e T 6 π ϕ , i l e T 12 π ( 1 σ ) ϕ , i l ) )
e i l C = ( ( λ ( 2 σ 1 ) 8 π μ ( 1 σ ) ) ( 1 6 π 1 12 π ( 1 σ ) ) ) e T ϕ , i l
e i l C = ( ( 3 λ ( 2 σ 1 ) 4 μ ( 1 σ ) + 2 μ ) 24 π μ ( 1 σ ) ) e T ϕ , i l = ( ( 3 λ ( 2 σ 1 ) 2 μ + 4 μ σ ) 24 π μ ( 1 σ ) ) e T ϕ , i l
Now the first lame constant λ can be written in terms of second lame constant μ and poisson’s ratio σ as
λ = 2 μ σ 1 2 σ
λ ( 2 σ 1 ) = 2 μ σ
3 λ ( 2 σ 1 ) = 6 μ σ
3 λ ( 2 σ 1 ) 2 μ + 4 μ σ = 6 μ σ 2 μ + 4 μ σ = 2 μ 2 μ σ = 2 μ ( 1 + σ )
3 λ ( 2 σ 1 ) 2 μ + 4 μ σ 24 π μ ( 1 σ ) = 2 μ ( 1 + σ ) 24 π μ ( 1 σ )
3 λ ( 2 σ 1 ) 2 μ + 4 μ σ 24 π μ ( 1 σ ) = ( 1 + σ ) 12 π ( 1 σ )
Therefore the strain e i l C an be written as
e i l C = ( ( 3 λ ( 2 σ 1 ) 4 μ ( 1 σ ) + 2 μ ) 24 π μ ( 1 σ ) ) e T ϕ , i l = ( ( 3 λ ( 2 σ 1 ) 2 μ + 4 μ σ ) 24 π μ ( 1 σ ) ) e T ϕ , i l = ( 1 + σ ) 12 π ( 1 σ ) e T ϕ , i l
e i l C = 1 4 π ( 1 + σ ) 3 ( 1 σ ) e T ϕ , i l

8. Discontinuities Across Inclusion Interface

The second derivative of a scalar function U in a given direction, with the direction specified by a unit vector n = ( n i , n j , n k ) , can be calculated by applying the directional derivative operator twice in that direction. Given that the direction vector is n = n i i ^ + n j j ^ + n k k ^ , the first derivative of U in the direction of n is given by:
U n = n i U x i + n j U x j + n k U x k
To find the second derivative, we take the derivative of the above expression again in the direction of n :
2 U n 2 = n n i U x i + n j U x j + n k U x k
Applying the directional derivative operator again:
2 U n 2 = n i n i 2 U x i 2 + n j 2 U x i x j + n k 2 U x i x k + n j n i 2 U x j x i + n j 2 U x j 2 + n k 2 U x j x k + n k n i 2 U x k x i + n j 2 U x k x j + n k 2 U x k 2
This can be simplified as:
2 U n 2 = p = 1 3 q = 1 3 n p n q 2 U x p x q
Where p and q run over the components i , j , k . This expression gives the second derivative of the scalar function U in the direction of the vector n . In indicial notation, the above formula can be written as
2 U n 2 = 2 U x i x j n i n j
To prove the theorem, we proceed step by step, ensuring that each part of the proof is mathematically rigorous. The theorem can be formally stated as:
Theorem (Poincaré 1899): Let U ( r ) be a potential function that satisfies Poisson’s equation in three-dimensional space,
2 U ( r ) = 4 π ρ ( r )
where ρ ( r ) is the mass density. If there is a surface S with normal n across which the density ρ ( r ) undergoes a discontinuity Δ ρ , then the second derivatives of the potential U ( r ) undergo a jump discontinuity given by
2 U x i x j = 4 π n i n j Δ ρ ,
where the bracket [ · ] denotes the difference in the second derivative as one crosses the surface S. Given Poisson’s equation:
2 U ( r ) = 4 π ρ ( r ) ,
where ρ ( r ) is the mass density distribution, the potential U ( r ) is related to the density ρ ( r ) by the integral:
U ( r ) = ρ ( r ) | r r | d 3 r .
Suppose that the density ρ ( r ) has a discontinuity across a surface S. Let ρ + and ρ denote the density values on either side of the surface. The discontinuity in ρ is given by:
Δ ρ = ρ + ρ .
Consider the surface S with a unit normal vector n at each point. We need to investigate the behavior of the second derivatives of U ( r ) as r crosses the surface S. The second derivative of U ( r ) in the direction n i n j is:
2 U x i x j = x j U x i .
To find the jump condition, we apply the divergence theorem over a small pillbox-shaped volume V that straddles the surface S. The pillbox has faces parallel and perpendicular to the surface S, with height ϵ in the direction normal to S. Integrating Poisson’s equation over V, we have:
V 2 U d V = 4 π V ρ d V
Using the divergence theorem, the left-hand side can be rewritten as:
V 2 U d V = V U · d A
where d A is the area element on the boundary V of the volume V. Let’s do the Evaluation on the Faces of the Pillbox. The contribution to the surface integral comes from the faces of the pillbox parallel to the surface S. Let the area of these faces be A, and the small height of the pillbox in the normal direction be ϵ . Then, the surface integral is dominated by the contributions from the faces perpendicular to n :
V U · d A = A U n ,
where U n is the derivative of U in the direction of n . Let us find the Relation Between the Jump in Density and the Jump in the Second Derivative. Now, consider the contribution of the jump in ρ to the integral:
V ρ d V = A ϵ ρ + + ρ 2 = A ϵ Δ ρ .
This leads to:
A U n = 4 π A ϵ Δ ρ .
Simplifying, we obtain:
2 U n 2 = 4 π Δ ρ
We have earlier proved that the second derivatives in arbitrary directions x i and x j are related to 2 U n 2 using the following relation
2 U n 2 = 2 U x i x j
Therefore the above boxed equation can be written as:
2 U x i x j = 4 π n i n j Δ ρ
where n i and n j are the components of the normal vector to the surface S. We just proved the Poincare 1899 theorem which states that the second derivatives of a potential function satisfying 2 U = 4 π ρ will have a jump
Δ U , i j = 4 π Δ ρ n i n j
on crossing a surface with normal n i across which there is a change of density Δ ρ . Let us consider the Potential function U = ϕ . Note that 2 U = 2 ϕ = 4 π inside S and 0 outside S. There is a jump of density Δ ρ = 1 across the surface S. Therefore we can write using the Poincare 1899 theorem that
ϕ , i j ( out ) ϕ , i j ( in ) = 4 π n i n j
Let us consider the Potential function U = ψ , i j . Note that 2 U = 2 ψ , i j = 2 ϕ , i j = 4 π ( ϕ , i j 2 π ) . Therefore the Potential function U = ψ , i j satisfies 2 U = 4 π ρ with ρ = ϕ , i j 2 π . There is a jump of density
Δ ρ = ( ϕ , i j ( in ) ϕ , i j ( out ) ) 2 π = 4 π n i n j 2 π = 2 n i n j
across the surface S. Therefore we can write using the Poincare 1899 theorem for the Potential function U = ψ , i j that
ψ , i j k l ( out ) ϕ , i j k l ( in ) = 4 π Δ ρ n k n l = 8 π n i n j n k n l
ψ , i j k l ( out ) ϕ , i j k l ( in ) = 8 π n i n j n k n l
We have earlier derived that the strain in the material at stage III shall be
e i l C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i l ) ( e i k T 4 π ϕ , k l + e l k T 4 π ϕ , k i e j k T 8 π ( 1 σ ) ψ , i j k l ) )
The trace of the strain tensor shall be
e i i C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i i ) ( e i k T 4 π ϕ , k i + e i k T 4 π ϕ , k i e j k T 8 π ( 1 σ ) ψ , i j k i ) )
e i i C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i i ) ( e i k T 2 π ϕ , k i e j k T 8 π ( 1 σ ) ψ , i j k i ) )
Since ϕ k i = ϕ i k and ψ i j k i = ψ i i j k = 2 ϕ j k , the trace of the strain tensor shall be
e i i C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i i ) ( e i k T 2 π ϕ , i k e j k T 4 π ( 1 σ ) ϕ , j k ) )
Repeated indices i , k in the second term can be renamed as j , k . Therefore we can write the equation as
e i i C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i i ) ( e j k T 2 π ϕ , j k e j k T 4 π ( 1 σ ) ϕ , j k ) )
e i i C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i i ) e j k T ϕ , j k ( 1 2 π 1 4 π ( 1 σ ) ) ) = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i i ) e j k T ϕ , j k ( 2 ( 1 σ ) 1 4 π ( 1 σ ) ) )
e i i C = ( λ ( 2 σ 1 ) 8 π μ ( 1 σ ) e T ϕ , i i ) e j k T ϕ , j k ( 1 2 σ 4 π ( 1 σ ) )
Now the first lame constant λ can be written in terms of second lame constant μ and poisson’s ratio σ as
λ = 2 μ σ 1 2 σ
λ ( 2 σ 1 ) = 2 μ σ
λ ( 2 σ 1 ) 8 π μ ( 1 σ ) = 2 μ σ 8 π μ ( 1 σ ) = σ 4 π ( 1 σ )
λ ( 2 σ 1 ) 8 π μ ( 1 σ ) = σ 4 π ( 1 σ )
Substituting the above equation in the trace of strain equation we get
e i i C = ( σ 4 π ( 1 σ ) e T ϕ , i i ) e j k T ϕ , j k ( 1 2 σ 4 π ( 1 σ ) )
We know that there is jump for both ϕ , i i and ϕ , j k (according to Poncaire’s 1899 Theorem). The jumps are
Δ ϕ , i i = ϕ , i i ( out ) ϕ , i i ( in ) = 4 π
Δ ϕ , j k = ϕ , j k ( out ) ϕ , j k ( in ) = 4 π n j n k
Substituting the above equations in the previous equation for the trace of strain we get
e i i C ( out ) e i i C ( in ) = ( σ 4 π ( 1 σ ) e T ) ( ϕ , i i ( out ) ϕ , i i ( in ) ) e j k T ( ϕ , j k ( out ) ϕ , j k ( in ) ) ( 1 2 σ 4 π ( 1 σ ) )
e i i C ( out ) e i i C ( in ) = ( σ ( 1 σ ) e T ) e j k T n j n k ( 1 2 σ ( 1 σ ) )
The strain tensor can be written in terms of hydrostatic and deviatoric components
e j k T = 1 3 e T δ j k + e j k T
Substituting the above equation into the previous boxed equation we get
e i i C ( out ) e i i C ( in ) = ( σ ( 1 σ ) e T ) ( 1 3 e T δ j k + e j k T ) n j n k ( 1 2 σ ( 1 σ ) )
e i i C ( out ) e i i C ( in ) = ( σ ( 1 σ ) e T ) ( 1 3 e T δ j k ) n j n k ( 1 2 σ ( 1 σ ) ) e j k T n j n k ( 1 2 σ ( 1 σ ) )
e i i C ( out ) e i i C ( in ) = ( σ ( 1 σ ) e T ) ( 1 3 e T ) n j n j ( 1 2 σ ( 1 σ ) ) e j k T n j n k ( 1 2 σ ( 1 σ ) )
e i i C ( out ) e i i C ( in ) = ( 3 σ ( 1 σ ) e T 3 ) ( e T 3 ) ( 1 2 σ ( 1 σ ) ) e j k T n j n k ( 1 2 σ ( 1 σ ) )
e i i C ( out ) e i i C ( in ) = e T 3 ( 3 σ ( 1 σ ) + 1 2 σ ( 1 σ ) ) e j k T n j n k ( 1 2 σ ( 1 σ ) )
e i i C ( out ) e i i C ( in ) = e T 3 ( 1 + σ 1 σ ) e j k T n j n k ( 1 2 σ 1 σ )
We have earlier derived that the strain in the material at stage III shall be
e i l C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i l ) ( e i k T 4 π ϕ , k l + e l k T 4 π ϕ , k i e j k T 8 π ( 1 σ ) ψ , i j k l ) )
We know that there is jump for both ϕ , i l , ϕ , k l , ϕ , k i and ψ i j k l (all 4 of them are according to Poncaire’s 1899 Theorem). The jumps are
Δ ϕ , i l = ϕ , i l ( out ) ϕ , i l ( in ) = 4 π n i n l
Δ ϕ , k l = ϕ , k l ( out ) ϕ , k l ( in ) = 4 π n k n l
Δ ϕ , k i = ϕ , k i ( out ) ϕ , k i ( in ) = 4 π n k n i
Δ ψ , i j k l = ϕ , i j k l ( out ) ϕ , i j k l ( in ) = 8 π n i n j n k n l
Therefore the jump in e i l C shall be
Δ e i l C = e i l C ( out ) e i l C ( in ) = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ( ϕ , i l ( out ) ϕ , i l ( in ) ) ) ( e i k T 4 π ( ϕ , k l ( out ) ϕ , k l ( in ) ) + e l k T 4 π ( ϕ , k i ( out ) ϕ , k i ( in ) )
+ e j k T 8 π ( 1 σ ) ( ψ , i j k l ( out ) ψ , i j k l ( in ) ) )
Substituting the boxed equations into the previous equations we get
Δ e i l C = e i l C ( out ) e i l C ( in ) = λ ( 2 σ 1 ) 2 μ ( 1 σ ) e T n i n l e i k T n k n l e l k T n k n i + e j k T ( 1 σ ) n i n j n k n l
We have earlier derived that
λ ( 2 σ 1 ) 8 π μ ( 1 σ ) = σ 4 π ( 1 σ )
λ ( 2 σ 1 ) 2 μ ( 1 σ ) = σ ( 1 σ )
Therefore we have
Δ e i l C = e i l C ( out ) e i l C ( in ) = σ ( 1 σ ) e T n i n l e i k T n k n l e l k T n k n i e j k T ( 1 σ ) n i n j n k n l
The strain tensor can be written in terms of hydrostatic and deviatoric components
e i k T = 1 3 e T δ i k + e i k T
e l k T = 1 3 e T δ l k + e l k T
e j k T = 1 3 e T δ j k + e j k T
Therefore we can write
Δ e i l C = e i l C ( out ) e i l C ( in ) = σ ( 1 σ ) e T n i n l ( 1 3 e T δ i k + e i k T ) n k n l ( 1 3 e T δ l k + e l k T ) n k n i + ( 1 3 e T δ j k + e j k T ) ( 1 σ ) n i n j n k n l
Δ e i l C = e i l C ( out ) e i l C ( in ) = σ ( 1 σ ) e T n i n l 1 3 e T δ i k n k n l e i k T n k n l 1 3 e T δ l k n k n i e l k T n k n i + e T δ j k 3 ( 1 σ ) n i n j n k n l + e j k T ( 1 σ ) n i n j n k n l
Δ e i l C = e i l C ( out ) e i l C ( in ) = σ ( 1 σ ) e T n i n l 1 3 e T n i n l e i k T n k n l 1 3 e T n l n i e l k T n k n i + e T 3 ( 1 σ ) n i n j n j n l + e j k T ( 1 σ ) n i n j n k n l
Δ e i l C = e i l C ( out ) e i l C ( in ) = σ ( 1 σ ) e T n i n l 1 3 e T n i n l e i k T n k n l 1 3 e T n l n i e l k T n k n i + e T 3 ( 1 σ ) n i n l + e j k T ( 1 σ ) n i n j n k n l
Δ e i l C = e i l C ( out ) e i l C ( in ) = e T n i n l ( σ ( 1 σ ) + 2 3 1 3 ( 1 σ ) ) e i k T n k n l e l k T n k n i + e j k T ( 1 σ ) n i n j n k n l
Δ e i l C = e i l C ( out ) e i l C ( in ) = e T n i n l ( ( 3 σ + 2 ( 1 σ ) 1 ) 3 ( 1 σ ) ) e i k T n k n l e l k T n k n i + e j k T ( 1 σ ) n i n j n k n l
Δ e i l C = e i l C ( out ) e i l C ( in ) = ( 1 + σ ) 3 ( 1 σ ) e T n i n l e i k T n k n l e l k T n k n i + e j k T ( 1 σ ) n i n j n k n l
Now note that Δ e i l C shall be
Δ e i l C = Δ ( 1 3 e C δ i l + e i l C ) = 1 3 δ i l Δ e C + Δ e i l C
We have earlier derived that Δ e C shall be
Δ e C = e i i C ( out ) e i i C ( in ) = e T 3 ( 1 + σ 1 σ ) e j k T n j n k ( 1 2 σ 1 σ )
Therefore Δ e i l C shall be
Δ e i l C = 1 3 δ i l ( e T 3 ( 1 + σ 1 σ ) e j k T n j n k ( 1 2 σ 1 σ ) ) + Δ e i l C
Δ e i l C = e T 9 ( 1 + σ 1 σ ) δ i l e j k T n j n k δ i l ( 1 2 σ 3 ( 1 σ ) ) + Δ e i l C
Therefore we can write the change in the deviatoric component of the strain tensor Δ e i l C as
Δ e i l C = Δ e i l C + e T 9 ( 1 + σ 1 σ ) δ i l + e j k T n j n k δ i l ( 1 2 σ 3 ( 1 σ ) )
Δ e i l C = ( 1 + σ ) 3 ( 1 σ ) e T n i n l e i k T n k n l e l k T n k n i + e j k T ( 1 σ ) n i n j n k n l + e T 9 ( 1 + σ 1 σ ) δ i l + e j k T n j n k δ i l ( 1 2 σ 3 ( 1 σ ) )
Δ e i l C = ( 1 + σ ) 3 ( 1 σ ) e T ( n i n l 1 3 δ i l ) e i k T n k n l e l k T n k n i + 1 ( 1 σ ) e j k T n i n j n k n l + 1 2 σ 3 ( 1 σ ) e j k T n j n k δ i l

9. Elastic Field in a Spherical and Ellipsoidal Inclusion

In this section, we shall analyze the elastic field due to the spontaneous change of form of an ellipsoidal inclusion within an isotropic elastic solid. We assume that the equation of ellipsoid is
X 2 a 2 + Y 2 b 2 + Z 2 c 2 = 1
Let’s assume a point inside the r = ( x , y , z ) within this ellipsoidal inclusion. We define the distance of the point r to the incremental surface d S in the direction l = ( l 1 , l 2 , l 3 ) = ( l , m , n ) as r ( l ) . Therefore r ( l ) is the positive root of
( x + r ( l ) l ) 2 a 2 + ( y + r ( l ) m ) 2 b 2 + ( z + r ( l ) n ) 2 c 2 = 1
( x 2 + r ( l ) 2 l 2 + 2 x r ( l ) l ) a 2 + ( y 2 + r ( l ) 2 m 2 + 2 y r ( l ) m ) b 2 + ( z 2 + r ( l ) 2 n 2 + 2 z r ( l ) n ) c 2 = 1
r ( l ) 2 ( l 2 a 2 + m 2 b 2 + n 2 c 2 ) + 2 r ( l ) ( x l a 2 + y m b 2 + z n c 2 ) ( 1 x 2 a 2 y 2 b 2 z 2 c 2 ) = 0
Let us define the following quantities
g = l 2 a 2 + m 2 b 2 + n 2 c 2
f = x l a 2 + y m b 2 + z n c 2
e = 1 x 2 a 2 y 2 b 2 z 2 c 2
The above boxed equation can therefore be written as
g r ( l ) 2 + 2 f r ( l ) e = 0
Therefore the solution of the above boxed equation shall be
r ( l ) = 2 f ± 4 f 2 + 4 g e 2 g
r ( l ) = f g ± f 2 g 2 + e g
In the earlier section we derived that
u i C = e j k T 8 π ( 1 σ ) 0 4 π r ( l ) ( ( 1 2 σ ) ( δ i j l k + ϵ j k δ i k l j δ j k l i ) + 3 l k l j l i ) d ω ( l ) = e j k T 8 π ( 1 σ ) 0 4 π r ( l ) g i j k ( l ) d ω ( l )
Substituting the expression r ( l ) = f g ± f 2 g 2 + e g in the above equation we get
u i C = e j k T 8 π ( 1 σ ) 0 4 π ( f g ± f 2 g 2 + e g ) g i j k ( l ) d ω ( l ) = e j k T 8 π ( 1 σ ) 0 4 π f g i j k ( l ) g d ω ( l ) ± e j k T 8 π ( 1 σ ) 0 4 π ( f 2 g 2 + e g ) g i j k ( l ) d ω ( l )
Note that ( f ( l ) 2 g ( l ) 2 + e ( l ) g ( l ) ) is an even function. To prove that we just need to show that ( f ( l ) 2 g ( l ) 2 + e ( l ) g ( l ) ) = ( f ( l ) 2 g ( l ) 2 + e ( l ) g ( l ) ) . That is easy to show once we realize that f ( l ) = f ( l ) , hence f ( l ) 2 = f ( l ) 2 . We can also easily show that g ( l ) = g ( l ) , e ( l ) = e ( l ) and hence g ( l ) 2 = g ( l ) 2 . But g i j k ( l ) is an odd function. Note that Integration of the Product of an even function ( f ( l ) 2 g ( l ) 2 + e ( l ) g ( l ) ) and odd function g i j k ( l ) shall always be zero. Therefore we can write
e j k T 8 π ( 1 σ ) 0 4 π ( f 2 g 2 + e g ) g i j k ( l ) d ω ( l ) = 0
Therefore we have
u i C = e j k T 8 π ( 1 σ ) 0 4 π f g i j k ( l ) g d ω ( l )
u i C = e j k T 8 π ( 1 σ ) 0 4 π f g i j k ( l ) g d ω ( l )
Defining λ 1 = l a 2 , λ 2 = m b 2 , λ 3 = n c 2 we can write f ( l ) as follows
f ( l ) = x l a 2 + y m b 2 + z n c 2 = x m λ m
Therefore we can write the above boxed equation as
u i C = e j k T 8 π ( 1 σ ) 0 4 π x m λ m g i j k ( l ) g d ω ( l )
Since the point r ( l ) = ( x , y , z ) inside the volume V is a fixed point, x m can be taken out of the integral. Therefore, we can write the above equation as
u i C = x m e j k T 8 π ( 1 σ ) 0 4 π λ m g i j k ( l ) g d ω ( l )
u l C = x m e j k T 8 π ( 1 σ ) 0 4 π λ m g l j k ( l ) g d ω ( l )
The strain e i l C shall be therefore
e i l C = 1 2 ( u i C x l + u l C x i ) = δ m l e j k T 16 π ( 1 σ ) 0 4 π λ m g i j k ( l ) g d ω ( l ) + δ m i e j k T 16 π ( 1 σ ) 0 4 π λ m g l j k ( l ) g d ω ( l )
e i l C = e j k T 16 π ( 1 σ ) 0 4 π ( λ l g i j k ( l ) + λ i g l j k ( l ) ) g d ω ( l )
We can write the relation between the constrained and stress-free strains in the inclusion as
e i l C = S i l j k e j k T
where S i l j k is the Eshelby Tensor defined as
S i l j k = 1 16 π ( 1 σ ) 0 4 π ( λ l g i j k ( l ) + λ i g l j k ( l ) ) g d ω ( l )

9.1. Spherical Inclusion

Let’s first compute the for a spherical inclusion where we have a = b = c
S i l j k = 1 16 π ( 1 σ ) 0 4 π ( λ l g i j k ( l ) + λ i g l j k ( l ) ) g d ω ( l ) = 1 16 π ( 1 σ ) 0 4 π ( l l g i j k ( l ) + l i g l j k ( l ) ) a 2 g d ω ( l )
Now note that we have
g i j k = ( 1 2 σ ) ( δ i j l k + δ i k l j δ j k l i ) + 3 l i l j l k
g l j k = ( 1 2 σ ) ( δ l j l k + δ l k l j δ j k l l ) + 3 l l l j l k
Therefore we can write
S i l j k = 1 16 π ( 1 σ ) 0 4 π ( l l ( ( 1 2 σ ) ( δ i j l k + δ i k l j δ j k l i ) + 3 l i l j l k ) + l i ( ( 1 2 σ ) ( δ l j l k + δ l k l j δ j k l l ) + 3 l l l j l k ) ) a 2 g d ω ( l )
S i l j k = 1 16 π ( 1 σ ) 0 4 π ( ( 1 2 σ ) ( δ i j l k l l + δ i k l j l l + δ l j l k l i + δ l k l j l i 2 δ j k l i l l ) + 6 l i l j l k l l ) a 2 g d ω ( l )
S i l j k = ( 1 2 σ ) 16 π ( 1 σ ) ( δ i j 0 4 π l k l l a 2 g d ω ( l ) + δ i k 0 4 π l j l l a 2 g d ω ( l ) + δ l j 0 4 π l k l i a 2 g d ω ( l ) + δ l k 0 4 π l j l i a 2 g d ω ( l ) 2 δ j k 0 4 π l i l l a 2 g d ω ( l ) )
+ 3 8 π ( 1 σ ) 0 4 π l i l j l k l l a 2 g d ω ( l )
Using the Routh Integrals we can reduce the solid angle integral 0 4 π l k l l a 2 g d ω ( l ) to simple integrals
0 4 π l k l l a 2 g d ω ( l ) = δ k l 0 4 π l l 2 a 2 g d ω ( l ) = δ k l ( 2 π a b c 0 d u ( a 2 + u ) Δ )
0 4 π l k l l a 2 g d ω ( l ) = δ k l ( 2 π a b c 0 d u ( a 2 + u ) Δ )
where Δ is given by the expression
Δ = ( a 2 + u ) 1 2 ( b 2 + u ) 1 2 ( c 2 + u ) 1 2
For a spherical inclusion where we have a = b = c , the Δ shall be
Δ = ( a 2 + u ) 3 2
Therefore the integral shall be
0 d u ( a 2 + u ) Δ = 0 d u ( a 2 + u ) 5 2
To solve the integral 0 1 ( a 2 + u ) 5 / 2 d u , we can use a standard technique for integrals of this form. Let’s set u = a 2 tan 2 ( θ ) . Then, d u = 2 a 2 tan ( θ ) sec 2 ( θ ) d θ . Substitute this into the integral:
0 1 ( a 2 + u ) 5 / 2 d u = 0 π 2 1 ( a 2 + a 2 tan 2 ( θ ) ) 5 / 2 · 2 a 2 tan ( θ ) sec 2 ( θ ) d θ
The term inside the integral simplifies as follows:
a 2 ( 1 + tan 2 ( θ ) ) = a 2 sec 2 ( θ )
Thus, the integral becomes:
0 1 ( a 2 + u ) 5 / 2 d u = 0 π 2 2 a 2 tan ( θ ) sec 2 ( θ ) ( a 2 sec 2 ( θ ) ) 5 / 2 d θ = 0 π 2 2 a 2 tan ( θ ) sec 2 ( θ ) a 5 sec 5 ( θ ) d θ = 2 a 3 0 π 2 tan ( θ ) sec 3 ( θ ) d θ = 2 a 3 0 π 2 sin ( θ ) cos 2 ( θ ) d θ
We can evaluate the integral 0 π 2 sin ( θ ) cos 2 ( θ ) d θ using the substitution x = cos ( θ ) , d x = sin ( θ ) d θ :
= 2 a 3 1 0 x 2 d x = 2 a 3 x 3 3 1 0 = 2 a 3 · 1 3 = 2 3 a 3
The value of the integral is:
0 1 ( a 2 + u ) 5 / 2 d u = 2 3 a 3
Therefore the solid angle integral 0 4 π l k l l a 2 g d ω ( l ) shall be
0 4 π l k l l a 2 g d ω ( l ) = δ k l ( 2 π a b c 0 d u ( a 2 + u ) Δ )
For spherical inclusion a = b = c , therefore the solid angle integral 0 4 π l k l l a 2 g d ω ( l ) shall be
0 4 π l k l l a 2 g d ω ( l ) = δ k l ( 2 π a 3 0 d u ( a 2 + u ) Δ ) = δ k l ( 2 π a 3 2 3 a 3 ) = 4 π 3 δ k l
0 4 π l k l l a 2 g d ω ( l ) = 4 π 3 δ k l
Similarly, we shall have
0 4 π l j l l a 2 g d ω ( l ) = 4 π 3 δ j l
0 4 π l k l i a 2 g d ω ( l ) = 4 π 3 δ k i
0 4 π l j l i a 2 g d ω ( l ) = 4 π 3 δ j i
0 4 π l i l l a 2 g d ω ( l ) = 4 π 3 δ i l
Using the Routh Integrals we can reduce the solid angle integral 0 4 π l i l j l k l l a 2 g d ω ( l ) to simple integrals
0 4 π l i l j l k l l a 2 g d ω ( l ) = 1 3 ( δ i k δ j l + δ i l δ j k + δ i j δ l k ) a 2 0 4 π l l 4 a 4 g d ω ( l ) = 1 3 ( δ i k δ j l + δ i l δ j k + δ i j δ l k ) a 2 ( 2 π a b c 0 d u ( a 2 + u ) 2 Δ )
0 4 π l i l l l j l k a 2 g d ω ( l ) = 2 ( δ i k δ j l + δ i l δ j k + δ i j δ l k ) π a 3 b c 3 0 d u ( a 2 + u ) 2 Δ
where Δ is given by the expression
Δ = ( a 2 + u ) 1 2 ( b 2 + u ) 1 2 ( c 2 + u ) 1 2
For a spherical inclusion where we have a = b = c , the Δ shall be
Δ = ( a 2 + u ) 3 2
Therefore the integral shall be
0 d u ( a 2 + u ) 2 Δ = 0 d u ( a 2 + u ) 7 2
To solve the integral 0 1 ( a 2 + u ) 7 / 2 d u , we can follow a similar approach to the one used for the previous integral. Let u = a 2 tan 2 ( θ ) . Then, d u = 2 a 2 tan ( θ ) sec 2 ( θ ) d θ . Substituting into the integral:
0 1 ( a 2 + u ) 7 / 2 d u = 0 π 2 1 ( a 2 + a 2 tan 2 ( θ ) ) 7 / 2 · 2 a 2 tan ( θ ) sec 2 ( θ ) d θ
The expression simplifies as follows:
0 1 ( a 2 + u ) 7 / 2 d u = 0 π 2 2 a 2 tan ( θ ) sec 2 ( θ ) ( a 2 sec 2 ( θ ) ) 7 / 2 d θ = 0 π 2 2 a 2 tan ( θ ) sec 2 ( θ ) a 7 sec 7 ( θ ) d θ = 2 a 5 0 π 2 tan ( θ ) sec 5 ( θ ) d θ = 2 a 5 0 π 2 sin ( θ ) cos 4 ( θ ) d θ
We can evaluate the integral 0 π 2 sin ( θ ) cos 4 ( θ ) d θ using the substitution x = cos ( θ ) , d x = sin ( θ ) d θ :
= 2 a 5 1 0 x 4 d x = 2 a 5 x 5 5 1 0 = 2 a 5 · 1 5 = 2 5 a 5
The value of the integral is:
0 d u ( a 2 + u ) 2 Δ = 0 1 ( a 2 + u ) 7 / 2 d u = 2 5 a 5
Therefore 0 4 π l i l j l k l l a 2 g d ω ( l ) can be written as
0 4 π l i l j l k l l a 2 g d ω ( l ) = 2 ( δ i k δ j l + δ i l δ j k + δ i j δ l k ) π a 3 b c 3 0 d u ( a 2 + u ) 2 Δ = ( δ i k δ j l + δ i l δ j k + δ i j δ l k ) 4 π b c 15 a 2
For a spherical inclusion where we have a = b = c , therefore we have
0 4 π l i l j l k l l a 2 g d ω ( l ) = 4 π 15 ( δ i k δ j l + δ i l δ j k + δ i j δ l k )
Therefore the Eshelby Tensor  S i l j k shall be
S i l j k = ( 1 2 σ ) 16 π ( 1 σ ) ( δ i j 0 4 π l k l l a 2 g d ω ( l ) + δ i k 0 4 π l j l l a 2 g d ω ( l ) + δ l j 0 4 π l k l i a 2 g d ω ( l ) + δ l k 0 4 π l j l i a 2 g d ω ( l ) 2 δ j k 0 4 π l i l l a 2 g d ω ( l ) )
+ 3 8 π ( 1 σ ) 0 4 π l i l j l k l l a 2 g d ω ( l )
S i l j k = ( 1 2 σ ) 16 π ( 1 σ ) ( δ i j δ k l 4 π 3 + δ i k δ j l 4 π 3 + δ l j δ k i 4 π 3 + δ l k δ j i 4 π 3 2 δ j k δ i l 4 π 3 ) + 1 10 ( 1 σ ) ( δ i k δ j l + δ i l δ j k + δ i j δ l k )
S i l j k = ( 1 2 σ ) 16 π ( 1 σ ) ( δ i j δ k l 8 π 3 + δ i k δ j l 8 π 3 δ j k δ i l 8 π 3 ) + 1 10 ( 1 σ ) ( δ i k δ j l + δ i l δ j k + δ i j δ l k )
S i l j k = ( 1 2 σ ) 6 ( 1 σ ) ( δ i j δ k l + δ i k δ j l δ j k δ i l ) + 1 10 ( 1 σ ) ( δ i k δ j l + δ i l δ j k + δ i j δ l k )
S i l j k = ( δ i j δ k l + δ i k δ j l ) ( ( 1 2 σ ) 6 ( 1 σ ) + 1 10 ( 1 σ ) ) + δ i l δ j k ( 1 10 ( 1 σ ) ( 1 2 σ ) 6 ( 1 σ ) )
S i l j k = ( δ i j δ k l + δ i k δ j l ) ( 5 ( 1 2 σ ) + 3 30 ( 1 σ ) ) + δ i l δ j k ( 5 ( 1 2 σ ) + 3 30 ( 1 σ ) )
S i l j k = ( δ i j δ k l + δ i k δ j l ) ( 8 10 σ 30 ( 1 σ ) ) + δ i l δ j k ( 10 σ 2 30 ( 1 σ ) )
S i l j k = ( δ i j δ k l + δ i k δ j l ) ( 4 5 σ 15 ( 1 σ ) ) + δ i l δ j k ( 5 σ 1 15 ( 1 σ ) )
Notice that there is no term of a , b , c in this equation. We can therefore say that the Eshelby tensor does not depend on the radius of the spherical inclusion.

9.2. Ellipsoidal Inclusion

Let’s now compute the Eshelby tensor for Ellipsoidal Inclusion. The S 1111 , S 1122 , S 1133 , S 1212 , S 1112 , S 1223 , S 1232 components of the Eshelby Tensor shall be
S 1111 = 1 16 π ( 1 σ ) 0 4 π ( λ 1 g 111 ( l ) + λ 1 g 111 ( l ) ) g d ω ( l ) = 0 4 π ( l g 111 ( l ) + l g 111 ( l ) ) a 2 g d ω ( l )
S 1122 = 1 16 π ( 1 σ ) 0 4 π ( λ 1 g 122 ( l ) + λ 1 g 122 ( l ) ) g d ω ( l ) = 0 4 π ( l g 122 ( l ) + l g 122 ( l ) ) a 2 g d ω ( l )
S 1133 = 1 16 π ( 1 σ ) 0 4 π ( λ 1 g 133 ( l ) + λ 1 g 133 ( l ) ) g d ω ( l ) = 0 4 π ( l g 133 ( l ) + l g 133 ( l ) ) a 2 g d ω ( l )
S 1212 = 1 16 π ( 1 σ ) 0 4 π ( λ 2 g 112 ( l ) + λ 1 g 212 ( l ) ) g d ω ( l ) = 0 4 π ( m b 2 g 112 ( l ) + l a 2 g 212 ( l ) ) g d ω ( l )
S 1112 = 1 16 π ( 1 σ ) 0 4 π ( λ 1 g 112 ( l ) + λ 1 g 112 ( l ) ) g d ω ( l ) = 1 16 π ( 1 σ ) 0 4 π ( l g 112 ( l ) + l g 112 ( l ) ) a 2 g d ω ( l )
S 1223 = 1 16 π ( 1 σ ) 0 4 π ( λ 2 g 123 ( l ) + λ 1 g 223 ( l ) ) g d ω ( l ) = 1 16 π ( 1 σ ) 0 4 π ( m b 2 g 123 ( l ) + l a 2 g 223 ( l ) ) g d ω ( l )
S 1232 = 1 16 π ( 1 σ ) 0 4 π ( λ 2 g 132 ( l ) + λ 1 g 232 ( l ) ) g d ω ( l ) = 1 16 π ( 1 σ ) 0 4 π ( m b 2 g 132 ( l ) + l a 2 g 232 ( l ) ) g d ω ( l )
Now note that we have
g i j k = ( 1 2 σ ) ( δ i j l k + δ i k l j δ j k l i ) + 3 l i l j l k
Therefore the components of the 3rd order tensor g 111 , g 112 , g 122 , g 123 , g 132 , g 223 , g 232 , g 133 , g 212 shall be
g 111 = ( 1 2 σ ) ( δ 11 l 1 + δ 11 l 1 δ 11 l 1 ) + 3 l 1 l 1 l 1 = ( 1 2 σ ) l 1 + 3 l 1 3
g 111 = ( 1 2 σ ) l + 3 l 3
g 112 = ( 1 2 σ ) ( δ 11 l 2 + δ 12 l 1 δ 12 l 1 ) + 3 l 1 2 l 2 = ( 1 2 σ ) l 2 + 3 l 1 2 l 2
g 112 = ( 1 2 σ ) m + 3 l 2 m
g 122 = ( 1 2 σ ) ( δ 12 l 2 + δ 12 l 2 δ 22 l 1 ) + 3 l 1 l 2 l 2 = ( 1 2 σ ) l 1 + 3 l 1 l 2 2
g 122 = ( 1 2 σ ) l + 3 l m 2
g 123 = ( 1 2 σ ) ( δ 12 l 3 + δ 13 l 2 δ 23 l 1 ) + 3 l 1 l 2 l 3 = 3 l 1 l 2 l 3
g 123 = 3 l m n
g 132 = ( 1 2 σ ) ( δ 13 l 2 + δ 12 l 3 δ 32 l 1 ) + 3 l 1 l 3 l 2 = 3 l 1 l 2 l 3
g 132 = 3 l m n
g 223 = ( 1 2 σ ) ( δ 22 l 3 + δ 23 l 2 δ 23 l 2 ) + 3 l 2 l 2 l 3 = ( 1 2 σ ) l 3 + 3 l 2 l 2 l 3
g 223 = ( 1 2 σ ) n + 3 m 2 n
g 232 = ( 1 2 σ ) ( δ 23 l 2 + δ 22 l 3 δ 32 l 2 ) + 3 l 2 l 3 l 2 = ( 1 2 σ ) l 3 + 3 l 2 l 3 l 2
g 232 = ( 1 2 σ ) n + 3 m 2 n
g 133 = ( 1 2 σ ) ( δ 13 l 3 + δ 13 l 3 δ 33 l 1 ) + 3 l 1 l 3 l 3 = ( 1 2 σ ) l 1 + 3 l 1 l 3 2
g 133 = ( 1 2 σ ) l + 3 l n 2
g 212 = ( 1 2 σ ) ( δ 21 l 2 + δ 22 l 1 δ 12 l 2 ) + 3 l 2 l 1 l 2 = ( 1 2 σ ) l 1 + 3 l 1 l 2 2
g 212 = ( 1 2 σ ) l + 3 l m 2
Before we compute the components of the Eshelby Tensor we would like to state that all integrals of the form 0 4 π l i m j n k d ω g shall vanish if any one of the i , j , k are odd, i.e we have
0 4 π l i m j n k d ω g = 0 if any of the i , j , k is odd
Therefore the S 1111 components of the Eshelby Tensor shall be
S 1111 = 1 16 π ( 1 σ ) 0 4 π ( l g 111 ( l ) + l g 111 ( l ) ) a 2 g d ω ( l ) = 1 8 π ( 1 σ ) 0 4 π l g 111 ( l ) a 2 g d ω ( l ) = 1 8 π ( 1 σ ) 0 4 π l ( ( 1 2 σ ) l + 3 l 3 ) a 2 g d ω ( l )
S 1111 = 1 8 π ( 1 σ ) 0 4 π ( 1 2 σ ) l 2 a 2 g d ω ( l ) + 1 8 π ( 1 σ ) 0 4 π 3 l 4 a 2 g d ω ( l )
S 1111 = ( 1 2 σ ) 8 π ( 1 σ ) 0 4 π l 2 a 2 d ω ( l ) g + 3 a 2 8 π ( 1 σ ) 0 4 π l 4 a 4 d ω ( l ) g
Therefore the S 1122 components of the Eshelby Tensor shall be
S 1122 = 1 16 π ( 1 σ ) 0 4 π ( l g 122 ( l ) + l g 122 ( l ) ) a 2 g d ω ( l ) = 1 8 π ( 1 σ ) 0 4 π l g 122 ( l ) a 2 g d ω ( l ) = 1 8 π ( 1 σ ) 0 4 π l ( ( 1 2 σ ) l + 3 l m 2 ) a 2 g d ω ( l )
S 1122 = 1 8 π ( 1 σ ) 0 4 π ( 1 2 σ ) l 2 a 2 g d ω ( l ) + 1 8 π ( 1 σ ) 0 4 π 3 l 2 m 2 a 2 g d ω ( l )
S 1122 = ( 1 2 σ ) 8 π ( 1 σ ) 0 4 π l 2 a 2 d ω ( l ) g + 3 b 2 8 π ( 1 σ ) 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g
Therefore the S 1133 components of the Eshelby Tensor shall be
S 1133 = 1 16 π ( 1 σ ) 0 4 π ( l g 133 ( l ) + l g 133 ( l ) ) a 2 g d ω ( l ) = 1 8 π ( 1 σ ) 0 4 π l g 133 ( l ) a 2 g d ω ( l ) = 1 8 π ( 1 σ ) 0 4 π l ( ( 1 2 σ ) l + 3 l n 2 ) a 2 g d ω ( l )
S 1133 = 1 8 π ( 1 σ ) 0 4 π ( 1 2 σ ) l 2 a 2 g d ω ( l ) + 1 8 π ( 1 σ ) 0 4 π 3 l 2 n 2 a 2 g d ω ( l )
S 1133 = ( 1 2 σ ) 8 π ( 1 σ ) 0 4 π l 2 a 2 d ω ( l ) g + 3 c 2 8 π ( 1 σ ) 0 4 π l 2 a 2 n 2 c 2 d ω ( l ) g
Therefore the S 1212 component of the Eshelby Tensor shall be
S 1212 = 1 16 π ( 1 σ ) 0 4 π m g 112 ( l ) b 2 g d ω ( l ) + 1 16 π ( 1 σ ) 0 4 π l g 212 ( l ) a 2 g d ω ( l )
S 1212 = 1 16 π ( 1 σ ) 0 4 π m ( ( 1 2 σ ) m + 3 l 2 m ) b 2 g d ω ( l ) + 1 16 π ( 1 σ ) 0 4 π l ( ( 1 2 σ ) l + 3 l m 2 ) a 2 g d ω ( l )
S 1212 = ( 1 2 σ ) 16 π ( 1 σ ) 0 4 π m 2 b 2 d ω ( l ) g + 3 a 2 16 π ( 1 σ ) 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g + ( 1 2 σ ) 16 π ( 1 σ ) 0 4 π l 2 a 2 d ω ( l ) g + 3 b 2 16 π ( 1 σ ) 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g
S 1212 = ( 1 2 σ ) 16 π ( 1 σ ) 0 4 π m 2 b 2 d ω ( l ) g + 3 ( a 2 + b 2 ) 16 π ( 1 σ ) 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g + ( 1 2 σ ) 16 π ( 1 σ ) 0 4 π l 2 a 2 d ω ( l ) g
Therefore the S 1112 component of the Eshelby Tensor shall be
S 1112 = 1 16 π ( 1 σ ) 0 4 π ( l g 112 ( l ) + l g 112 ( l ) ) a 2 g d ω ( l ) = 1 8 π ( 1 σ ) 0 4 π l g 112 ( l ) a 2 g d ω ( l )
S 1112 = 1 8 π ( 1 σ ) 0 4 π l ( ( 1 2 σ ) m + 3 l 2 m ) a 2 g d ω ( l ) = ( 1 2 σ ) 8 π ( 1 σ ) 0 4 π l m a 2 g d ω ( l ) + 3 8 π ( 1 σ ) 0 4 π l 3 m a 2 g d ω ( l )
Now note that both the integrals on the RHS of the previous equation have l , m raised to odd powers. As we previously stated, all integrals of the form 0 4 π l i m j n k d ω g shall vanish if any one of the i , j , k are odd. Therefore we can write
0 4 π l m a 2 g d ω ( l ) = 0 , 0 4 π l 3 m a 2 g d ω ( l ) = 0
Therefore the S 1112 component of the Eshelby Tensor shall be
S 1112 = ( 1 2 σ ) 8 π ( 1 σ ) 0 4 π l m a 2 g d ω ( l ) + 3 8 π ( 1 σ ) 0 4 π l 3 m a 2 g d ω ( l ) = 0
S 1112 = 0
Therefore the S 1223 component of the Eshelby Tensor shall be
S 1223 = 1 16 π ( 1 σ ) 0 4 π ( λ 2 g 123 ( l ) + λ 1 g 223 ( l ) ) g d ω ( l ) = 1 16 π ( 1 σ ) 0 4 π ( m b 2 g 123 ( l ) + l a 2 g 223 ( l ) ) g d ω ( l )
S 1223 = 1 16 π ( 1 σ ) 0 4 π m g 123 ( l ) b 2 g d ω ( l ) + 1 16 π ( 1 σ ) 0 4 π l g 223 ( l ) ) a 2 g d ω ( l )
S 1223 = 3 16 π ( 1 σ ) 0 4 π l m 2 n b 2 g d ω ( l ) + 1 16 π ( 1 σ ) 0 4 π l ( ( 1 2 σ ) n + 3 m 2 n ) ) a 2 g d ω ( l )
S 1223 = 3 16 π ( 1 σ ) 0 4 π l m 2 n b 2 g d ω ( l ) + ( 1 2 σ ) 16 π ( 1 σ ) 0 4 π l n a 2 g d ω ( l ) + 3 16 π ( 1 σ ) 0 4 π m 2 n l a 2 g d ω ( l )
Now note that all the 3 integrals on the RHS of the previous equation have l , m raised to odd powers. As we previously stated, all integrals of the form 0 4 π l i m j n k d ω g shall vanish if any one of the i , j , k are odd. Therefore we can write
0 4 π l m 2 n b 2 g d ω ( l ) = 0 , 0 4 π l n a 2 g d ω ( l ) = 0 , 0 4 π m 2 n l a 2 g d ω ( l ) = 0
Therefore the S 1223 component of the Eshelby Tensor shall be
S 1223 = 3 16 π ( 1 σ ) 0 4 π l m 2 n b 2 g d ω ( l ) + ( 1 2 σ ) 16 π ( 1 σ ) 0 4 π l n a 2 g d ω ( l ) + 3 16 π ( 1 σ ) 0 4 π m 2 n l a 2 g d ω ( l ) = 0
S 1223 = 0
Therefore the S 1232 component of the Eshelby Tensor shall be
S 1232 = 1 16 π ( 1 σ ) 0 4 π ( λ 2 g 132 ( l ) + λ 1 g 232 ( l ) ) g d ω ( l ) = 1 16 π ( 1 σ ) 0 4 π ( m b 2 g 132 ( l ) + l a 2 g 232 ( l ) ) g d ω ( l )
S 1232 = 1 16 π ( 1 σ ) 0 4 π m g 132 ( l ) b 2 g d ω ( l ) + 1 16 π ( 1 σ ) 0 4 π l g 232 ( l ) a 2 g d ω ( l )
S 1232 = 3 16 π ( 1 σ ) 0 4 π l m 2 n b 2 g d ω ( l ) + 1 16 π ( 1 σ ) 0 4 π l ( ( 1 2 σ ) n + 3 m 2 n ) a 2 g d ω ( l )
S 1232 = 3 16 π ( 1 σ ) 0 4 π l m 2 n b 2 g d ω ( l ) + ( 1 2 σ ) 16 π ( 1 σ ) 0 4 π l n a 2 g d ω ( l ) + 3 16 π ( 1 σ ) 0 4 π m 2 n l a 2 g d ω ( l )
Now note that all the 3 integrals on the RHS of the previous equation have l , m raised to odd powers. As we previously stated, all integrals of the form 0 4 π l i m j n k d ω g shall vanish if any one of the i , j , k are odd. Therefore we can write
0 4 π l m 2 n b 2 g d ω ( l ) = 0 , 0 4 π l n a 2 g d ω ( l ) = 0 , 0 4 π m 2 n l a 2 g d ω ( l ) = 0
Therefore we can write
S 1232 = 3 16 π ( 1 σ ) 0 4 π l m 2 n b 2 g d ω ( l ) + ( 1 2 σ ) 16 π ( 1 σ ) 0 4 π l n a 2 g d ω ( l ) + 3 16 π ( 1 σ ) 0 4 π m 2 n l a 2 g d ω ( l ) = 0
S 1232 = 0
Let us define the terms
Q = ( 1 2 σ ) 8 π ( 1 σ ) , R = 3 8 π ( 1 σ )
I a = 0 4 π l 2 a 2 d ω ( l ) g , I a a = 0 4 π l 4 a 4 d ω ( l ) g , I a b = 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g , I a c = 0 4 π l 2 a 2 n 2 c 2 d ω ( l ) g
Therefore the S 1111 , S 1122 , S 1133 , S 1212 components of the Eshelby Tensor can be written as
S 1111 = R I a + Q a 2 I a a
S 1122 = R I a + Q b 2 I a b
S 1133 = R I a + Q c 2 I a c
S 1212 = R 2 ( I a + I b ) + Q 2 ( a 2 + b 2 ) I a b
S 1112 = 0
S 1223 = 0
S 1232 = 0
The rest of the nonzero terms can be found by cyclic permutation of the above formulas. Notice that we should also let a b c together with 1 2 3 . Using the Routh Integrals, The I a , I b , I c terms can be written in terms of standard elliptic integrals,
I a = 0 4 π l 2 a 2 d ω ( l ) g = 2 π a b c 0 d u ( a 2 + u ) Δ = 4 π a b c ( a 2 b 2 ) ( a 2 c 2 ) 1 / 2 F ( θ , k ) E ( θ , k ) , where θ = arcsin a 2 c 2 a 2 , k = a 2 b 2 a 2 c 2 .
I b = 0 4 π m 2 b 2 d ω ( l ) g = 2 π a b c 0 d u ( b 2 + u ) Δ = 4 π a b c ( a 2 b 2 ) ( b 2 c 2 ) 1 / 2 F ( θ , k ) E ( θ , k ) , where θ = arcsin a 2 c 2 b 2 , k = a 2 b 2 b 2 c 2 .
I c = 0 4 π n 2 c 2 d ω ( l ) g = 2 π a b c 0 d u ( c 2 + u ) Δ = 4 π a b c ( a 2 c 2 ) ( b 2 c 2 ) 1 / 2 F ( θ , k ) E ( θ , k ) , where θ = arcsin a 2 b 2 c 2 , k = a 2 c 2 b 2 c 2 .
Using the Routh Integrals, The I a a , I b b , I c c terms can be written as
I a a = 0 4 π l 4 a 4 d ω ( l ) g = 2 π a b c 0 d u ( a 2 + u ) 2 Δ , I b b = 0 4 π m 4 b 4 d ω ( l ) g = 2 π a b c 0 d u ( b 2 + u ) 2 Δ
I c c = 0 4 π n 4 c 4 d ω ( l ) g = 2 π a b c 0 d u ( c 2 + u ) 2 Δ
Using the Routh Integrals, The I a b , I a c , I b c terms can be written as
I a b = 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g = 2 3 π a b c 0 d u ( a 2 + u ) ( b 2 + u ) Δ , I a c = 0 4 π l 2 a 2 n 2 c 2 d ω ( l ) g = 2 3 π a b c 0 d u ( a 2 + u ) ( c 2 + u ) Δ
I b c = 0 4 π m 2 b 2 n 2 c 2 d ω ( l ) g = 2 3 π a b c 0 d u ( b 2 + u ) ( c 2 + u ) Δ
The I terms also satisfy the following properties:
I a + I b + I c = 0 4 π l 2 a 2 d ω ( l ) g + 0 4 π m 2 b 2 d ω ( l ) g + 0 4 π n 2 c 2 d ω ( l ) g = 0 4 π ( l 2 a 2 + m 2 b 2 + l 2 a 2 ) d ω ( l ) g = 0 4 π g d ω ( l ) g = 0 4 π ω ( l ) = 4 π
I a + I b + I c = 4 π
a 2 I a a + b 2 I a b + c 2 I a c = a 2 0 4 π l 4 a 4 d ω ( l ) g + b 2 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g + c 2 0 4 π l 2 a 2 n 2 c 2 d ω ( l ) g = 0 4 π l 2 a 2 l 2 d ω ( l ) g + 0 4 π l 2 a 2 m 2 d ω ( l ) g + 0 4 π l 2 a 2 n 2 d ω ( l ) g
a 2 I a a + b 2 I a b + c 2 I a c = 0 4 π l 2 a 2 ( l 2 + m 2 + n 2 ) d ω ( l ) g = 0 4 π l 2 a 2 d ω ( l ) g = I a
a 2 I a a + b 2 I a b + c 2 I a c = I a
I a a + I a b + I a c = 0 4 π l 4 a 4 d ω ( l ) g + 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g + 0 4 π l 2 a 2 n 2 c 2 d ω ( l ) g = 0 4 π l 2 a 2 ( l 2 a 2 + m 2 b 2 + n 2 c 2 ) d ω ( l ) g = 0 4 π l 2 a 2 ( g ) d ω ( l ) g
I a a + I a b + I a c = 0 4 π l 2 a 2 d ω ( l ) = 4 π 3 a 2
Proof of 0 4 π l 2 a 2 d ω ( l ) = 4 π 3 a 2 :
For direction cosines, the integral over the solid angle can be expressed in spherical coordinates where l = cos θ , with θ being the polar angle. The differential solid angle element d ω in spherical coordinates is given by: d ω = sin θ d θ d ϕ where θ varies from 0 to π and ϕ varies from 0 to 2 π . Thus, the integral can be rewritten as:
0 4 π l 2 a 2 d ω = 1 a 2 0 2 π 0 π cos 2 θ sin θ d θ d ϕ
Let’s evaluate the integral step by step. First, integrate with respect to θ :
0 π cos 2 θ sin θ d θ
Let u = cos θ , so d u = sin θ d θ . The limits of integration change from θ = 0 to θ = π to u = 1 to u = 1 . The integral becomes: 1 1 u 2 d u = 1 1 u 2 d u . This integral evaluates to: u 3 3 1 1 = 1 3 1 3 = 2 3 . Now, integrate with respect to ϕ : 0 2 π d ϕ = 2 π Multiply the results of the two integrals: 0 4 π l 2 a 2 d ω = 1 a 2 × 2 3 × 2 π = 4 π 3 a 2 . So, the result of the integral is:
0 4 π l 2 a 2 d ω = 4 π 3 a 2
Therefore we have
I 11 + I 12 + I 13 = 0 4 π l 2 a 2 d ω ( l ) = 4 π 3 a 2
2 π a b c 0 d u ( a 2 + u ) 2 Δ + 2 3 π a b c 0 d u ( a 2 + u ) ( b 2 + u ) Δ
+ 2 3 π a b c 0 d u ( a 2 + u ) ( c 2 + u ) Δ = π a b c 0 ( 6 ( b 2 + u ) ( c 2 + u ) + 2 ( a 2 + u ) ( c 2 + u ) + 2 ( a 2 + u ) ( b 2 + u ) ) 3 ( a 2 + u ) 2 Δ
3 a 2 I 11 + b 2 I 12 + c 2 I 13 = 3 I 1
I 12 = I 2 I 1 a 2 b 2
and the standard elliptic integrals are defined as
F ( θ , k ) = 0 θ d w ( 1 k 2 sin 2 w ) 1 / 2
E ( θ , k ) = 0 θ ( 1 k 2 sin 2 w ) 1 / 2 d w

9.3. Elliptic Cylinder Inclusion

For an elliptic cylinder inclusion ( c ). The components of Eshelby tensor are
S 1111 = 1 2 ( 1 ν ) b 2 + 2 a b ( a + b ) 2 + ( 1 2 ν ) b a + b
S 2222 = 1 2 ( 1 ν ) a 2 + 2 a b ( a + b ) 2 + ( 1 2 ν ) a a + b
S 3333 = 0
S 1122 = 1 2 ( 1 ν ) b 2 ( a + b ) 2 ( 1 2 ν ) b a + b
S 2233 = 1 2 ( 1 ν ) 2 ν a a + b
S 2211 = 1 2 ( 1 ν ) a 2 ( a + b ) 2 ( 1 2 ν ) a a + b
S 3311 = S 3322 = 0
S 1212 = 1 2 ( 1 ν ) a 2 + b 2 2 ( a + b ) 2 + ( 1 2 ν ) 2
S 1133 = 1 2 ( 1 ν ) 2 ν b a + b
S 2323 = a 2 ( a + b )
S 3131 = b 2 ( a + b )

9.4. Flat Ellipsoid Inclusion

For a flat ellipsoid ( a > b c ). The I integrals in this limiting case reduce to
I 1 = 4 π ( F ( k ) E ( k ) ) b c a 2 b 2
I 2 = 4 π E ( k ) c ( F ( k ) E ( k ) ) b c a 2 b 2
I 3 = 4 π 1 E ( k ) c
I 12 = 4 π a 2 b 2 E ( k ) c 2 ( F ( k ) E ( k ) ) b c a 2 b 2
I 23 = 4 π b 2 1 2 E ( k ) c + ( F ( k ) E ( k ) ) b c a 2 b 2
I 31 = 4 π a 2 1 E ( k ) c ( F ( k ) E ( k ) ) b c a 2 b 2
I 33 = 4 π 3 c 2
where E ( k ) and F ( k ) are complete elliptic integrals defined as
F ( k ) = 0 π 2 d w 1 k 2 sin 2 w
E ( k ) = 0 π 2 1 k 2 sin 2 w d w

9.5. Penny Shaped Inclusion

We have a penny shaped inclusion when a = b . We can compute the S 1111 , S 1122 , S 1133 , S 1212 , S 1112 , S 1223 , S 1232 components of the Eshelby Tensor by computing the values of I a , I b , I a a , I a b , I a c which were defined in the previous subsection.
Using Routh Integrals we can reduce the solid angle integral 0 4 π l 2 i m 2 j n 2 k d ω g to simple integrals.
I a = 0 4 π l 2 a 2 d ω ( l ) g = 2 π a b c 0 d u ( a 2 + u ) Δ , I a a = 0 4 π l 4 a 4 d ω ( l ) g = 2 π a b c 0 d u ( a 2 + u ) 2 Δ
I a b = 0 4 π l 2 a 2 m 2 b 2 d ω ( l ) g = 2 3 π a b c 0 d u ( a 2 + u ) ( b 2 + u ) Δ , I a c = 0 4 π l 2 a 2 n 2 c 2 d ω ( l ) g = 2 3 π a b c 0 d u ( a 2 + u ) ( c 2 + u ) Δ
where Δ is given by the expression
Δ = ( a 2 + u ) 1 2 ( b 2 + u ) 1 2 ( c 2 + u ) 1 2 = ( a 2 + u ) ( c 2 + u ) 1 2
Δ = ( a 2 + u ) ( c 2 + u ) 1 2
The S 1111 , S 1122 , S 1133 , S 1212 , S 1112 , S 1223 , S 1232 components of the Eshelby Tensor for Penny Shaped Inclusion are
S 1111 = S 2222 = π ( 13 8 ν ) 32 ( 1 ν ) c a
S 3333 = 1 π ( 1 2 ν ) 4 ( 1 ν ) c a
S 1122 = S 2211 = π ( 8 ν 1 ) 32 ( 1 ν ) c a
S 1133 = S 2233 = π ( 2 ν 1 ) 8 ( 1 ν ) c a
S 3311 = S 3322 = v 1 v 1 π ( 4 ν + 1 ) 8 ν c a
S 1212 = π ( 7 8 ν ) 32 ( 1 ν ) c a
S 3131 = S 2323 = 1 2 1 + π ( ν 2 ) 4 ( 1 ν ) c a

10. Eshelby’s Tensor for 2-Dimensional Inclusions

This derivation focuses on Eshelby’s tensor in the context of 2-dimensional (2D) elliptical inclusions in an isotropic elastic medium under plane strain conditions. Consider a 2D elliptical inclusion in an infinite, homogeneous, isotropic elastic medium. The inclusion is characterized by an eigenstrain e * , which is the strain the inclusion would undergo if isolated from the surrounding matrix. The 2D elliptical inclusion is defined by:
x 1 2 a 2 + x 2 2 b 2 1
where a and b are the semi-axes of the ellipse along the x 1 and x 2 directions, respectively.
For an isotropic elastic medium, the stress-strain relation (Hooke’s law) in 2D under plane strain conditions is:
σ i j = C i j k l e k l
where σ i j is the stress tensor, e k l is the strain tensor, and C i j k l is the elastic constant tensor given by:
C i j k l = λ δ i j δ k l + μ ( δ i k δ j l + δ i l δ j k )
Here, λ and μ are the Lamé constants. The inclusion undergoes an eigenstrain e i j * , leading to the following stress and strain fields inside the inclusion:
σ i j = C i j k l ( e k l e k l * )
The goal is to determine the Eshelby tensor S i j k l , which relates the elastic strain e i j inside the inclusion to the eigenstrain e k l * :
e i j = S i j k l e k l *
The Green’s function G i j ( x ) represents the displacement at point x due to a unit force applied at x . For an isotropic, infinite elastic medium in 2D, the Green’s function in Fourier space is:
g i j ( k ) = ( z z ) i j 1 k 2
where z = k | k | . The Green’s function in real space is then:
G i j ( x ) = 1 4 π 2 exp ( i k · x ) k 2 ( z z ) i j 1 d k
The auxiliary tensor D i j k l relates the gradient of the Green’s function to the eigenstress inside the inclusion:
D i j k l ( x ) = V 0 G i j , k l ( x x ) d V ( x )
Using the Fourier representation of the Green’s function:
D i j k l ( x ) = 1 ( 2 π ) 2 exp ( i k · x ) ( z z ) i j 1 z k z l k 2 Q ( k ) d k
where the integral Q ( k ) is defined as:
Q ( k ) = V 0 exp ( i k · x ) d V ( x )
where V 0 is the volume (in this case, area since we are in 2D) of the inclusion. For an elliptical inclusion, the domain V 0 is given by:
x 1 2 a 2 + x 2 2 b 2 1
Here, x 1 and x 2 are the coordinates within the inclusion, and a, b are the semi-major and semi-minor axes of the ellipse, respectively. We perform a change of variables to transform the elliptical domain into a circular domain, making the integral easier to handle. Define new variables:
X = x 1 a , Y = x 2 b
In these new variables, the elliptical domain becomes:
X 2 + Y 2 1
which represents a unit circle in the X Y -plane. The differential volume element transforms as:
d V ( x ) = a · b d X d Y
Thus, the integral Q ( k ) can be rewritten in terms of X and Y :
Q ( k ) = a b X 2 + Y 2 1 exp i k 1 · ( a X ) + i k 2 · ( b Y ) d X d Y
Here, k 1 and k 2 are the components of the wavevector k in the x 1 and x 2 directions, respectively. Expressing k · x in polar coordinates:
k · x = k 1 x 1 + k 2 x 2 = a k 1 X + b k 2 Y
Let us define
λ x = a k 1 , λ y = b k 2
Then we have
k · x = λ x X + λ y Y
We now express X and Y in polar coordinates ( R , θ ) :
X = R cos θ , Y = R sin θ
Therefore we can write
k · x = R ( λ x cos θ + λ y sin θ ) = R λ cos ϕ
where ϕ is the angle between λ = ( λ x , λ y ) and ( X , Y ) . Integrating Over the Elliptical Domain. The integral Q ( k ) now becomes:
Q ( k ) = a b 0 1 0 2 π exp ( i R λ cos ϕ ) R d θ d R
The integral over θ can be recognized as a standard Bessel function of the first kind J 0 :
0 2 π exp ( i R λ cos ϕ ) d θ = 2 π J 0 ( R λ )
Thus, Q ( k ) simplifies to:
Q ( k ) = 2 π a b 0 1 R J 0 ( R λ ) d R
The integral over R can be evaluated using the known result for the integral involving the Bessel function:
0 1 R J 0 ( R λ ) d R = 1 λ 2 ( sin λ λ cos λ )
Substituting back, we get:
Q ( k ) = 2 π a b · 1 λ 2 ( sin λ λ cos λ )
Finally, the expression for Q ( k ) is:
Q ( k ) = 2 π a b λ 3 ( sin λ λ cos λ )
Here, λ = ( a k 1 ) 2 + ( b k 2 ) 2 .
Let’s now do the Eshelby Tensor Derivation. Substituting Q ( k ) back into the expression for D i j k l and using polar coordinates to evaluate the integrals:
D i j k l ( x ) = a b c 2 π 2 0 2 π ( z z ) i j 1 z k z l λ 3 sin λ λ cos λ d θ
For a circular inclusion (where a = b ):
D i j k l = 1 2 π 0 2 π ( z z ) i j 1 z k z l μ d θ
The Eshelby tensor S i j k l for Circular Inclusion is obtained using:
S i j m n = 1 2 C k l m n D i k l j + D j k l i
Substituting the expression for D i j k l :
S i j m n = ν 2 ( 1 ν ) δ i j δ m n + 1 16 ( 1 ν ) ( 6 8 ν ) ( δ i n δ j m + δ j n δ i m ) 2 δ i j δ m n
For a circular inclusion, this simplifies to:
S i j m n = 4 ν 1 8 ( 1 ν ) δ i j δ m n + 3 4 ν 8 ( 1 ν ) ( δ i n δ j m + δ j n δ i m )
The Eshelby tensor S i j m n for a 2D inclusion (circular or elliptical) in an isotropic elastic medium has been derived rigorously. This tensor allows us to determine the strain field inside the inclusion given the eigenstrain e k l * . The derivation used the Green’s function approach, combined with Fourier transforms, change of variables, and integration techniques to arrive at the final expressions.

11. Derivation of the Green’s Function for Anisotropic Medium

To begin, we derive the Green’s function G i j ( x x ) for an infinite, homogeneous, and anisotropic elastic medium. This Green’s function represents the displacement in the i-th direction at point x due to a unit point force applied in the j-th direction at point x . We start with the Equilibrium Equation for an Infinite Body. The Green’s function satisfies the following equilibrium equation:
C i j k l 2 G k m ( x ) x j x l + δ i m δ ( x ) = 0
where C i j k l is the fourth-order stiffness (elasticity) tensor, and δ ( x ) is the Dirac delta function representing the point force at the origin. To solve this equation, we take the Fourier transform. The Fourier transform f ^ ( k ) of a function f ( x ) is defined as:
f ^ ( k ) = R 3 f ( x ) e i k · x d x ,
with the inverse Fourier transform given by:
f ( x ) = 1 ( 2 π ) 3 R 3 f ^ ( k ) e i k · x d k .
Applying the Fourier transform to the equilibrium equation, we have:
C i j k l k j k l g k m ( k ) = δ i m
where g i j ( k ) is the Fourier transform of G i j ( x ) . The equation simplifies to:
( z z ) i k g k m ( k ) = δ i m
where ( z z ) i k = C i j k l k j k l . Taking the inverse of ( z z ) i k , we find:
g i j ( k ) = ( z z ) i j 1 .
Thus, the Green’s function in Fourier space is:
g i j ( k ) = ( z z ) i j 1 k 2 .
The real-space Green’s function G i j ( x ) is obtained by taking the inverse Fourier transform:
G i j ( x ) = 1 ( 2 π ) 3 R 3 g i j ( k ) e i k · x d k .
Substituting the expression for g i j ( k ) :
G i j ( x ) = 1 ( 2 π ) 3 R 3 ( z z ) i j 1 k 2 e i k · x d k
To evaluate this integral, we switch to spherical coordinates in Fourier space:
k x = k sin θ cos ϕ , k y = k sin θ sin ϕ , k z = k cos θ ,
with the volume element given by k 2 sin θ d k d θ d ϕ . The integral becomes:
G i j ( x ) = 1 ( 2 π ) 2 0 sin ( k R ) k R ( z z ) i j 1 d k
where R = | x | is the distance from the source point. After evaluating the integral, the Green’s function in real space is:
G i j ( x ) = 1 8 π R ( z z ) i j 1
For isotropic materials, this can be further simplified using the Lame constants λ and μ :
G i j ( x ) = 1 8 π μ R λ + 3 μ λ + 2 μ δ i j + λ + μ λ + 2 μ x i x j R 2
This Green’s function represents the displacement field in an infinite, homogeneous, and isotropic elastic medium due to a point force.

12. Derivation of the Auxiliary Tensor D i j k l

Using the derived Green’s function, we now derive the auxiliary tensor D i j k l that relates the constrained displacement gradients to the eigenstress inside an inclusion. The auxiliary tensor D i j k l is defined by:
u i , l c ( x ) = σ k j * D i j k l ( x )
where u i , l c ( x ) is the constrained displacement gradient, and σ k j * is the eigenstress. The tensor D i j k l is related to Eshelby’s tensor S i j k l :
S i j m n e m n * = e i j c ,
where the constrained strain e i j c is given by:
e i j c = 1 2 u i , j c + u j , i c .
The constrained displacement gradient can also be expressed using the Green’s function as:
u i , j c ( x ) = S 0 σ l k * n k G i l , j ( x x ) d S ( x )
By comparing this with the earlier definition, we identify:
D i l k j ( x ) = S 0 n k G i l , j ( x x ) d S ( x ) ,
or equivalently:
D i j k l ( x ) = S 0 G i j , l ( x x ) n k d S ( x )
Using Gauss’s theorem, this surface integral is converted into a volume integral:
D i j k l ( x ) = V 0 x k G i j , l ( x x ) d V ( x ) .
This simplifies to:
D i j k l ( x ) = V 0 G i j , k l ( x x ) d V ( x )
where G i j , k l ( x x ) is the second derivative of the Green’s function. Substituting the Fourier-transformed Green’s function into the expression for D i j k l ( x ) , we get:
D i j k l ( x ) = V 0 2 x k x l 1 ( 2 π ) 3 exp i k · ( x x ) ( z z ) i j 1 k 2 d k d V ( x ) .
Taking the second derivatives with respect to x k and x l , we obtain:
D i j k l ( x ) = 1 ( 2 π ) 3 V 0 i k k i k l exp i k · ( x x ) ( z z ) i j 1 k 2 d k d V ( x ) .
Simplifying, this becomes:
D i j k l ( x ) = 1 ( 2 π ) 3 V 0 exp i k · ( x x ) ( z z ) i j 1 k k k l k 2 d k d V ( x )
Recognizing that k k k l = | k | 2 z k z l , where z k and z l are the components of the unit vector in the direction of k, we rewrite the equation as:
D i j k l ( x ) = 1 ( 2 π ) 3 V 0 exp i k · ( x x ) ( z z ) i j 1 z k z l k 2 d k d V ( x )
The function Q ( k ) is defined as the Fourier transform of the characteristic function of the inclusion V 0 :
Q ( k ) = V 0 exp i k · x d V ( x )
For an ellipsoidal inclusion, the region V 0 is defined by:
x 2 a 2 + y 2 b 2 + z 2 c 2 1 ,
where a, b, and c are the semi-axes of the ellipsoid. The Fourier transform Q ( k ) can be computed as follows. Let:
x = a X , y = b Y , z = c Z ,
where X, Y, and Z are dimensionless variables that satisfy X 2 + Y 2 + Z 2 1 . Thus, the region V 0 can be expressed as:
V 0 = ( X , Y , Z ) R 3 : X 2 + Y 2 + Z 2 1
The volume element transforms as:
d V ( x ) = a b c d X d Y d Z .
The Fourier transform Q ( k ) is then given by:
Q ( k ) = a b c V 0 exp i a k x X + b k y Y + c k z Z d X d Y d Z
This integral can be separated into three independent integrals over X, Y, and Z:
Q ( k ) = a b c 1 1 exp ( i a k x X ) d X 1 1 exp ( i b k y Y ) d Y 1 1 exp ( i c k z Z ) d Z .
Each integral is of the form:
1 1 exp ( i λ ξ ) d ξ = 2 sin ( λ ) λ ,
where λ = a k x , b k y , or c k z . Thus:
Q ( k ) = a b c 2 sin ( a k x ) a k x 2 sin ( b k y ) b k y 2 sin ( c k z ) c k z .
Alternatively, the Fourier transform can be expressed in spherical coordinates:
Q ( k ) = a b c λ 3 4 π sin λ λ cos λ
where λ = a 2 k x 2 + b 2 k y 2 + c 2 k z 2 . Substituting the expression for Q ( k ) into the formula for D i j k l ( x ) , we get:
D i j k l ( x ) = a b c 2 π 2 ( z z ) i j 1 z k z l λ 3 sin λ λ cos λ exp i k · x d k .
To simplify the integral, we switch to spherical coordinates in Fourier space:
k x = k sin Φ cos Θ , k y = k sin Φ sin Θ , k z = k cos Φ ,
where λ becomes:
λ = k a 2 cos 2 Θ + b 2 sin 2 Θ sin Φ + c 2 cos 2 Φ
The integral expression for D i j k l ( x ) then becomes:
D i j k l ( x ) = a b c 2 π 2 0 π 0 2 π ( z z ) i j 1 z k z l β 3 / 2 sin λ λ cos λ e i k · x k 2 sin Φ d k d Θ d Φ
where β = a 2 cos 2 Θ + b 2 sin 2 Θ sin 2 Φ + c 2 cos 2 Φ . To further simplify the expression for D i j k l , we integrate over the angular coordinates Θ and Φ . This yields:
D i j k l ( x ) = a b c ( 2 π ) 2 0 0 π 0 2 π ( z z ) i j 1 z k z l sin Φ β 3 / 2 d Θ d Φ k 2 sin ( k λ ) k λ cos ( k λ ) e i k · x d k .
This integral simplifies further when evaluated within the ellipsoid due to its symmetry, leading to:
D i j k l = a b c 4 π 0 π 0 2 π ( z z ) i j 1 z k z l β 3 / 2 sin Φ d Θ d Φ
After evaluating the integral, D i j k l inside the ellipsoid becomes a constant tensor, and the final expression is given by:
D i j k l = 1 ( 2 π ) 3 exp i k · x ( z z ) i j 1 z k z l Q ( k ) d k ,
Thus, the auxiliary tensor D i j k l for an ellipsoidal inclusion is a constant value inside the inclusion, given by:
D i j k l = a b c 4 π 0 π 0 2 π ( z z ) i j 1 z k z l β 3 / 2 sin Φ d Θ d Φ
This derivation shows that the auxiliary tensor becomes a constant within an ellipsoidal inclusion, providing insights into the behavior of materials with embedded inclusions and forming the basis for understanding internal stress and strain fields within such materials.

13. Derivation of Inclusion Energy in an Infinite Solid

We start by defining the elastic fields inside the inclusion (denoted by superscript I) and outside the inclusion in the matrix (denoted by superscript M). These fields include stress ( σ i j ), strain ( e i j ), and displacement ( u i ). For a homogeneous infinite solid, these fields are given by:
For Matrix : e i j M = e i j c , σ i j M = σ i j c , u i M = u i c
For Inclusion : e i j I = e i j c e i j * , σ i j I = σ i j c σ i j * , u i I = u i c e i j * x j
The total elastic energy E stored in the solid due to the inclusion is given by the sum of the energies inside and outside the inclusion:
E = 1 2 V 0 σ i j I e i j I d V + 1 2 V V 0 σ i j M e i j M d V
Rewriting E in terms of displacements:
E = 1 4 V 0 σ i j I ( u i , j I + u j , i I ) d V + 1 4 V V 0 σ i j M ( u i , j M + u j , i M ) d V
Due to the symmetry of the stress tensor:
E = 1 2 V 0 σ i j I u j , i I d V + 1 2 V V 0 σ i j M u j , i M d V
Using the identity:
σ i j u i , j = ( σ i j u j ) , i σ i j , i u j
And assuming there are no body forces (so σ i j , i = 0 ):
E = 1 2 V 0 ( σ i j I u j I ) , i d V + 1 2 V V 0 ( σ i j M u j M ) , i d V
Using Gauss’s theorem, these volume integrals are converted to surface integrals:
E = 1 2 S 0 σ i j I u j I n i out d S 1 2 S 0 σ i j M u j M n i out d S + 1 2 S σ i j M u j M n i d S
As S tends to infinity, the surface integral over S vanishes:
E = 1 2 S 0 σ i j I u j I σ i j M u j M n i out d S
Since the traction across the interface S 0 must be continuous ( σ i j I n i out = σ i j M n i out ):
E = 1 2 S 0 σ i j I ( u j I u j M ) n i out d S
Given that u j I u j M = e j k * x k :
E = 1 2 S 0 σ i j I n i out e j k * x k d S
Transforming the surface integral back into a volume integral:
E = 1 2 V 0 e j k * σ i j I d V = 1 2 e i j * V 0 σ i j I d V
For an ellipsoidal inclusion, where the stress inside is uniform:
E = 1 2 σ i j I e i j * V 0
Thus, the total inclusion energy in an infinite solid is:
E = 1 2 σ i j I e i j * V 0
We know that the stress σ i j I in the inclusion shall be
σ i j I = σ i j C σ i j T = λ ( e C e T ) δ i j + 2 μ ( e i j C e i j T )
As we derived earlier, the strain e i j C in the material is given by
e i j C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i j ) ( e i k T 4 π ϕ , k j + e j k T 4 π ϕ , k i e l k T 8 π ( 1 σ ) ψ , i j k l ) )
As we derived earlier, the trace of the strain tensor e i j C shall be
e i i C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i i ) ( e i k T 4 π ϕ , k i + e i k T 4 π ϕ , k i e l k T 8 π ( 1 σ ) ψ , i i k l ) )
e i i C = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i i ) ( e i k T 2 π ϕ , k i e l k T 8 π ( 1 σ ) ψ , i i k l ) )
We have earlier derived that ϕ , i i = 4 π inside the inclusion. We also derived that ψ , i i k l = 2 ϕ , k l . Therefore the above boxed equation can be written as
e i i C = ( λ ( 1 2 σ ) e T 2 μ ( 1 σ ) ) ( e i k T 2 π ϕ , k i e l k T 4 π ( 1 σ ) ϕ , k l ) )
Repeated indices k , l in the third term of the previous equation can be written as k , i .
e i i C = ( λ ( 1 2 σ ) e T 2 μ ( 1 σ ) ) ( e i k T 2 π ϕ , k i e i k T 4 π ( 1 σ ) ϕ , k i ) )
e i i C = ( λ ( 1 2 σ ) 2 μ ( 1 σ ) ) e T e i k T ϕ , k i ( 1 2 π 1 4 π ( 1 σ ) )
e i i C = ( λ ( 1 2 σ ) 2 μ ( 1 σ ) ) e T e i k T ϕ , k i ( ( 2 ( 1 σ ) 1 ) 4 π ( 1 σ ) )
e i i C = ( λ ( 1 2 σ ) 2 μ ( 1 σ ) ) e T e i k T ϕ , k i ( 1 2 σ 4 π ( 1 σ ) )
Note that we have
λ ( 1 2 σ ) = 2 μ σ
λ ( 1 2 σ ) 2 μ ( 1 σ ) = σ ( 1 σ )
Therefore the trace of the strain tensor e i j C shall be
e C = e i i C = σ ( 1 σ ) e T e l k T ϕ , k l ( 1 2 σ 4 π ( 1 σ ) )
Therefore we have
e C e T = ( σ ( 1 σ ) 1 ) e T e l k T ϕ , k l ( 1 2 σ 4 π ( 1 σ ) )
e C e T = ( 2 σ 1 ) ( 1 σ ) e T e l k T ϕ , k l ( 1 2 σ 4 π ( 1 σ ) )
( e C e T ) δ i j = ( 2 σ 1 ) ( 1 σ ) e T δ i j e l k T ϕ , k l δ i j ( 1 2 σ 4 π ( 1 σ ) )
For uniform expansion, we have e l k T = 1 3 e T δ l k . Therefore we can write the above boxed equation as
( e C e T ) δ i j = ( 2 σ 1 ) ( 1 σ ) e T δ i j 1 3 e T δ l k ϕ , k l δ i j ( 1 2 σ 4 π ( 1 σ ) )
( e C e T ) δ i j = ( 2 σ 1 ) ( 1 σ ) e T δ i j 1 3 e T ϕ , k k δ i j ( 1 2 σ 4 π ( 1 σ ) )
Inside the inclusion, we have ϕ , k k = 4 π . Therefore we can write the above equation as
( e C e T ) δ i j = ( 2 σ 1 ) ( 1 σ ) e T δ i j + e T δ i j ( 1 2 σ 3 ( 1 σ ) )
( e C e T ) δ i j = ( ( 2 σ 1 ) ( 1 σ ) + ( 1 2 σ ) 3 ( 1 σ ) ) e T δ i j
( e C e T ) δ i j = ( ( 6 σ 3 ) + ( 1 2 σ ) 3 ( 1 σ ) ) e T δ i j = ( 4 σ 2 ) 3 ( 1 σ ) e T δ i j = 2 ( 2 σ 1 ) 3 ( 1 σ ) e T δ i j
( e C e T ) δ i j = 2 ( 2 σ 1 ) 3 ( 1 σ ) e T δ i j
Similarly, we also have
e i j C e i j T = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i j ) ( e i k T 4 π ϕ , k j + e j k T 4 π ϕ , k i e l k T 8 π ( 1 σ ) ψ , i j k l ) e i j T
For uniform expansion, we have e l k T = 1 3 e T δ l k . Therefore we can write the above equation as
e i j C e i j T = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i j ) ( 1 12 π δ i k ϕ , k j e T + 1 12 π δ j k ϕ , k i e T 1 24 π ( 1 σ ) ψ , i j k l δ l k e T ) 1 3 e T δ i j
e i j C e i j T = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i j ) ( 1 12 π ϕ , i j e T + 1 12 π ϕ , j i e T 1 24 π ( 1 σ ) ψ , i j k k e T ) 1 3 e T δ i j
Now note that we have ϕ , i j = ϕ , j i and ψ , i j k k = 2 ϕ i j , therefore we can write the above equation as
e i j C e i j T = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i j ) ( 1 6 π ϕ , i j e T 1 12 π ( 1 σ ) ϕ , i j e T ) 1 3 e T δ i j
e i j C e i j T = ( λ ( 2 σ 1 ) e T 8 π μ ( 1 σ ) ϕ , i j ) ( 1 6 π 1 12 π ( 1 σ ) ) ϕ , i j e T 1 3 e T δ i j
e i j C e i j T = ( λ ( 2 σ 1 ) 8 π μ ( 1 σ ) e T ϕ , i j ) ( ( 2 ( 1 σ ) 1 ) 12 π ( 1 σ ) ) ϕ , i j e T 1 3 e T δ i j
e i j C e i j T = ( λ ( 2 σ 1 ) 8 π μ ( 1 σ ) e T ϕ , i j ) ( 1 2 σ 12 π ( 1 σ ) ) ϕ , i j e T 1 3 e T δ i j
Note that we have
λ ( 2 σ 1 ) = 2 μ σ
λ ( 2 σ 1 ) 8 π μ ( 1 σ ) = σ 4 π ( 1 σ )
Therefore the above boxed equation can be written as
e i j C e i j T = ( σ 4 π ( 1 σ ) ) e T ϕ , i j ( 1 2 σ 12 π ( 1 σ ) ) e T ϕ , i j 1 3 e T δ i j
e i j C e i j T = ( 3 σ + 1 2 σ 12 π ( 1 σ ) ) e T ϕ , i j 1 3 e T δ i j = ( 1 + σ 12 π ( 1 σ ) ) e T ϕ , i j 1 3 e T δ i j
e i j C e i j T = ( 1 + σ 12 π ( 1 σ ) ) e T ϕ , i j 1 3 e T δ i j
Therefore the stress at the inclusion σ i j I shall be
σ i j I = σ i j C σ i j T = λ ( e C e T ) δ i j + 2 μ ( e i j C e i j T )
σ i j I = λ ( 2 ( 2 σ 1 ) 3 ( 1 σ ) e T δ i j ) δ i j + 2 μ ( ( 1 + σ 12 π ( 1 σ ) ) e T ϕ , i j 1 3 e T δ i j )
σ i j I = 2 λ ( 2 σ 1 ) 3 ( 1 σ ) e T δ i j μ ( 1 + σ ) 6 π ( 1 σ ) e T ϕ , i j 2 μ 3 e T δ i j
Note that we have
λ ( 2 σ 1 ) = 2 μ σ
2 λ ( 2 σ 1 ) 3 ( 1 σ ) = 4 μ σ 3 ( 1 σ )
Therefore we can write the stress at the inclusion σ i j I as
σ i j I = 4 μ σ 3 ( 1 σ ) e T δ i j μ ( 1 + σ ) 6 π ( 1 σ ) e T ϕ , i j 2 μ 3 e T δ i j
σ i j I = ( 4 σ 3 ( 1 σ ) + 2 3 ) μ e T δ i j μ ( 1 + σ ) 6 π ( 1 σ ) e T ϕ , i j
σ i j I = ( 4 σ + 2 ( 1 σ ) 3 ( 1 σ ) ) μ e T δ i j μ ( 1 + σ ) 6 π ( 1 σ ) e T ϕ , i j
σ i j I = 2 ( 1 + σ ) 3 ( 1 σ ) μ e T δ i j μ ( 1 + σ ) 6 π ( 1 σ ) e T ϕ , i j
We derived earlier that the total inclusion energy in an infinite solid is:
E = 1 2 σ i j I e i j T V 0
For uniform expansion, we have e i j T = 1 3 e T δ i j . Therefore the total inclusion energy in an infinite solid shall be
E = 1 2 σ i j I ( 1 3 e T δ i j ) V 0 = 1 6 σ i j I δ i j e T V 0 = 1 6 σ i i I e T V 0
Now note that
σ i i I = 2 ( 1 + σ ) 3 ( 1 σ ) μ e T δ i i μ ( 1 + σ ) 6 π ( 1 σ ) e T ϕ , i i
Now we know that δ i i = 3 and ϕ i i = 4 π . Therefore σ i i I shall be
σ i i I = 2 ( 1 + σ ) ( 1 σ ) μ e T + 2 ( 1 + σ ) 3 ( 1 σ ) μ e T
σ i i I = ( 2 ( 1 + σ ) ( 1 σ ) 2 ( 1 + σ ) 3 ( 1 σ ) ) μ e T
σ i i I = ( 6 ( 1 + σ ) 2 ( 1 + σ ) 3 ( 1 σ ) ) μ e T
σ i i I = 4 ( 1 + σ ) 3 ( 1 σ ) μ e T
Therefore the total inclusion energy in an infinite solid shall be
E = 1 6 σ i i I e T V 0 = 1 6 ( 4 ( 1 + σ ) 3 ( 1 σ ) μ e T ) e T V 0 = 2 ( 1 + σ ) 9 ( 1 σ ) μ ( e T ) 2 V 0
E = 2 ( 1 + σ ) 9 ( 1 σ ) μ ( e T ) 2 V 0

14. Derivation of Inclusion Energy in a Finite Solid

We begin by considering an inclusion in a finite solid. The stress-strain fields in this case can be solved by superposition. Suppose the finite solid assumes the stress-strain fields of an infinite solid containing an inclusion. Then, to maintain equilibrium, a set of traction forces T ˜ j must be applied to the outer surface S e x t of the solid.
To obtain the solution for a finite solid with zero traction on its outer surface, we need to remove T ˜ j on S e x t . This is equivalent to applying a canceling traction force F ˜ j = T ˜ j on S e x t . The resulting elastic fields due to this canceling traction force are called image fields, denoted by the strain e i j i m , stress σ i j i m ¯ , and displacement u i i m ¯ fields. Thus, the elastic fields inside the matrix and the inclusion are given by:
e i j M = e i j c + e i j i m , e i j I = e i j c e i j * + e i j i m , σ i j M = σ i j c + σ i j i m ¯ , σ i j I = σ i j c σ i j * + σ i j i m ¯ , u i M = u i c + u i i m ¯ , u i I = u i c e i j * x j + u i i m ¯ .
The image fields satisfy the following conditions:
e i j i m ( x ) = 1 2 u i , j i m ( x ) + u j , i i m ( x )
σ i j i m ¯ ( x ) = C i j k l e k l i m ( x )
where C i j k l is the stiffness tensor of the material. Similar to the infinite solid case, the total elastic energy E in the solid can be expressed in terms of surface integrals:
E = 1 2 S 0 σ i j I u j I σ i j M u j M n i o u t d S + S e x t σ i j M u j M n i e x t d S .
Using the free traction boundary condition on the outer surface S e x t , σ i j M n i e x t = 0 , the second integral vanishes, leaving:
E = 1 2 S 0 σ i j I u j I u j M n i o u t d S
Substituting u j I u j M = e j k * x k , we obtain:
E = 1 2 S 0 σ i j I e j k * x k n i o u t d S .
This is the same as the energy expression for an infinite solid except that the stress field inside the inclusion now contains the image component. Let σ i j I , denote the stress field inside the inclusion in an infinite medium. Then the stress inside the inclusion in the finite solid is:
σ i j I ( x ) = σ i j I , + σ i j i m ¯ ( x )
The inclusion energy in an infinite solid is given by:
E = 1 2 σ i j I , e i j * V 0
where V 0 is the volume of the inclusion. In the case of a finite solid, the inclusion energy E becomes:
E = E 1 2 V 0 σ i j i m ¯ e i j * d V .
Converting the surface integral into a volume integral and averaging the image stress inside the inclusion:
E = E 1 2 σ i j i m ¯ e i j * V 0 = E + E i m
where E i m is the image contribution to the total inclusion energy, defined as:
E i m = 1 2 σ i j i m ¯ e i j * V 0
and σ i j i m ¯ is the averaged image stress inside the inclusion:
σ i j i m ¯ = 1 V 0 V 0 σ i j i m ( x ) d V ( x ) .
The total inclusion energy in a finite solid can thus be expressed as:
E = 1 2 σ i j I e i j * V 0
where σ i j I = σ i j I , + σ i j i m ¯ .

15. Derivation of Inclusion Energy of Finite solid with Applied Tractions

Before deriving the Inclusion Energy of Finite solid with applied tractions, we first need to prove the Colonetti’s Theorem. Let us consider a solid with volume V and outer surface S. We define two stress states:
  • State 1: Purely internal, generated by an eigenstrain (or some inhomogeneity) inside the solid.
  • State 2: Purely applied, generated by external tractions on the surface of the solid without any internal eigenstrain.
The total elastic energy in each of these states can be expressed as:
E ( 1 ) = 1 2 V σ i j ( 1 ) ϵ i j ( 1 ) d V
E ( 2 ) = 1 2 V σ i j ( 2 ) ϵ i j ( 2 ) d V
where σ i j and ϵ i j denote the stress and strain tensors, respectively. When both states are present, the combined stress and strain fields can be written as:
σ i j ( 1 + 2 ) = σ i j ( 1 ) + σ i j ( 2 )
ϵ i j ( 1 + 2 ) = ϵ i j ( 1 ) + ϵ i j ( 2 )
Thus, the total elastic energy in the combined state is:
E ( 1 + 2 ) = 1 2 V σ i j ( 1 + 2 ) ϵ i j ( 1 + 2 ) d V
Expanding the integrand:
E ( 1 + 2 ) = 1 2 V σ i j ( 1 ) ϵ i j ( 1 ) + σ i j ( 1 ) ϵ i j ( 2 ) + σ i j ( 2 ) ϵ i j ( 1 ) + σ i j ( 2 ) ϵ i j ( 2 ) d V
This can be separated into the sum of energies for individual states and an interaction term:
E ( 1 + 2 ) = E ( 1 ) + E ( 2 ) + E ( 1 2 )
where the interaction term is:
E ( 1 2 ) = 1 2 V σ i j ( 1 ) ϵ i j ( 2 ) + σ i j ( 2 ) ϵ i j ( 1 ) d V
Colonetti’s theorem:
Colonetti’s theorem asserts that this interaction energy E ( 1 2 ) is zero. To prove this, we proceed as follows. We first express Stress in Terms of Strain
σ i j = C i j k l ϵ k l
where C i j k l is the stiffness tensor of the material. This allows us to write:
σ i j ( 2 ) ϵ i j ( 1 ) = C i j k l ϵ k l ( 2 ) ϵ i j ( 1 )
and similarly:
σ i j ( 1 ) ϵ i j ( 2 ) = C i j k l ϵ k l ( 1 ) ϵ i j ( 2 )
Since C i j k l is symmetric in ( i j ) and ( k l ) , it follows that:
σ i j ( 1 ) ϵ i j ( 2 ) = σ i j ( 2 ) ϵ i j ( 1 )
We shall now Simplify Interaction Energy. Using the above equality, the interaction energy becomes:
E ( 1 2 ) = V σ i j ( 1 ) ϵ i j ( 2 ) d V
Given that σ i j ( 1 ) is purely internal (i.e., it satisfies σ i j , i ( 1 ) = 0 ), we can use Gauss’s theorem to convert the volume integral into a surface integral:
E ( 1 2 ) = V ( σ i j ( 1 ) u j ( 2 ) ) , i d V = S n i σ i j ( 1 ) u j ( 2 ) d S
Here, n i is the outward normal to the surface S. Since σ i j ( 1 ) is internal, the traction t j = σ i j ( 1 ) n i on the external surface is zero:
σ i j ( 1 ) n i = 0 on S
Thus:
E ( 1 2 ) = 0
This concludes the proof of Colonetti’s theorem, showing that the interaction term between internal and applied stress fields vanishes, meaning there is no cross-term in the total elastic energy between these fields.
The total elastic energy E in a finite solid due to an inclusion under applied tractions can be expressed as the sum of the elastic energies due to the applied tractions E A and the elastic energy due to the eigenstrain in the finite solid E F :
E = E A + E F
The elastic energy due to the applied tractions is given by:
E A = 1 2 V σ i j A e i j A d V
Here:
1.
σ i j A is the stress field due to the applied tractions.
2.
e i j A is the strain field corresponding to σ i j A .
This expression comes from the standard formula for elastic energy density 1 2 σ i j e i j , integrated over the volume of the solid. The elastic energy due to the eigenstrain in the finite solid is given by:
E F = 1 2 σ i j I , + σ i j i m ¯ e i j * V 0
Where:
1.
σ i j I , is the stress field due to the inclusion in an infinite medium.
2.
σ i j i m ¯ is the averaged image stress field due to the boundary effects of the finite solid.
3.
e i j * is the eigenstrain in the inclusion.
4.
V 0 is the volume of the inclusion.
This expression represents the interaction energy between the eigenstrain and the stress fields in both the infinite medium and the finite boundary. The enthalpy H is defined as the difference between the total elastic energy E and the work done by the loading mechanism Δ W L M :
H = E Δ W L M
The work done by the loading mechanism can be split into two parts: the work done due to the applied tractions alone Δ W L M A and the interaction term Δ W L M A F :
Δ W L M = Δ W L M A + Δ W L M A F
The work done by the loading mechanism due to the applied tractions is:
Δ W L M A = S e x t σ i j A u j A n i e x t d S
By applying the divergence theorem, this surface integral can be converted into a volume integral:
Δ W L M A = V σ i j A e i j A d V
Recognizing that this is twice the elastic energy E A stored in the system due to the applied tractions, we have:
Δ W L M A = 2 E A
The interaction term between the applied tractions and the eigenstrain-induced fields is given by:
Δ W L M A F = S e x t σ i j A u j F n i e x t d S
This represents the work done by the applied tractions on the displacements generated by the eigenstrain u j F . This interaction term can also be expressed as an integral over the inclusion volume V 0 :
Δ W L M A F = V 0 σ i j A e i j * d V
Since the stress field σ i j A is assumed to be uniform over the inclusion, this integral simplifies to:
Δ W L M A F = σ i j A e i j * V 0
To further explore the interaction energy, consider the integral over the matrix volume V M = V V 0
Δ W L M A F = V M σ i j A e i j * d V
Using Gauss’s theorem, the volume integral over V M can be converted into a surface integral over the inclusion boundary S 0 :
Δ W L M A F = S 0 σ i j A u j F , M σ i j F , M u j A n i o u t d S
Here, u j F , M and σ i j F , M are the displacement and stress fields in the matrix due to the inclusion. The integral in the previous equation can be transformed back into a volume integral over V 0 :
Δ W L M A F = V 0 σ i j A e i j * d V
Recognizing that this is negative of Equation 3.64, we conclude:
Δ W L M A F = e i j * σ i j A V 0

16. Ellipsoidal Inhomogenity

16.1. Application of Eshelby’s Inclusion Solution to Inhomogeneities

Eshelby’s solution, originally developed for inclusions, is applicable to various problems such as inhomogeneities, cracks, and dislocations. This is achieved using the Equivalent Inclusion Method, where an eigenstrain is selected to model the specific problem. This method is particularly effective for ellipsoidal inhomogeneities, where the stress and strain inside these inclusions remain constant.
Example Problem: Liquid-Filled Void in an Infinite Solid
Consider a situation where a volume V 0 is excised from an infinite solid and replaced with a liquid under pressure p 0 . The task is to determine the stress, strain, and displacement fields within the surrounding matrix.
In principle, the problem could be addressed using Green’s functions, as shown in the equation:
u i ( x ) = p 0 δ k j n k G ˜ i j ( x , x ) d x
where G ˜ i j ( x , x ) is the Green’s function for an infinite body with a cavity. However, since the exact expression for G ˜ i j ( x , x ) is unknown, this direct method is impractical.
Eshelby’s Equivalent Inclusion Method
A more feasible approach involves replacing the liquid with an inclusion whose eigenstrain e i j * is chosen so that the internal stress within the inclusion matches that within the liquid (i.e., σ i j I = p 0 δ i j ). Given the constant stress and strain in both the inclusion and the liquid, the required eigenstrain can be determined using Eshelby’s tensor S i j k l . The stress inside the equivalent inclusion can be expressed as:
σ i j I = σ i j c σ i j * = C i j k l ( e k l c e k l * ) = C i j k l ( S k l m n e m n * e k l * )
This leads to the equation:
C i j k l ( S k l m n δ k m δ l n ) e m n * = p 0 δ i j
From this set of six equations, the six unknown components of the equivalent eigenstrain e i j * can be solved.
Elastic Energy Considerations
Once the eigenstrain is known, the displacement on the void surface S 0 can be computed:
u i = u i c = S i j k l e k l * x j
The elastic energy inside the matrix, which must be identical to that in the case where the equivalent inclusion replaces the liquid, is given by:
E = E I + E M = 1 2 σ i j I e i j * V 0
where E I is the energy in the inclusion and E M is the energy in the matrix. Specifically:
E I = 1 2 σ i j I e i j I V 0 = 1 2 σ i j I ( e i j c e i j * ) V 0
And the matrix energy is:
E M = E E I = 1 2 σ i j I e i j c V 0 = 1 2 p 0 S i j k l e k l * V 0
This section thus lays the groundwork for extending Eshelby’s methods to more complex scenarios involving inhomogeneities, using the equivalent inclusion approach to simplify the analysis of elastic fields.

16.2. Transformed Inhomogeneity

We consider an inhomogeneity within an elastic matrix, where the inhomogeneity has different material properties (represented by a stiffness tensor C i j k l ) compared to the matrix (with stiffness tensor C i j k l ). The inhomogeneity undergoes a permanent transformation described by an eigenstrain e i j * . The goal is to determine the stress, strain distribution, and total elastic energy in the solid. The stress inside the inhomogeneity is given by:
σ i j = σ i j c σ i j * = C i j k l ( e k l c e k l * )
where e k l c is the total strain inside the inhomogeneity, and e k l * is the eigenstrain specific to the inhomogeneity. To simplify the problem, we introduce an equivalent homogeneous inclusion, which is assumed to be ellipsoidal, having the same material properties as the matrix C i j k l , but subjected to an effective eigenstrain e i j e f f that ensures the stress and strain inside the equivalent inclusion are identical to those in the inhomogeneity. The stress inside the equivalent inclusion is given by:
σ i j = σ i j c σ i j * = C i j k l ( e k l c e k l * )
where e k l * is the eigenstrain for the equivalent inclusion. For the equivalent inclusion to correctly represent the inhomogeneity, both the stress and total strain must match between the two systems:
σ i j = σ i j and e i j c = e i j c
Thus, substituting the stress and strain equations into these conditions yields:
C i j k l ( e k l c e k l * ) = C i j k l ( e k l c e k l * )
Since e i j c = e i j c , we substitute e k l c = S k l m n e m n * into the equation to get:
C i j k l ( S k l m n e m n * e k l * ) = C i j k l ( S k l m n e m n * e k l * )
Rearranging the equation, we obtain the relation between the actual eigenstrain e k l * and the effective eigenstrain e k l * :
( C i j k l C i j k l ) S k l m n + C i j k l ) e m n * = C i j k l e k l *
This equation allows us to solve for the equivalent eigenstrain e m n * in terms of the eigenstrain e k l * of the transformed inhomogeneity. The stress and strain fields in the matrix and inclusion can now be computed once the effective eigenstrain e m n * is known. The total strain inside the inhomogeneity and the equivalent inclusion are identical:
e i j c = e i j c = S i j k l e k l *
The stress inside the inhomogeneity and the equivalent inclusion are also identical:
σ i j = σ i j = C i j k l ( S k l m n e m n * e k l * )
The elastic energy inside the matrix must be identical in both the transformed inhomogeneity problem and the equivalent inclusion problem:
E M = 1 2 σ i j e i j c V 0
The elastic energy inside the inhomogeneity E I and the equivalent inclusion E I are given by:
E I = 1 2 σ i j e i j I V 0 = 1 2 σ i j ( e i j c e i j * ) V 0
E I = 1 2 σ i j e i j I V 0 = 1 2 σ i j ( e i j c e i j * ) V 0
Thus, the total energy for the solid with a transformed inhomogeneity is:
E = E I + E M = 1 2 σ i j e i j * V 0
This derivation rigorously connects the stress, strain, and energy fields for a transformed inhomogeneity with those of an equivalent inclusion, providing a comprehensive framework for analyzing such problems in elasticity.

16.3. Inhomogeneity under Uniform Applied Loads

Consider a solid containing an inhomogeneity with no eigenstrain. The solid is subjected to external loads, and if it were homogeneous (without the inhomogeneity), the stress and strain fields would be uniform throughout the solid. The primary question is how the presence of the inhomogeneity affects the stress and strain fields. To solve this problem, we construct the stress and strain fields by superimposing two sets of fields:
1.
First Set: Suppose the solid with the inhomogeneity is subjected to a uniform strain e i j A . The stress fields inside the matrix and inhomogeneity are given by:
σ i j A = C i j k l e k l A and σ i j A = C i j k l e k l A
However, this stress field does not satisfy equilibrium conditions unless a body force T j = ( σ i j A σ i j A ) n i is applied on the surface S 0 of the inhomogeneity.
2.
Second Set: To restore equilibrium, apply a body force F j = T j on S 0 . The corresponding stress and strain fields due to this body force are σ i j c and e i j c .
The elastic stress field inside the inhomogeneity, resulting from the superposition of these two sets of fields, is:
σ i j I = σ i j A + σ i j c = C i j k l ( e k l A + e k l c )
The total strain field inside the inhomogeneity is:
e i j I = e i j A + e i j c
At the same time, consider an equivalent inclusion with eigenstrain e i j * in a solid under the same uniform applied load. The elastic stress field inside this inclusion is:
σ i j I = σ i j A + σ i j c σ i j * = C i j k l ( e k l A + e k l c e k l * )
The total strain field inside the inclusion is:
e i j I = e i j A + e i j c
To ensure that the elastic stress and total strain match between the inhomogeneity and the inclusion problems, the following conditions must hold:
C i j k l ( e k l A + e k l c ) = C i j k l ( e k l A + e k l c e k l * )
e i j A + e i j c = e i j A + e i j c
From the second equation, e i j c = e i j c , which when substituted into the first equation gives:
C i j k l ( e k l A + e k l c ) = C i j k l ( e k l A + e k l c e k l * )
This can be rearranged to solve for the effective eigenstrain e k l * :
( C i j k l C i j k l ) S k l m n + C i j k l ) e m n * = ( C i j k l C i j k l ) e k l A
This equation expresses the equivalent eigenstrain e m n * in terms of the applied strain e k l A and the difference in stiffness tensors C i j k l and C i j k l .
Application of the Feynman-Hellmann Theorem:
The Feynman-Hellmann theorem is applied in the context of deriving the total elastic energy E and enthalpy H of the inhomogeneous solid. This theorem relates the variation in the total energy of the system to the variation in the applied field.
Specifically, the Feynman-Hellmann theorem leads to the following result for the change in enthalpy Δ H due to the presence of the inhomogeneity:
Δ H = 1 2 ( C i j k l C i j k l ) e i j A e k l I V 0
where V 0 is the volume of the inhomogeneity, and e k l I is the strain field within the inhomogeneity, which includes contributions from both the applied strain e k l A and the correction due to the presence of the inhomogeneity itself. The Feynman-Hellmann theorem in this context is derived by considering the total energy E of the system as a function of the stiffness tensors and the applied strain. The theorem states:
E C i j k l = E e k l A = 1 2 σ i j e k l A
This expression allows us to compute the change in energy due to a small perturbation in the applied strain or stiffness tensor. In this case, the enthalpy H is related to the total energy by:
Δ H = E C i j k l Δ C i j k l
Substituting the stress-strain relation into this expression gives the final form of the change in enthalpy as:
Δ H = 1 2 δ C i j k l e i j A e k l A V 0
where δ C i j k l = C i j k l C i j k l is the difference in the stiffness tensors of the inhomogeneity and the matrix. To compute the total elastic energy of the system, consider a reversible path where the inhomogeneity is subjected to a uniform strain e i j A . The elastic energy for this state is given by:
E 1 = 1 2 σ i j A e i j A V M + 1 2 σ i j A e i j A V 0 = 1 2 σ i j A e i j A V + 1 2 ( σ i j A σ i j A ) e i j A V 0
Where V M is the volume of the matrix, V 0 is the volume of the inhomogeneity, and V is the total volume of the solid. Gradually removing the body force results in a final energy E 2 , which is the desired solution. The total elastic energy E of the inhomogeneous solid is given by:
E = E 2 = E 1 + Δ W 12
where Δ W 12 represents the work done during the transformation. This accounts for both internal and external work contributions.
Using the relation between the equivalent eigenstrain and the applied strain, the total elastic energy and enthalpy of the system can be derived, yielding the following final expressions:
E = 1 2 σ i j A e i j A V 1 2 ( σ i j A σ i j A ) e i j I V 0
The enthalpy H is obtained by subtracting the work done by the external loading mechanism from the internal energy:
H = E Δ W L M
This rigorous derivation, combined with the application of the Feynman-Hellmann theorem, provides a comprehensive framework for understanding the stress, strain, and energy fields in an inhomogeneous material subjected to uniform loads.

17. Cracks

17.1. Ellipsoidal Void

17.1.1. Stress and Strain Relations for an Ellipsoidal Void

Given that the stiffness tensor C i j k l of the inhomogeneity approaches zero, the inhomogeneity becomes a void. The stress field inside the void must be zero, so the stress-strain relation becomes:
0 = C i j k l ( e k l A + e k l c ) = C i j k l ( e k l A + e k l c e k l * )
This equation reduces to:
e k l A + e k l c = e k l *
The applied strain e k l A plus the strain inside the void e k l c must balance with the eigenstrain e k l * of the equivalent inclusion.

17.1.2. Eigenstrain in the Void

The equivalent eigenstrain, which generates no stress inside the void, is related to the applied stress:
e i j * = 1 C i j k l σ k l A
The total strain in the void, given by the eigenstrain of the equivalent inclusion, is:
e i j I = e i j A + e i j c = e i j *
This ensures zero stress within the void since the total strain equals the eigenstrain.

17.1.3. Enthalpy Calculation

The change in enthalpy Δ H of the system due to the presence of the void is given by:
Δ H = 1 2 ( C i j k l C i j k l ) e i j A e k l I V 0 = 1 2 σ i j A e i j * V 0
The eigenstrain e i j * can be explicitly solved using the applied stress σ i j A , reinforcing that the stress inside the void cancels the applied stress.

17.2. Penny-Shaped Crack

A penny-shaped crack is modeled as an ellipsoidal void where the axis c approaches zero while the other two axes, a and b, are equal, i.e., a = b and c 0 . This configuration leads to a circular disk-like crack in an infinite elastic medium.

17.2.1. Eshelby’s Tensor for Penny-Shaped Crack

Derivation of Eshelby’s Tensor Components:
Eshelby’s tensor S i j k l relates the eigenstrain e k l * within the inclusion to the resulting strain e i j c in the material:
e i j c = S i j k l e k l *
For an ellipsoidal inclusion, Eshelby’s tensor is generally a function of the aspect ratios of the inclusion. In our case, for a penny-shaped crack where c 0 and a = b , the non-zero components of the Eshelby tensor are determined as follows:
1. Component S 1111 : Consider the geometry where the crack lies in the x 1 - x 2 plane. The component S 1111 is given by:
S 1111 = π ( 13 8 ν ) 32 ( 1 ν ) c a
This expression is derived by evaluating the Eshelby tensor in the limit as c 0 , which simplifies the general expressions for the tensor components.
2. Component S 1133 : Similarly, for the component S 1133 , which relates the eigenstrain in the x 3 direction to the resulting strain in the x 1 direction, we have:
S 1133 = π ( 7 8 ν ) 16 ( 1 ν ) c a
The derivation follows by considering the effect of the eigenstrain in the x 3 direction on the strain in the x 1 direction for the penny-shaped crack.
3. Component S 1212 : The component S 1212 , which is related to shear deformation, is given by:
S 1212 = π ( 3 4 ν ) 8 ( 1 ν ) c a
This component is derived by analyzing the shear response of the material due to the eigenstrain in the x 1 - x 2 plane.

17.2.2. Eigenstrain and Stress in Penny-Shaped Crack

Derivation of Eigenstrain e i j * :
For a penny-shaped crack under an applied tensile stress σ 33 A , the eigenstrain components are related to the applied stress and Eshelby’s tensor:
σ 11 A = 2 μ 1 ν + 13 μ π c 16 ( 1 ν ) a e 11 * +
To derive this, we start by considering the stress-strain relationship in the material:
σ i j = C i j k l e k l
where C i j k l is the stiffness tensor for the isotropic material. The total strain e i j in the material is the sum of the applied strain e i j A and the eigenstrain e i j * . Therefore, the stress is:
σ i j = C i j k l ( e k l A + e k l * )
Given that the stress inside the crack is zero, the applied stress must be balanced by the eigenstrain-induced stress:
σ i j A = C i j k l e k l *
Now, for the penny-shaped crack, the eigenstrain component e 11 * is derived considering the symmetry and the specific components of Eshelby’s tensor:
σ 11 A = 2 μ 1 ν + 13 μ π c 16 ( 1 ν ) a e 11 * + ( other terms )
where μ is the shear modulus, and the terms involving Eshelby’s tensor components are included to account for the interaction between the applied stress and the eigenstrain.

17.2.3. Limiting Behavior of Eigenstrain e 33 *

As c 0 , the eigenstrain e 33 * in the direction normal to the crack tends to infinity, but the product e 33 * c remains finite. This leads to:
e 33 * = 2 ( 1 ν ) a μ π σ 33 A
This expression is derived by considering the balance of forces and the boundary conditions on the crack surface, ensuring that the stress inside the crack is zero.

17.3. Energy Considerations and Griffith Criterion

Derivation of Enthalpy Change Δ H :
The change in enthalpy Δ H due to the presence of the crack is calculated using the energy associated with the eigenstrain:
Δ H = 1 2 V σ i j e i j d V
Substituting the expression for σ i j and integrating over the volume of the crack, we get:
Δ H = 4 ( 1 ν ) 3 μ ( σ 33 A ) 2 a 3
This result quantifies the energy difference between the cracked and uncracked states.
Derivation of Griffith Criterion:
The Griffith criterion for crack growth is derived from the Gibbs free energy Δ G , which includes the enthalpy change Δ H and the surface energy γ associated with the crack surfaces:
Δ G = Δ H + 2 π γ a 2
To find the condition for crack growth, we set the derivative of Δ G with respect to a to zero:
d Δ G d a = d Δ H d a + d ( 2 π γ a 2 ) d a = 0
d Δ H d a + 4 π γ a = 0
Substituting the expression for Δ H , we obtain:
8 ( 1 ν ) 3 μ ( σ 33 A ) 2 a 2 + 4 π γ a = 0
Simplifying, we get the critical stress σ 33 A required for crack propagation:
σ 33 A = π μ γ ( 1 ν ) a
This is the Griffith criterion, which determines the stress at which the crack will grow, leading to material failure.

17.4. Slit-Like Crack

17.4.1. Derivation of Eshelby’s Tensor in the Slit-Like Crack Limit

Given the geometry of the slit-like crack, the limits c and b 0 are applied to Eshelby’s tensor for an ellipsoidal inclusion in an isotropic medium. The Eshelby tensor components are given as follows:
S 1111 = 1 2 ( 1 ν ) b 2 + 2 a b ( a + b ) 2 + ( 1 2 ν ) b a + b
S 2222 = 1 2 ( 1 ν ) a 2 + 2 a b ( a + b ) 2 + ( 1 2 ν ) a a + b
S 1122 = 1 2 ( 1 ν ) b 2 ( a + b ) 2 ( 1 2 ν ) b a + b
S 2211 = 1 2 ( 1 ν ) a 2 ( a + b ) 2 ( 1 2 ν ) a a + b
As b 0 , these simplify to:
S 1111 1 2 ( 1 ν )
S 2222 a 2 2 ( 1 ν ) ( a + b ) 2
S 1122 1 2 ν 2 ( 1 ν )
S 2211 1 2 ν 2 ( 1 ν )

17.4.2. Equivalent Eigenstrain

For the slit-like crack, we assume that the eigenstrain tensor components are given by e 11 * and e 22 * , with all other components being zero due to the geometry and loading conditions. The equivalent eigenstrain in the limit b 0 and c is obtained by solving the following system of equations:
σ 11 A = 2 a 2 + a b ( 1 ν ) ( a + b ) 2 μ e 11 * a b ( 1 ν ) ( a + b ) 2 μ e 22 *
σ 22 A = a b ( 1 ν ) ( a + b ) 2 μ e 11 * a b + 2 b 2 ( 1 ν ) ( a + b ) 2 μ e 22 *
Taking the limit b 0 and assuming σ 11 A = 0 , we have:
0 = 2 μ 1 ν e 11 * b μ ( 1 ν ) a e 22 *
σ 22 A = b μ ( 1 ν ) a e 11 * b μ ( 1 ν ) a e 22 *
Defining e * = lim b 0 e 22 * b and allowing e 11 * to remain finite, we solve these equations to obtain:
e * = ( 1 ν ) a μ σ 22 A
e 11 * = ( 1 ν ) σ 22 A 2 μ

17.4.3. Griffith Criterion for Slit-Like Crack Growth

Using the derived eigenstrain, the total enthalpy per unit length of the crack can be computed as:
Δ H / c = ( 1 ν ) π ( σ 22 A ) 2 a 2 2 μ
The driving force per unit length for crack growth is then:
f t o t a = Δ G / c a = ( 1 ν ) π ( σ 22 A ) 2 a μ 4 γ
Setting the critical condition f t o t a = 0 gives the Griffith criterion:
σ 22 A = 4 μ γ ( 1 ν ) π a

17.4.4. Stress Intensity Factors and Crack Tip Fields

To evaluate the stress intensity factors, we consider the stress field near the crack tip, denoted by r (distance from the crack tip) and θ (polar angle). The stress field is singular as r 0 and follows:
σ r r = K I 2 π r 5 4 cos θ 2 1 4 cos 3 θ 2
σ θ θ = K I 2 π r 3 4 cos θ 2 + 1 4 cos 3 θ 2
σ r θ = K I 2 π r 1 4 sin θ 2 + 1 4 sin 3 θ 2
The stress intensity factor K I is then defined as:
K I = σ r r 2 π r as r 0
For the slit-like crack under uniform tension, the stress intensity factor simplifies to:
K I = π a σ 22 A

17.5. Flat Ellipsoidal Crack

A flat ellipsoidal crack represents a case between the two extremes of penny-shaped and slit-shaped cracks. This type of crack has an ellipsoidal shape where a > b and c 0 . The goal is to understand whether the crack will tend to become more elongated (slit-like) or less elongated (penny-shaped).

17.5.1. Eigenstrain Calculation

Let us consider a simple tensile stress applied in the σ 33 A direction, with all other components of the applied stress being zero. The key idea is to keep the product e 33 * c constant as c 0 . The solution for the eigenstrain e * is given by:
e * = ( 1 ν ) b μ E ( k ) σ 33 A
where E ( k ) is the elliptic integral of the second kind, defined as:
E ( k ) = 0 π / 2 1 k 2 sin 2 w d w
and k is given by:
k = 1 b 2 a 2
The eigenstrain e * is a function of the applied stress σ 33 A , Poisson’s ratio ν , the semi-minor axis b, and the elliptic integral E ( k ) .

17.5.2. Enthalpy Change Calculation

The extra enthalpy due to the presence of the crack is calculated as:
Δ H = 1 2 σ 33 A e 33 * · 4 π 3 a b c = 2 π 3 σ 33 A e * a b
Substituting the expression for e * :
Δ H = 2 π ( 1 ν ) 3 μ a b 2 E ( k ) ( σ 33 A ) 2
This expression captures the change in enthalpy due to the crack and shows its dependence on the crack dimensions a and b, the applied stress σ 33 A , and the material properties.

17.5.3. Gibbs Free Energy and Griffith Criterion

The Gibbs free energy, which includes the surface energy of the crack, is given by:
Δ G = 2 π ( 1 ν ) 3 μ a b 2 E ( k ) ( σ 33 A ) 2 + 2 π γ a b
To find the conditions for crack growth, we differentiate the Gibbs free energy with respect to a and b:
Δ G a = 0
Δ G b = 0
These conditions provide the critical stresses required for crack growth in the a and b directions:
σ 33 A , a = 3 μ γ k 2 E 2 ( k ) b ( 1 ν ) ( 1 + 2 k 2 ) E ( k ) + ( 1 k 2 ) F ( k )
σ 33 A , b = 3 μ γ k 2 E 2 ( k ) b ( 1 ν ) ( 1 + k 2 ) E ( k ) ( 1 k 2 ) F ( k )
Here, F ( k ) is the elliptic integral of the first kind:
F ( k ) = 0 π / 2 d w 1 k 2 sin 2 w
These expressions indicate whether the crack will grow in the a or b direction, depending on which stress component reaches its critical value first. If σ 33 A , b < σ 33 A , a , the crack will tend to become more penny-shaped. Otherwise, it will become more slit-like.

17.6. Crack Opening Displacement: Rigorous Derivation

We consider the elastic fields (displacement, strain, and stress) of a slit-like crack under tensile loading stress σ 22 A . The goal is to determine the crack opening displacement d ( x ) as a function of position x.
Let d ( x ) be the distance between the crack faces as a function of x. In a purely elastic model, d ( ± a ) = 0 , i.e., the crack tip opening displacement is zero. To determine the displacements along the crack face, we consider an equivalent inclusion problem. The displacement field u j ( x ) for an inclusion is given by:
u j ( x ) = e i j * x j
where e i j * is the eigenstrain of the inclusion. Since we are dealing with a slit-like crack, the displacement in the x-direction ( u 1 ) is zero, and the displacement in the y-direction ( u 2 ) on the crack face is given by:
u 2 = e 22 * y
Here, e 22 * is the eigenstrain component associated with the opening of the crack. The equivalent inclusion is an ellipse with semi-axes a and b (with b 0 in the slit-like crack case). The relationship between x and y on the crack surface is:
x 2 a 2 + y 2 b 2 = 1
The displacement field on the upper surface of the crack at x [ a , a ] is:
u 2 ( x ) = e 22 * b · y b = e 22 * 1 x 2 a 2
The eigenstrain e 22 * is related to the applied stress σ 22 A by:
e 22 * = σ 22 A a ( 1 ν ) μ
where μ is the shear modulus and ν is Poisson’s ratio. Substituting the eigenstrain into the displacement expression:
u 2 ( x ) = σ 22 A a ( 1 ν ) μ 1 x 2 a 2 = σ 22 A ( 1 ν ) μ a 2 x 2
u 2 ( x ) = σ 22 A ( 1 ν ) μ a 2 x 2
The crack opening displacement d ( x ) is twice the displacement u 2 ( x ) :
d ( x ) = 2 σ 22 A ( 1 ν ) μ a 2 x 2
This is the crack opening displacement in plane strain. In plane stress conditions, the displacement is modified as:
d ( x ) = 2 σ 22 A μ ( 1 + ν ) a 2 x 2
Using the expression for d ( x ) , we calculate the enthalpy of the crack by measuring the work done while opening up the crack. In plane stress, the enthalpy change is given by:
Δ H c = 1 2 a a d ( x ) σ 22 A d x
Substituting the expression for d ( x ) and evaluating the integral:
Δ H c = 1 2 σ 22 A · 2 σ 22 A ( 1 ν ) μ a a a 2 x 2 d x
The integral evaluates to:
a a a 2 x 2 d x = π a 2 2
Thus, the enthalpy is:
Δ H c = 1 ν 2 μ ( σ 22 A ) 2 π a 2
This matches the previously calculated enthalpy, confirming the correctness of the derived crack opening displacement.

17.7. Stress Intensity Factors

Let r be the distance to the crack tip. The stress field in the vicinity of the crack tip exhibits a singularity of the form:
σ ( r ) 1 r .
The stress intensity factor K I is defined as:
K I = lim r 0 σ ( r ) 2 π r
The three modes of crack opening are:
  • Mode I: Tensile mode, K I ,
  • Mode II: In-plane shear mode, K I I ,
  • Mode III: Out-of-plane shear mode, K I I I .
Using Eshelby’s tensor, the stress intensity factors can be related to the eigenstrain inside an inclusion. The auxiliary tensor D i j k l is:
D i j k l = a b 2 π 0 2 π ( zz ) i j 1 z k z l κ ( γ ) β 2 d θ
with κ ( γ ) , β , and γ defined as functions of the geometry. For x > a , the stress field on the crack plane is derived using the Eshelby tensor and auxiliary tensor components. The stress intensity factor for Mode I is:
K I = π a σ 22 A ,
where σ 22 A is the applied stress normal to the crack plane, and a is the half-length of the crack.

17.8. Another Derivation of Crack Extension Force

Let’s Define the Problem and Initial Conditions. Consider a two-dimensional crack under uniform tension σ 22 = σ A . The crack half-size is a, and we analyze the situation where the crack extends by a small amount δ a , making the new crack half-size a + δ a .
Initially, additional traction forces T j ± are applied on the surfaces of the crack in the region [ a , a + δ a ] and [ a δ a , a ] to keep the crack shape unchanged. The traction forces are then removed gradually, allowing the crack to extend freely. The work done by these forces corresponds to the change in system enthalpy δ H . Let’s now Compute the Work Done by Traction Forces. The applied traction forces on the surfaces of the crack are given by:
T j + ( x ) = σ j 2 ( x ) , T j ( x ) = σ j 2 ( x )
The crack opening displacement is defined as:
d ( x ) = u 2 u 2 +
The change in enthalpy δ H is computed by the work done by the traction forces over the region [ a , a + δ a ] :
δ H = a a + δ a ( T j + u j + + T j u j ) d x
Since T j = T j + , the equation simplifies to:
δ H = 2 a a + δ a T j + u j + d x
Substituting T j + = σ 22 ( x ) and using the expression for d ( x ) :
δ H = 2 a a + δ a σ 22 ( x ) d ( x ) d x
We shall now Evaluate the Crack Opening Displacement d ( x ) . The crack opening displacement for a two-dimensional crack under uniform tension is:
d ( x ) = 2 σ A ( 1 ν ) μ a 2 x 2
For small δ a , d ( x ) near x = a can be approximated as:
d ( x ) 2 σ A ( 1 ν ) μ 2 a δ a
Using the above approximation in the expression for δ H :
δ H = 2 a a + δ a σ A 2 σ A ( 1 ν ) μ 2 a δ a d x
Simplifying the integral and keeping only terms linear in δ a :
δ H = 4 σ A 2 ( 1 ν ) μ δ a a a + δ a x x 2 a 2 d x
The integral can be evaluated, and in the limit δ a a :
δ H = 4 σ A 2 ( 1 ν ) μ · π δ a ( 2 a + δ a ) 4 = π σ A 2 ( 1 ν ) μ δ a
Let’s Derive the Crack Extension Force. The crack extension force is defined as:
f = δ H δ a = π σ A 2 ( 1 ν ) μ
This result matches the crack extension force obtained by previous methods, confirming the correctness of this alternative derivation. Thus, the very rigorous derivation leads to the final expression for the crack extension force:
f = ( 1 ν ) σ A 2 π a μ
This matches the previously derived expressions, demonstrating consistency across different approaches.

17.9. J-integral as Driving Force

The J-integral, denoted as J i , in a three-dimensional elastic medium represents the force on an elastic singularity in the i-th direction:
J i = S w n i T j u j , i d S
where w is the strain energy density, n i is the unit normal vector to the surface S, T j is the traction vector, and u j , i is the displacement gradient. In two dimensions, for a crack along the x-axis, the J-integral simplifies to:
J = Γ w d y T u x d s
where Γ is a contour encircling the crack tip. The strain energy density w is:
w = 0 e i j σ i j d e i j
The total enthalpy H of the system is:
H = E S T T j ext u j d S
where E is the total strain energy:
E = V w d V
The driving force f i on the singularity at ξ i is:
f i = δ H δ ξ i
To compute f i , we determine the variation of total enthalpy δ H as the crack tip moves by δ ξ i :
δ H = V δ w d V S T T j ext δ u j d S
Consider a sub-volume V 0 with surface S 0 , and V E = V V 0 :
V E δ w d V = V E σ i j δ e i j d V
V E δ w d V = V E σ i j δ u j , i d V = V E ( σ i j δ u j ) , i d V
Using Gauss’s Theorem:
V E δ w d V = S T T j ext δ u j d S S 0 T j δ u j d S
Substituting into δ H :
δ H = V δ w d V S T T j ext δ u j d S = V 0 w ξ i δ ξ i d V + S 0 T j u j ξ i δ ξ i d S
Therefore, the driving force f i is:
f i = δ H δ ξ i = S 0 w n i T j u j , i d S = J i

17.10. Invariance of J-Integral

The J-integral is a fundamental quantity in fracture mechanics, representing the driving force on a crack. The invariance of the J-integral with respect to the surface or contour on which it is evaluated is a critical property that makes it a powerful tool in the analysis of crack problems. The J-integral in its general three-dimensional form is defined as:
J k = V 0 w x k d V S 0 T j u j x k d S
where:
  • w is the strain energy density,
  • T j is the traction vector,
  • u j is the displacement field,
  • x k is the spatial coordinate, and
  • S 0 is a surface surrounding the crack tip.
To prove the invariance of the J-integral, we start by considering the derivative of the strain energy density with respect to the spatial coordinate x k :
w x k = w e i j e i j x k = σ i j e i j x k
Using the relationship between strain and displacement gradients:
e i j = 1 2 u i x j + u j x i
we can express the derivative of the strain with respect to x k as:
e i j x k = 1 2 2 u i x j x k + 2 u j x i x k
Substituting this into the expression for w x k , we obtain:
w x k = σ i j 2 u j x k x i = x i σ i j u j x k
The equilibrium condition σ i j , i = 0 has been used in the last step. We now consider the J-integral over a closed surface S 0 containing no defects. Using the expression for w x k , the J-integral becomes:
J k = V 0 x i σ i j u j x k d V S 0 T j u j x k d S
Applying the divergence theorem to the volume integral:
J k = S 0 n i σ i j u j x k T j u j x k d S
where n i is the outward normal to the surface S 0 . To prove that the J-integral is invariant with respect to the contour or surface used in its evaluation, consider two contour lines Γ 1 and Γ 2 around the crack tip in a 2-dimensional problem. If we take a complete contour Γ = Γ 1 + B + Γ 2 + B that encloses no singularities, the J-integral over this contour must be zero:
J ( Γ ) = J ( Γ 1 ) J ( Γ 2 ) + J ( B + ) + J ( B )
Since J ( B + ) = J ( B ) (because d y = 0 and T = 0 on the crack faces), we conclude that:
J ( Γ 1 ) = J ( Γ 2 )
This demonstrates the invariance of the J-integral.

17.11. Applications of J-Integral

Consider a very long solid slab with a crack in the middle. The top and bottom surfaces are subjected to constant displacement boundary conditions, and the left and right ends are subjected to zero surface traction boundary conditions. The J-integral in two dimensions is expressed as:
J = Γ w d y T u x d s
where Γ is the contour surrounding the crack tip, w is the strain energy density, and T is the traction vector.
1.
On S 2 and S 4 , d y = 0 and u x = 0 , hence the contributions to J are zero.
2.
On S 1 and S 5 , w = 0 and u x = 0 , leading to zero contributions as well.
3.
On S 3 , w = w and u x = 0 .
Therefore, the total J-integral becomes:
J = w h
where h is the height of the slab.
Consider a two-dimensional crack with a blunt tip. The J-integral for this configuration simplifies to:
J = Γ w d y
This integral represents the average strain energy density around the crack tip.
Consider a mode-I crack with stress intensity factor K I . The J-integral is evaluated over a circular contour Γ with radius r in the limit r 0 . The stress fields around the crack are given by the leading singular terms in polar coordinates ( r , θ ) :
σ r r = K I 2 π r 5 4 cos θ 2 1 4 cos 3 θ 2 +
σ θ θ = K I 2 π r 3 4 cos θ 2 + 1 4 cos 3 θ 2 +
σ r θ = K I 2 π r 1 4 sin θ 2 + 1 4 sin 3 θ 2 +
The strain energy density w is:
w = 1 2 σ θ θ e θ θ + σ r r e r r + 2 σ r θ e r θ
Substituting the stress fields into the expression for w, the J-integral can be evaluated as:
J = lim r 0 Γ w d y
The final result is:
J = K I 2 E
where E is the effective modulus.

18. Dislocations

18.1. Introduction to Dislocations

The concept of dislocations was introduced by Volterra in 1907 as a mathematical construct to model discontinuities in a solid material. Dislocations are line defects within a crystal structure, where atoms are misaligned. These defects are crucial for understanding the mechanical behavior of materials, particularly their plasticity.
Dislocations remained a purely theoretical construct until the 1930s, when Taylor, Orowan, and Polanyi independently proposed that dislocations are responsible for crystal plasticity. They suggested that the motion of dislocations under stress could explain the actual yield stress observed in metals, which was much lower than previous theoretical predictions.
The theoretical strength of a perfect crystal, τ th , is the stress required to cause plastic shear deformation across an entire slip plane. This theoretical stress is much higher than the experimentally observed yield stress, which is due to the presence of dislocations.
Let’s explore a Mathematical Representation of Dislocations. Consider a perfect crystal subject to shear stress τ along a plane A, as illustrated in Figure 7.1. The shear stress τ ( x ) required to displace the upper half of the crystal by a distance x relative to the lower half is a periodic function due to the crystal’s atomic structure:
τ ( x ) = μ b 2 π a sin 2 π x b
where:
  • μ is the shear modulus,
  • b is the magnitude of the Burgers vector (which represents the magnitude of lattice distortion),
  • a is a constant related to the atomic spacing.
The maximum shear stress, known as the theoretical critical shear stress τ th , occurs when x = b / 2 :
τ th = μ b 2 π a
This theoretical critical shear stress τ th is significantly higher than the experimentally observed yield stress in metals. The discrepancy arises because real crystals contain dislocations, which lower the stress required to move atomic planes relative to each other. The experimentally measured yield stress is much lower than the theoretical prediction because dislocations provide a mechanism for plastic deformation at much lower stress levels. The movement of dislocations through the crystal lattice under applied stress enables plastic deformation to occur more easily, thus reducing the yield stress.
This understanding revolutionized the field of materials science, providing insights into why materials deform plastically under much lower stresses than would be expected from a perfect crystal model.

18.2. Dislocation’s Effects on Mechanical Properties

Dislocations play a critical role in the mechanical behavior of materials, especially in plastic deformation. When a material is subjected to stress, dislocations move, enabling the material to deform plastically. This section rigorously examines how dislocations affect the mechanical properties of crystals, particularly metals and semiconductors.
The stress-strain curve of a crystal is linear up to the yield stress, beyond which dislocations begin to move, and plastic deformation occurs. As plastic deformation progresses, the length of dislocations within the crystal increases, necessitating higher stresses for continued deformation. This phenomenon is known as work hardening.

18.2.1. Orowan’s Law

One of the key relationships describing the plastic deformation due to dislocations is Orowan’s law, which relates the plastic strain rate ϵ ˙ pl to the dislocation density ρ , the Burgers vector b, and the average dislocation velocity v:
ϵ ˙ pl = ρ b v
where:
  • ρ is the mobile dislocation density (in units of m 2 ),
  • b is the Burgers vector,
  • v is the average dislocation velocity.
Derivation of Orowan’s Law:
Orowan’s law can be derived using Betti’s theorem. The plastic strain rate ϵ ˙ pl is proportional to the rate at which dislocations traverse a given area. Consider a volume element of area A through which dislocations move. If n dislocations pass through A per unit time, the plastic strain rate is given by:
ϵ ˙ pl = b × n A
Since n = ρ v A , where ρ is the dislocation density and v is the velocity, we have:
ϵ ˙ pl = ρ b v
This is Orowan’s law. Let’s analyze the Stress-Strain Curve Behavior in BCC Metals. For body-centered cubic (BCC) metals such as molybdenum, the stress-strain curve under uniaxial tension at a constant strain rate typically shows three stages of deformation:
  • Stage I: Immediately after yielding, plastic deformation occurs with little increase in applied stress. Dislocations primarily glide on parallel planes with minimal interaction.
  • Stage II: At higher deformation, the slope of the stress-strain curve increases, indicating work hardening. Dislocations on several non-parallel slip planes interact, blocking each other’s motion and forming dense, entangled structures. The dislocation density increases significantly.
  • Stage III: The hardening rate decreases as recovery mechanisms begin to annihilate dislocations, leading to a saturation in dislocation density.
Dislocations also influence fracture behavior. In ductile materials, a crack tip can nucleate many dislocations, which shield and blunt the crack tip, leading to a higher critical strain energy release rate J c for crack propagation and higher fracture toughness. Additionally, dislocations can initiate fracture, particularly during fatigue processes. Under cyclic loading, dislocations multiply and can form pile-ups with high local stresses, leading to crack nucleation even in ductile materials.

18.3. Elastic Fields of a Dislocation Loop

A dislocation loop is a closed dislocation line in a crystal lattice that generates elastic fields within the material. The elastic fields associated with the dislocation loop include stress, strain, and displacement fields. These fields can be derived using continuum mechanics and elasticity theory.
Consider a dislocation loop L in an elastic medium, characterized by a Burgers vector b . The displacement field u ( x ) at a point x due to the dislocation loop can be derived using the Green’s function approach:
u i ( x ) = L b j G i j ( x x ) d L ( x )
where G i j ( x x ) is the Green’s function, representing the displacement at point x due to a unit force applied at point x . The Green’s function for an infinite isotropic medium is expressed as:
G i j ( x x ) = 1 8 π ( 1 ν ) μ δ i j 1 r + ( x i x i ) ( x j x j ) r 3
where r = | x x | is the distance between the points x and x . The strain field e i j ( x ) associated with the dislocation loop is obtained by differentiating the displacement field:
e i j ( x ) = 1 2 u i ( x ) x j + u j ( x ) x i
Substituting the displacement field expression, the strain field becomes:
e i j ( x ) = 1 2 L b k G i k ( x x ) x j + G j k ( x x ) x i d L ( x )
Let’s analyze the Stress Field Due to a Dislocation Loop. The stress field σ i j ( x ) is related to the strain field through Hooke’s law:
σ i j ( x ) = C i j k l e k l ( x )
where C i j k l is the fourth-order elasticity tensor for an isotropic material:
C i j k l = λ δ i j δ k l + μ ( δ i k δ j l + δ i l δ j k )
Substituting the expression for the strain field into Hooke’s law, we get:
σ i j ( x ) = 1 2 L C i j k l b k G k l ( x x ) x l + G j l ( x x ) x k d L ( x )
Due to the symmetry of the problem and the properties of the Green’s function, the expressions for the strain and stress fields can be further simplified.
For complex dislocation loop geometries, the integrals in the expressions for the displacement, strain, and stress fields are often evaluated numerically. The loop is discretized into segments, and the fields are computed as the sum of contributions from each segment:
σ i j ( x ) n = 1 N σ i j seg ( x , x n )

18.4. Self Energy of a Dislocation Loop

The self-energy of a dislocation loop refers to the energy stored in the elastic fields due to the dislocation itself. It is an important quantity as it influences the mechanical behavior and stability of dislocations within the material. This derivation rigorously follows the principles of elasticity theory to compute the self-energy of a dislocation loop.
The self-energy of a dislocation loop E can be evaluated by integrating the strain energy density w over the volume V of the material:
E = V w d V
For a linear elastic material, the strain energy density is given by:
w = 1 2 σ i j e i j
where σ i j and e i j are the stress and strain tensors, respectively.
An alternative and more elegant method to calculate the self-energy is by considering the reversible work done to create the dislocation loop. Imagine creating the dislocation loop by applying traction forces F j + and F j on the surfaces S + and S of the loop, respectively. These surfaces are displaced by b , the Burgers vector, relative to each other. The work done W to create the dislocation loop is given by:
W = 1 2 S + F j + u j + d S + 1 2 S F j u j d S
Using the relation u j + u j = b j , the above expression simplifies to:
W = 1 2 S σ k j n k + b j d S
The energy of a dislocation loop obtained from linear elasticity theory is actually singular (infinite) without a proper truncation scheme. This is because, at the core of the dislocation, the strain fields become very large, leading to a divergent integral for the self-energy. To address this, a core cutoff radius is introduced, truncating the fields at a small distance from the dislocation line to avoid the singularity.
For practical computations, especially in numerical simulations, dislocation loops are often represented by a set of connected straight dislocation segments. The stress field from each segment only has physical meaning when summed over the entire loop. The total stress field of the dislocation loop is obtained by summing over the stress fields of individual segments:
σ i j Loop = n = 1 N σ i j seg ( x ( n ) , x ( n + 1 ) , b )
The stress and displacement fields of a dislocation loop in isotropic elasticity can be reduced to line integrals over the dislocation line. For example, the displacement gradients can be expressed as:
u i , j elastic ( x ) = L ϵ j n h C k l m n b m v h G i k , l ( x x ) d S ( x )
This approach, known as Mura’s formula, is valid for evaluating fields around a complete loop and represents the continuous distribution of the dislocation’s influence.

18.5. Force on a Dislocation

The force acting on a dislocation line is a fundamental concept in dislocation theory, as it determines how dislocations move within a crystal lattice under applied stresses. This section rigorously derives the force on a dislocation using the principles of energy variation and the Peach-Koehler force formulation.
Let’s now analyze the Energy Variation and Virtual Displacement. Consider a dislocation loop L with line direction ξ . Let the loop undergo a small virtual displacement δ r ( x ) , where δ r ( x ) · ξ ( x ) = 0 because a line moving along itself has no physical consequence. The energy change δ E due to this displacement can be expressed as:
δ E = L f ( x ) · δ r ( x ) d L ( x )
where f ( x ) is the line force (per unit length) on the dislocation loop L. The force f ( x ) can be found by differentiating the total energy E of the system with respect to the virtual displacement δ r ( x ) .
The total energy E of a system of N dislocation loops can be written as the sum of the loop self-energies E i and the interaction energies W i j between the loops:
E = i = 1 N E i + i = 1 N j = i + 1 N W i j
To calculate the force on a particular loop L 1 , we need to compute the variation of the total energy with respect to the virtual displacement δ r 1 ( x ) of loop L 1 :
f 1 = δ E δ r 1 ( x ) = δ E 1 δ r 1 ( x ) j = 2 N δ W 1 j δ r 1 ( x )
The first term δ E 1 δ r 1 ( x ) corresponds to the self-force, while the second term j = 2 N δ W 1 j δ r 1 ( x ) corresponds to the interaction force.
Let’s analyze the Interaction Energy and Peach-Koehler Force. For simplicity, consider a system with only two dislocations, so that we only have one interaction term W 12 :
W 12 = S 1 σ i j ( 2 ) ( x ) n i ( 1 ) b j ( 1 ) d S ( x )
where σ i j ( 2 ) ( x ) is the stress field due to the second dislocation loop, n i ( 1 ) is the normal to the surface of the first loop, and b j ( 1 ) is the Burgers vector of the first loop.
The variation of W 12 due to the virtual displacement δ r ( x ) of the first loop is given by:
δ W 12 = δ S 1 σ i j ( 2 ) ( x ) n i ( 1 ) b j ( 1 ) d S ( x )
Using the relation n δ S = δ r × v d L , we can express the variation as:
δ W 12 = L σ i j ( 2 ) ( x ) b j ( 1 ) ϵ i m n δ r m ( x ) v n ( 1 ) ( x ) d L ( x )
This leads to the expression for the force per unit length on the dislocation loop:
f m ( x ) = ϵ i n m σ i j ( 2 ) ( x ) b j ( 1 ) v n ( 1 ) ( x )
In vector notation, this is written as:
f = ( σ · b ) × ξ
This expression is known as the Peach-Koehler force. It describes the force on a dislocation due to an applied stress field, which can originate from other dislocations, external stresses, or any other source. The self-force contribution δ E 1 δ r 1 ( x ) is generally divergent because the self-energy E 1 of the dislocation loop is singular. This divergence is typically handled by introducing a truncation scheme or a non-singular dislocation model, which is discussed in the subsequent sections of the document.

18.6. Non-Singular Dislocation Model

In classical dislocation theory, the stress field and self-energy associated with a dislocation are singular at the dislocation core. This presents difficulties in calculating the self-force on a dislocation. The non-singular dislocation model aims to remove these singularities while maintaining the analytical structure of the original theory. This derivation will rigorously follow the non-singular dislocation model’s development, which involves distributing the dislocation core over a finite region.
The stress field for a dislocation loop in the non-singular model is obtained by convolving the classical (singular) stress field with a spreading function w ( x ) . Consider the classical stress field given by Mura’s formula:
σ α β ( x ) = L C α β k l ϵ l n h C p q m n b m v n ( x ) G k p , q ( x x ) d L ( x )
In the non-singular theory, the stress field is obtained by convolving this expression with a spreading function w ( x ) , which spreads the dislocation core over a finite region:
σ α β ns ( x ) = σ α β ( x ) w ( x )
A commonly used spreading function is:
w ( x ) = 15 a 4 8 π ( | x | 2 + a 2 ) 7 / 2
where a is a small parameter characterizing the spread of the dislocation core. The convolution of G k p , q ( x x ) with w ( x ) modifies the singularity at the core:
G k p , q ns ( x x ) = G k p , q ( x x ) w ( x x ) d x
This convolution results in a non-singular Green’s function:
R w ( x ) = R a R 2 + a 2
Thus, the non-singular stress field becomes:
σ α β ns ( x ) = μ 8 π L i p p R a b m ϵ i m α d x β + b m ϵ i m β d x α + μ 4 ( 1 ν ) L b m ϵ i m k i α β R a δ α β i p p R a d x k
The self-energy of the dislocation loop in the non-singular model is derived similarly to the classical model but using the non-singular Green’s function R a :
E = L L μ 16 π b i b j R a , p p d x i d x j + μ 8 π ( 1 ν ) ϵ i k l ϵ j m n b k b m R a , i j d x l d x n
The interaction energy between two dislocations in the non-singular model is given by:
W 12 = μ 4 π L 1 L 2 ( b 1 × b 2 ) · ( d L 1 × d L 2 ) 2 R a + μ 8 π L 1 L 2 ( b 1 · d L 1 ) ( b 2 · d L 2 ) 2 R a + μ 4 π L 1 L 2 ( b 1 × d L 1 ) · R a · ( b 2 × d L 2 )
With the non-singular stress field, the Peach-Koehler force can now be safely applied without ambiguity:
f = ( σ ns · b ) × ξ

References

  1. W. M. L. Challis, D. Rubin and E. Krempl, Introduction to Continuum Mechanics, (Butterworth-Heinemann, 1999).
  2. L. D. Landau and E. M. Lifshits, Theory of Elasticity, 2nd English ed., (Pergamon Press, Oxford, 1970).
  3. T. Mura, Micromechanics of Defects in Solids, 2nd rev. ed. (Kluwer Academic Publishers, 1991). [CrossRef]
  4. J. R. Barber, Elasticity, 2nd ed. (Kluwer Academic Publishers, 2002).
  5. L. Challis and L. Sheard, The Green of Green Functions, Physics Today, Dec. 2003.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2025 MDPI (Basel, Switzerland) unless otherwise stated