Two-dimensional small-strain elasticity. Crystal axes: x = [10], y = [01]. Shape and eigenstrain are prescribed; these experiments calculate energies and orientations, not a complete evolving microstructure. All numerical inputs use consistent dimensionless units.
Each module has its own inputs. Changing or resetting one module does not change another. Try these values recalculates that module; Reset this module restores its initial values.
Elastic constants and mechanical stability
Elastic constants relate elastic strain to stress: \(\sigma_{ij}=C_{ijkl}e_{kl}\). Here \(\sigma_{ij}\) is stress in pascals, \(e_{kl}\) is dimensionless elastic strain, and \(C_{ijkl}\) is stiffness in pascals. The indices identify Cartesian components. The elastic energy per volume is \(w=\tfrac12 C_{ijkl}e_{ij}e_{kl}\), in joules per cubic metre. The numerical examples below use a common scaled energy-density unit rather than a specified material.
Positive definite means positive elastic energy for every nonzero elastic strain. A symmetric stiffness matrix must have strictly positive eigenvalues. A zero eigenvalue is a marginal mode; a negative eigenvalue permits a strain that lowers the reference-state energy. These are small-strain elastic stability conditions for an unstressed crystal, not a complete test of lattice dynamical stability or stability under finite applied stress.
Cubic crystals: three independent stiffness constants
The first condition resists a volume-preserving tetragonal distortion. The second gives a positive bulk modulus \(K=(C_{11}+2C_{12})/3\), which resists uniform compression. The third resists shear in the principal planes. The normal-strain block has eigenvalues \(C_{11}+2C_{12}\) and \(C_{11}-C_{12}\), the latter occurring twice. The three engineering-shear entries are \(C_{44}\). Stiffness eigenvalues and \(K\) use the same units as the stiffnesses.
The dimensionless Zener ratio measures cubic elastic anisotropy. The equality gives an isotropic elastic tensor within the cubic law. A cubic crystal is not automatically isotropic.
Convert between stiffnesses and Schmidt–Gross parameters
Following I. Schmidt and D. Gross, Sandeep’s thesis uses the average in-plane shear modulus \(G\), the effective Poisson parameter \(\nu\), and \(A_Z\). Here the same average shear parameter is labeled \(\mu=G\). Specifically, \(\mu\) is the angular average of the rotated \(C_{1212}\) in the (001) plane, and \(\nu=\Lambda/[2(\Lambda+\mu)]\), where \(\Lambda\) is the angular average of the rotated \(C_{1122}\). These are averaged parameters, not the shear modulus and Poisson ratio for every direction in an anisotropic crystal.
At \(A_Z=1\), the mapping becomes the usual isotropic law: \(C_{44}=\mu\), \(C_{12}=2\mu\nu/(1-2\nu)\), and \(C_{11}=2\mu(1-\nu)/(1-2\nu)\). For the anisotropic mapping with \(\mu>0\), \(A_Z>0\), and \(\nu<1/2\), cubic stability additionally requires \(1+3A_Z+4\nu>0\). Always check the converted stiffnesses. Module IV retains the requested fixed \(\nu=1/3\); the other conversion panels permit changing \(\nu\).
Change this calculator’s constants directly or use the equivalent parameters. It does not change another module. Try an unstable case to see which deformation loses its restoring energy.
The graph uses strain amplitude \(q\), which is dimensionless. Uniform dilation has \(e_{xx}=e_{yy}=e_{zz}=q\); tetragonal distortion has \(e_{xx}=q,e_{yy}=-q,e_{zz}=0\); simple shear has engineering shear \(\gamma_{xy}=q\), so tensor shear is \(e_{xy}=q/2\). All other components vanish in each example. The plotted energies are respectively \(\tfrac32(C_{11}+2C_{12})q^2\), \((C_{11}-C_{12})q^2\), and \(\tfrac12C_{44}q^2\).
Other crystal symmetries
The following are three-dimensional stiffness conditions. The later widgets use the in-plane cubic sub-block; they are not hexagonal or tetragonal solvers.
| Crystal symmetry | Independent constants | Stability conditions |
|---|---|---|
| Hexagonal | 5 | \(C_{11}>|C_{12}|\), \(C_{44}>0\), and \((C_{11}+C_{12})C_{33}>2C_{13}^2\). Here \(C_{66}=(C_{11}-C_{12})/2>0\) follows from the first condition. |
| Tetragonal I (4/mmm) | 6 | \(C_{11}>|C_{12}|\), \((C_{11}+C_{12})C_{33}>2C_{13}^2\), \(C_{44}>0\), \(C_{66}>0\). Here \(C_{66}\) is independent. Lower-symmetry tetragonal classes require an additional condition involving \(C_{16}\). |
| Orthorhombic | 9 | \(C_{44},C_{55},C_{66}>0\); for the normal block, \(C_{11}>0\), \(C_{11}C_{22}-C_{12}^2>0\), and \(\det\begin{pmatrix}C_{11}&C_{12}&C_{13}\\C_{12}&C_{22}&C_{23}\\C_{13}&C_{23}&C_{33}\end{pmatrix}>0\). Positive diagonal entries alone are insufficient. |
References: F. Mouhat and F.-X. Coudert, Necessary and sufficient elastic stability conditions in various crystal systems, Physical Review B 90, 224104 (2014); I. Schmidt and D. Gross, The equilibrium shape of an elastically inhomogeneous inclusion, Journal of the Mechanics and Physics of Solids 45, 1521–1549 (1997); Sandeep Sugathan’s thesis, Chapter 3, Eqs. 3.23–3.29, printed pp. 29–30. For tensor notation, see J. F. Nye, Physical Properties of Crystals: Their Representation by Tensors and Matrices.
Misfit and the habit line.
Eigenstrain is the deformation a phase would adopt if unconstrained. Elastic strain is the part left after subtracting that eigenstrain from the compatible total strain.
Construct Mohr’s circle of strain.
Axes below are strain divided by \(\epsilon_\eta\). The vertical coordinate is tensor shear \(\epsilon^0_{nq}\), not engineering shear \(\gamma^0_{nq}=2\epsilon^0_{nq}\). While θ follows the zero-extension line, P stays on a green intersection as t changes. Type a measurement angle above to release it and move the point around the circle.
Why did this happen?
Physics: eigenstrain, Eshelby’s inclusion and the habit line
Eigenstrain
Eigenstrain is the stress-free transformation strain \(\epsilon^0\). When a region transforms with no external constraint, its stress stays zero:
where \(C\) is the elastic stiffness tensor and \(\epsilon\) the total strain. In the free region the total strain equals the transformation strain, \(\epsilon=\epsilon^0\). Examples are the lattice-misfit dilatation of a precipitate, \(\epsilon^0_{ij}=\delta_{ij}(a_p-a_m)/a_m\), thermal expansion, and the spontaneous distortion of a lattice during a phase change. In this module \(\epsilon^0=\epsilon_\eta\,\mathrm{diag}(1,t)\) in the crystal axes: principal strains \(\epsilon_1=\epsilon_\eta\) along [10] and \(\epsilon_2=t\,\epsilon_\eta\) along [01].
Eshelby’s four steps
To find the stress when the transforming region is constrained by a surrounding elastic matrix, Eshelby (1957) used a thought experiment:
- Cut the region Ω (the inclusion) out of the matrix and let it undergo its stress-free transformation \(\epsilon^0\).
- Apply surface tractions to force it back into its original hole. This imposes elastic strain \(-\epsilon^0\), i.e. stress \(-C:\epsilon^0\).
- Weld the interface back together.
- Release the tractions. The inclusion and the matrix relax together to mechanical equilibrium.
The relaxation leaves a constrained total strain \(\epsilon^c\) inside the inclusion. The elastic strain is the difference between the constrained and the stress-free geometry, and only it stores energy:
For an ellipsoidal inclusion (an ellipse in 2-D) in a homogeneous infinite matrix, Eshelby proved that \(\epsilon^c\) is uniform and linear in the eigenstrain, \(\epsilon^c=S:\epsilon^0\). The Eshelby tensor S depends on the inclusion’s shape and orientation and on the elastic constants. This is where orientation dependence comes from: the same \(\epsilon^0\) costs different energy in differently oriented inclusions.
Thin-plate limit: a negotiation at the interface
Eshelby’s theory is a negotiation of space. The transformed region wants a new shape; the surrounding matrix resists. The theory decides who yields, and where stress is left trapped. For a very thin plate (a martensite lath, a plate-shaped precipitate) the whole negotiation happens at the flat interface. The plate is long compared with its thickness, so the matrix, which extends far beyond the plate and is unstrained away from its ends, dictates what happens along the plate. Across its thickness, the plate is free to relax.
Why jump conditions appear at an inclusion
The eigenstrain switches on abruptly at the interface: it is zero in the matrix and \(\epsilon^0\) in the plate. Suppose every strain and stress component were continuous across the interface. Hooke’s law, \(\epsilon^0=\epsilon-C^{-1}:\sigma\), would then make the eigenstrain continuous as well. It is not. So some components of strain or stress must jump. Two physical laws decide which ones:
- No tearing, no overlapping. The interface is welded, so the displacement is continuous. This constrains the strain.
- Newton’s third law. The force that the plate exerts on the matrix across the interface equals and opposes the force the matrix exerts on the plate. This constrains the stress.
The interface handshake
1. The matrix rules the length (Hadamard kinematics). Plate and matrix cannot separate, so their displacements match everywhere along the interface direction \(\mathbf s\). Their derivatives along \(\mathbf s\) then match too, and the plate’s total strain along the interface equals the matrix value: \(\epsilon_{ss}=0\). The plate is still free to thicken, thin or shear across its thickness (direction \(\mathbf m\)). That freedom is a jump in the displacement gradient of the form \(\mathbf b\otimes\mathbf m\), where \(\mathbf b\) is not yet known.
2. The forces rule the thickness (traction continuity). The forces pushing, pulling and sliding across the interface must balance. So the traction vector \(\sigma\mathbf m\), which contains the normal stress \(\sigma_{mm}\) and the shear stress \(\sigma_{sm}\), is continuous. No balance law acts on the stress that stretches the plate along its length, \(\sigma_{ss}\): it may differ completely from the matrix value.
Across the welded interface, each component is fixed by exactly one law. When the strain component is fixed, the matching stress component may jump, and the other way round:
| component | strain \(\epsilon\) | stress \(\sigma\) |
|---|---|---|
| ss (along the plate) | continuous: Hadamard (matches the matrix) | free to jump: remnant stress |
| sm (shear) | free to jump: set by \(\mathbf b\) | continuous: traction |
| mm (across the plate) | free to jump: set by \(\mathbf b\) | continuous: traction |
Counting. In 2-D a symmetric strain has 3 components. Hadamard fixes 1 of them (\(\epsilon_{ss}\)); the other 2 are carried by the 2 unknown components of \(\mathbf b\). Traction continuity supplies 2 equations (\(\sigma_{sm},\sigma_{mm}\)), exactly enough to find \(\mathbf b\). In 3-D: 6 strain components, 3 fixed in-plane by Hadamard, 3 carried by \(\mathbf b\), and 3 traction equations to find them.
Duality and the remnant stress
The unknown deformation \(\mathbf b\) and the interface traction \(\sigma\mathbf m\) are work conjugates: a movement and the force that resists it. Changing \(\mathbf b\) by \(\delta\mathbf b\) changes the plate’s energy density by \(\delta w=(\sigma\mathbf m)\cdot\delta\mathbf b\). So when the plate chooses \(\mathbf b\) to minimize its energy, the traction on its faces is driven to zero: \(\sigma_{mm}=\sigma_{sm}=0\). The plate relaxes completely across its thickness.
What is left is the stress along the plate. The plate is still held to zero length change along \(\mathbf s\), although it wanted to change by \(\epsilon^0_{ss}\). It carries the elastic strain \(e_{ss}=-\epsilon^0_{ss}\) and, because its faces are free, a single uniaxial stress
where \(E_s\) is the (plane-strain) stiffness of the crystal along \(\mathbf s\). This is the whole elastic penalty of a thin plate: the remnant stress it pays for staying compatible with the matrix. It vanishes only when \(\epsilon^0_{ss}=0\), on the habit line. In the worked case of Step 7(vii) the plate wanted to shrink by 0.005 along \(\mathbf s\), was held at its original length, and carries a tension of 0.0133. Figure 3 shows every component across that interface.
Names
- \(\mathbf n,\ \mathbf q\): an arbitrary pair of orthogonal unit vectors defining a rotated frame; \(\mathbf n\) is at angle θ from [10].
- \(\mathbf s\): the invariant (habit) line direction, the value of \(\mathbf n\) where \(\epsilon^0_{nn}=0\). Green line on the plot.
- \(\mathbf m\): the habit normal, orthogonal to \(\mathbf s\); the value of \(\mathbf q\) at that θ. Orange arrow on the plot.
“Habit angle” is used in the literature for both the line and the normal, which differ by 90°. This page always says which.
Limits
The ideal plate theory assumes an infinitely long plate, a sharp interface, a uniform transformation strain, linear small-strain elasticity and a homogeneous stiffness. Real plates have ends, finite thickness, interfacial energy and diffuse interfaces, and plasticity can destroy coherency. In this laboratory, Module II replaces the ideal plate by a finite elliptical inclusion with Eshelby’s full solution, Module III adds elastic anisotropy, Module IV two phases with different eigenstrains, and Module V plates stiffer or softer than the matrix.
Planned extensions. Future modules will relax the remaining assumptions: inclusions stiffer or softer than the matrix through Eshelby’s equivalent-inclusion method, interfacial energy and phase-field gradients for diffuse interfaces and evolving shapes, dislocation arrays for the loss of coherency, and the 3-D invariant-plane analysis of the phenomenological theory of martensite crystallography.
Mathematics: derivation, step by step
Small strain; crystal axes x = [10], y = [01]; repeated indices summed. Vectors are columns.
Step 1. Rotate \(\epsilon^0\) into the \((\mathbf n,\mathbf q)\) frame.
In its principal frame the eigenstrain is diagonal, \(\epsilon^0=\begin{pmatrix}\epsilon_1&0\\0&\epsilon_2\end{pmatrix}\). A frame at angle θ has basis vectors and rotation matrix (rows = new basis vectors)
Each new component is the tensor sandwiched between two basis vectors, e.g. \(\epsilon^0_{nq}=\mathbf n^{\mathsf T}\epsilon^0\,\mathbf q=n_i\,\epsilon^0_{ij}\,q_j\). All four at once are \(\epsilon'=R\,\epsilon^0R^{\mathsf T}\). Multiply in full:
So the entries are
Check: the trace \(\epsilon^0_{nn}+\epsilon^0_{qq}=\epsilon_1+\epsilon_2\) and the determinant \(\epsilon^0_{nn}\epsilon^0_{qq}-(\epsilon^0_{nq})^2=\epsilon_1\epsilon_2\) do not change with θ. For this module put \(\epsilon_1=\epsilon_\eta\), \(\epsilon_2=t\,\epsilon_\eta\): \(\epsilon^0_{nn}=\epsilon_\eta(\cos^2\theta+t\sin^2\theta)\), the curve in the “Normal eigenstrain along n(θ)” plot. At a habit root, \(\mathbf n\) becomes \(\mathbf s\) and \(\mathbf q\) becomes \(\mathbf m\).
Step 2. Eshelby’s thin plate.
(a) Hadamard compatibility. Let the interface be the plane with normal \(\mathbf m\). The displacement \(\mathbf u\) is continuous across it, so its derivative along any direction \(\mathbf s\) lying in the interface (\(\mathbf s\cdot\mathbf m=0\)) is the same on both sides. The jump in the gradient therefore kills every in-plane vector:
The matrix far from the plate ends is unstrained, so inside the plate \(\nabla\mathbf u=\mathbf b\otimes\mathbf m\). Here \(\mathbf b\) is an unknown vector (a displacement per unit distance across the plate). In index notation,
where i labels the displacement component and j the direction of the derivative. In 3-D, with \(\mathbf u=(u_1,u_2,u_3)^{\mathsf T}\), \(\mathbf b=(b_1,b_2,b_3)^{\mathsf T}\) and \(\mathbf m=(m_1,m_2,m_3)^{\mathsf T}\), the outer product is the column \(\mathbf b\) times the row \(\mathbf m^{\mathsf T}\):
Every entry of the gradient is fixed by one component of \(\mathbf b\) and one of \(\mathbf m\); for example the shear gradient \(\partial u_1/\partial x_2=b_1m_2\). The matrix has rank one: every row is a multiple of \(\mathbf m^{\mathsf T}\), so \(\nabla\mathbf u\,\mathbf s=\mathbf b\,(\mathbf m\cdot\mathbf s)=\mathbf 0\) for any in-plane \(\mathbf s\). Nothing changes along the interface. In this 2-D laboratory \(m_3=b_3=0\) and the upper-left block remains:
(b) Energy. The elastic strain is \(e=\operatorname{sym}(\mathbf b\otimes\mathbf m)-\epsilon^0\). With \(\sigma^0=C:\epsilon^0\) and the acoustic tensor \(Q_{ik}(\mathbf m)=C_{ijkl}m_jm_l\), and using the symmetries of C,
(c) Traction balance. Minimising over \(\mathbf b\), \(\partial w/\partial\mathbf b=Q(\mathbf m)\mathbf b-\sigma^0\mathbf m=\mathbf 0\). This is exactly zero traction on the faces, \(\sigma\mathbf m=(C:e)\,\mathbf m=\mathbf 0\). Substituting \(\mathbf b=Q(\mathbf m)^{-1}\sigma^0\mathbf m\) back:
The first term is the cost of forcing the plate fully back (Eshelby’s step 2). The second is what it recovers by relaxing across its thickness (step 4). Their difference is the elastic energy kernel \(B(\mathbf m)\) of Module III and of Appendix A. In Eshelby’s language this is the thin-plate limit, with \(S:\epsilon^0=\operatorname{sym}(\mathbf b\otimes\mathbf m)\).
(d) The two jump conditions, derived. Let \([\![f]\!]=f^{\rm plate}-f^{\rm matrix}\) denote a jump across the interface. Kinematic: \(\mathbf u\) is continuous, so for any in-plane \(\mathbf s\) the directional derivative is continuous, \([\![\nabla\mathbf u]\!]\mathbf s=\mathbf 0\), which gives \([\![\nabla\mathbf u]\!]=\mathbf b\otimes\mathbf m\) as in (a). Static: apply \(\int_{\partial P}\sigma\,\mathbf n\,dA=\mathbf 0\) to a pillbox P of face area A straddling the interface. As its height goes to zero the side faces vanish and only the two flat faces remain:
With a homogeneous C, \([\![\sigma]\!]=C:\big(\operatorname{sym}(\mathbf b\otimes\mathbf m)-[\![\epsilon^0]\!]\big)\), so the static condition reads \(Q(\mathbf m)\,\mathbf b=(C:[\![\epsilon^0]\!])\,\mathbf m\). These hold at every point of every coherent interface. For the thin plate the matrix side is unstrained and \([\![\epsilon^0]\!]=\epsilon^0\), which returns \(Q(\mathbf m)\mathbf b=\sigma^0\mathbf m\) of (c): the thin-plate solution is the pair of jump conditions.
(e) Duality. Vary \(\mathbf b\): \(\delta w=\sigma:\delta e=\sigma:\operatorname{sym}(\delta\mathbf b\otimes\mathbf m) =\delta b_i\,\sigma_{ij}m_j=(\sigma\mathbf m)\cdot\delta\mathbf b\). The traction is the work conjugate of \(\mathbf b\), so \(\partial w/\partial\mathbf b=\mathbf 0\) and \(\sigma\mathbf m=\mathbf 0\) are the same statement.
(f) The remnant stress and a one-line energy. With \(\sigma\mathbf m=\mathbf 0\), the only stress left in the \((\mathbf s,\mathbf m)\) frame is \(\sigma_{ss}\): \(\sigma=\sigma_{ss}\,\mathbf s\otimes\mathbf s\). Hadamard fixes \(\epsilon_{ss}=0\), so \(e_{ss}=-\epsilon^0_{ss}\). Using the compliance \(S=C^{-1}\) (plane strain), \(e_{ss}=S_{ssss}\,\sigma_{ss}\), hence
This is exact for any homogeneous stiffness, isotropic or anisotropic, and it makes Step 3 immediate. For the isotropic law of this page \(E_s=(C_{11}^2-C_{12}^2)/C_{11}=8/3\) in every direction, so \(B=\tfrac83(0.005)^2=6.67\times10^{-5}\) for a plate along [01] and \(\tfrac83(0.01)^2=2.67\times10^{-4}\) along [10], matching Step 7 and the Module III curve. In an anisotropic crystal \(E_s\) depends on the direction, which is how elastic anisotropy enters in Module III.
(g) Any coherent interface: the complete set of jump conditions. Steps (a)–(f) used a thin plate, but nothing in the jump conditions depends on the plate. They hold at every point of any perfectly bonded interface between an inclusion with uniform eigenstrain \(\epsilon^0_{ij}\) (often written \(\epsilon^*_{ij}\)) and a matrix with the same stiffness \(C_{ijkl}\). Let \(n_i\) be the outward unit normal of the inclusion and \([\![A]\!]=A^{\rm in}-A^{\rm out}\), the convention used on this page.
1. Displacement gradient. \(u_i\) is continuous, so only the derivative across the interface can jump (Hadamard):
2. Strain. The total strain is the symmetric part of the gradient, and the eigenstrain jumps by its full value:
3. Stress. With the same \(C_{ijkl}\) on both sides and its minor symmetry,
4. Solve for the amplitude. Substitute the stress jump into traction continuity and use the acoustic (Christoffel) tensor \(K_{ik}=C_{ijkl}n_jn_l\):
So the jumps are fixed locally by the normal, the stiffness and the eigenstrain alone, whatever the shape of the inclusion. The shape only decides the field on one side (Module II computes it for an ellipse). Texts that define the jump as outside minus inside write \(\lambda_i=-g_i\), and every jump below changes sign.
| field | across a coherent interface | jump \([\![A]\!]=A^{\rm in}-A^{\rm out}\) |
|---|---|---|
| displacement | continuous | \([\![u_i]\!]=0\) |
| traction | continuous | \([\![\sigma_{ij}]\!]n_j=0\) |
| displacement gradient | jumps | \([\![u_{i,j}]\!]=g_in_j\) |
| total strain | jumps (in-plane part continuous) | \([\![\epsilon_{ij}]\!]=\tfrac12(g_in_j+g_jn_i)\) |
| eigenstrain | jumps (prescribed) | \([\![\epsilon^0_{ij}]\!]=\epsilon^0_{ij}\) |
| elastic strain | jumps | \([\![e_{ij}]\!]=\tfrac12(g_in_j+g_jn_i)-\epsilon^0_{ij}\) |
| stress | jumps (in-plane part only) | \([\![\sigma_{ij}]\!]=C_{ijkl}(g_kn_l-\epsilon^0_{kl})\) |
| amplitude | solved from traction | \(g_i=K^{-1}_{ik}C_{kjml}\epsilon^0_{ml}n_j\) |
These are Hill’s discontinuity relations for a coherent interface (R. Hill, “Discontinuity relations in mechanics of solids”, Progress in Solid Mechanics 2, 247–276, 1961; see also Mura, 1987). For the thin plate, the matrix side is unstrained, so the plate’s own field equals the jump and \(\mathbf b=\mathbf g\).
Step 3. \(B(\mathbf m)=0\iff\epsilon^0_{ss}=0\).
(⇒) If \(B(\mathbf m)=0\) then \(w=0\), and because C is positive definite, \(e=0\): the total strain equals the transformation strain, \(\operatorname{sym}(\mathbf b\otimes\mathbf m)=\epsilon^0\). Project along the interface direction \(\mathbf s\) (\(\mathbf s\cdot\mathbf m=0\)):
(⇐) If \(\epsilon^0_{ss}=0\), write both sides in the \((\mathbf s,\mathbf m)\) frame, where \(\mathbf m=(0,1)^{\mathsf T}\):
They match with \(a_s=2\epsilon^0_{sm}\), \(a_m=\epsilon^0_{mm}\). Then \(\mathbf b=\mathbf a\) gives \(e=0\) and \(B(\mathbf m)=0\). The condition is kinematic: it holds for every stiffness C. The elastic constants only decide how costly the other orientations are.
Step 4. Solve for θ and the plate kinematics.
With \(\epsilon_1=\epsilon_\eta\), \(\epsilon_2=t\,\epsilon_\eta\) and \(\cos\theta\ne0\):
A real root requires \(t\le0\) (\(t=0\) is the limit θ = 90°). On the root \(\cos^2\theta=-t/(1-t)\), \(\sin^2\theta=1/(1-t)\) and \(\sin\theta\cos\theta=\pm\sqrt{-t}/(1-t)\). Substituting into Step 1:
The upper sign is for \(\theta_1\), the lower for \(\theta_2\). The rotation is the antisymmetric part of \(\mathbf a\otimes\mathbf m\): \(\mathbf a\otimes\mathbf m=\epsilon^0+\epsilon^0_{sm}(\mathbf s\otimes\mathbf m-\mathbf m\otimes\mathbf s)\). In units of \(\epsilon_\eta\):
| \(t=\epsilon_2/\epsilon_1\) | line θ₁ (θ₂) | thickening \(\epsilon^0_{mm}\) | shear \(2\epsilon^0_{sm}\) | rotation ω | ω for \(\epsilon_\eta=0.01\) |
|---|---|---|---|---|---|
| −1 | 45° (135°) | 0 | ∓2.00 | ∓1.00 | 0.57° |
| −0.5 | 54.74° (125.26°) | 0.50 | ∓1.41 | ∓0.71 | 0.41° |
| −0.25 | 63.43° (116.57°) | 0.75 | ∓1.00 | ∓0.50 | 0.29° |
| 0 | 90° | 1.00 | 0 | 0 | 0 |
Step 5. Where \(\alpha=\tan^{-1}\sqrt{-t}\) comes from.
The habit condition can be written for the line or for the normal. Both routes give the same plate.
Route A: start from the line. Step 4 gave the line angle through \(\tan\theta_1=1/\sqrt{-t}\). The normal is the line turned by 90°:
The acute angle α between the normal and [10] satisfies \(\tan\alpha=|m_y/m_x|=\cos\theta_1/\sin\theta_1=1/\tan\theta_1\). So
Route B: start from the normal. Call the normal angle α, so \(\mathbf m=(\cos\alpha,\sin\alpha)\). The tangent is \(\mathbf s=(-\sin\alpha,\cos\alpha)\). Put \(\mathbf s\) into the stretch formula of Step 1 (the same sandwich \(\mathbf s\cdot\epsilon^0\mathbf s\)):
Set it to zero and divide by \(\cos^2\alpha\):
Route B is the form you meet when the energy is written in terms of the normal, as in Module III and in the thesis. The same distinction between the interface line and its normal applies to habit planes in 3-D, a planned extension of this laboratory. The line and the normal pair up as follows (angles from [10]; a direction and its opposite are the same line):
| line θ | normal (direction) | acute normal angle α |
|---|---|---|
| \(\theta_1=\tan^{-1}(1/\sqrt{-t})\) | \(\theta_1+90^\circ\equiv-\alpha\) | \(\tan^{-1}\sqrt{-t}\) |
| \(\theta_2=180^\circ-\theta_1\) | \(\theta_2-90^\circ=+\alpha\) | \(\tan^{-1}\sqrt{-t}\) |
For t = −0.5: lines at 54.74° and 125.26°; normals at 144.74° and 35.26°. Module III lists its zero minima of \(B(\mathbf n)\) at these two normal directions (sampled at 1°, so near 145° and 35°).
Step 6. Mohr’s circle.
Apply \(\cos^2\theta=\tfrac12(1+\cos2\theta)\), \(\sin^2\theta=\tfrac12(1-\cos2\theta)\) and \(2\sin\theta\cos\theta=\sin2\theta\) to the Step 1 entries:
As θ grows, \((\epsilon^0_{nn},\epsilon^0_{nq})\) moves round a circle with center \(\bar c\) and radius R, turning by −2θ. The vertical axis is tensor shear; engineering shear is twice it. The circle meets \(\epsilon^0_{nn}=0\) only if \(R\ge\bar c\), i.e. \(t\le0\), at
the green dots: the same shear \(\epsilon^0_{sm}\) as in Step 4. (The plot labels the center c; it is \(\bar c\) here to avoid clashing with \(\cos\theta\).)
Step 7. One case solved completely.
Take the page defaults: \(t=-0.5\), \(\epsilon_\eta=0.01\), and the isotropic stiffness of Module III, \(C_{11}=3\), \(C_{12}=1\), \(C_{44}=1\) (dimensionless). We find the habit line, then compute the plate energy twice: at the habit normal and at a poor orientation.
(i) Eigenstrain.
(ii) Habit line. \(\tan^2\theta=-1/t=2\), so \(\tan\theta_1=\sqrt2\), \(\theta_1=54.7356^\circ\). Then \(\cos^2\theta_1=1/3\), \(\sin^2\theta_1=2/3\):
(iii) Rotate the eigenstrain (Step 1 with numbers).
Entry by entry: \(\epsilon^0_{ss}=0.01\cdot\tfrac13-0.005\cdot\tfrac23=0\); \(\epsilon^0_{mm}=0.01\cdot\tfrac23-0.005\cdot\tfrac13=0.005\); \(\epsilon^0_{sm}=-0.015\times0.57735\times0.81650=-0.0070711\). Trace 0.005 and determinant \(-5\times10^{-5}\) agree with the original diagonal matrix. The Mohr readout shows the last entry as \(-0.7071\,\epsilon_\eta\).
(iv) The frame vectors and the compatible vector a.
Check that it reproduces the eigenstrain in crystal axes:
The skew part is a rotation by 0.0070711 rad (0.41°). The eigenstrain is exactly a compatible plate deformation plus this rigid rotation.
(v) Stiffness in matrix form. For the in-plane cubic law with tensor shear, \(\sigma_{xx}=C_{11}\epsilon_{xx}+C_{12}\epsilon_{yy}\), \(\sigma_{yy}=C_{12}\epsilon_{xx}+C_{11}\epsilon_{yy}\), \(\sigma_{xy}=2C_{44}\epsilon_{xy}\). The eigenstress is
The acoustic tensor \(Q_{ik}(\mathbf m)=C_{ijkl}m_jm_l\) for this law is
(vi) Energy at the habit normal. With \(m_x^2=2/3\), \(m_y^2=1/3\), \(m_xm_y=-0.47140\):
The plate relaxes all of its eigenstrain: \(\mathbf b=\mathbf a\), \(e=0\), no stress, no energy.
(vii) Energy at a poor orientation. Put the plate along [01], so \(\mathbf m=(1,0)\):
Read this physically. \(\sigma_{xx}=0\): the faces (normal x) carry no traction, as Step 2 demands. \(e_{yy}=0.005\): along the plate the particle wanted to shrink by 0.005 but the matrix holds its length, so it is stretched elastically and carries tension 0.0133. The energy density is \(\tfrac12 e:\sigma=\tfrac12(0.005)(0.013333)=3.33\times10^{-5} =\tfrac12B\). Module III’s curve gives \(B=6.67\times10^{-5}\) at normal 0° and \(2.67\times10^{-4}\) at 90°, and zero at the habit normals.
References. J. D. Eshelby, Proc. R. Soc. Lond. A 241, 376–396 (1957). T. Mura, Micromechanics of Defects in Solids, 2nd ed., Kluwer (1987). A. G. Khachaturyan, Theory of Structural Transformations in Solids, Wiley (1983).
Try it yourself
- Press “t = −1”. Check that the lines are at 45° and 135°, and that the Mohr circle is centered on the origin.
- Drag the t slider slowly from −1 toward 0. Watch θ follow the line and turn toward [01], and the green dots move together.
- Push t above 0. The circle leaves the vertical axis and the habit line disappears.
- Type your own θ (for example 30°). The angle is now free. Move it and watch P travel round the circle at twice the rate of the plate.
- Set t = −0.5, then open Module III and compare its minimum normal with α = 35.26° from Step 5.
Computation
The page evaluates \(\epsilon^0_{nn}(\theta)\) directly and solves its zero condition analytically; no remembered angle is inserted. Positive t returns no roots and t = 0 is handled as the limit θ = 90°. When θ follows t, the calculation uses the exact root, and the input box shows it rounded to four decimals. The plate sketch rotates a prescribed shape; it is not a kinetic simulation.
Convention note: my thesis Eq. (3.2), p. 20, uses diag(1, t). Eq. (4.10), p. 48, gives \(\tan^{-1}\sqrt{-t}\) as a line angle from x. Direct projection under Eq. (3.2) gives that value for the interface normal, and its complement for the line. This laboratory labels the line and the normal separately.
Eshelby’s inclusion: what jumps at the interface.
Put the eigenstrain of Module I into an elliptical inclusion inside a homogeneous matrix. Eshelby’s solution gives the uniform strain inside. Pick any point on the interface to see which strain and stress components are continuous there, which jump, and by how much. Then flatten the ellipse into a thin plate and rotate it.
Jump conditions at the chosen point, (s, m) frame. Green rows are continuous; orange rows jump.
See the whole field.
The map uses the selectable periodic FFT grid. The profiles below instead evaluate the uniform interior and the spatially varying exterior of an isolated ellipse in an infinite matrix separately. Both branches end at their analytical one-sided interface values. Dashed vertical connectors mark actual jumps; no interpolation or smoothing crosses the interface. Periodic map values can differ from these infinite-matrix profiles.
Why did this happen?
Physics: the inclusion, its interface and its limits
Eshelby’s uniform interior
Module I described Eshelby’s cut, force back, weld and release experiment. Its remarkable result is that for an ellipse (an ellipsoid in 3-D) in a homogeneous infinite matrix, the constrained strain inside is uniform, whatever the aspect ratio and orientation. The inclusion also turns rigidly by a small angle ω. The panel on the left shows three outlines: the hole, the shape the inclusion would take if it were free, and the shape it actually takes.
The interface is local
At every point of the interface the two jump conditions of Module I hold, with the local normal \(\mathbf m\) and tangent \(\mathbf s\): welding keeps the displacement continuous, so \(\epsilon_{ss}\) is continuous; force balance keeps the traction \(\sigma\mathbf m\) continuous, so \(\sigma_{sm}\) and \(\sigma_{mm}\) are continuous. The other components, \(\epsilon_{sm}\), \(\epsilon_{mm}\) and \(\sigma_{ss}\), jump. The size of the jump depends only on \(\mathbf m\), the eigenstrain and the stiffness. It does not depend on the inclusion’s shape. The shape only sets the uniform field inside; the field just outside is that interior field minus the jump.
From circle to thin plate
For a circle the matrix pushes back from every direction, and the interior keeps a mix of all components. As the ellipse flattens, its long faces dominate. The matrix beside those faces becomes unstrained, so the interior field becomes exactly the jump across the faces: the thin plate of Module I. Its energy per unit area approaches \(\tfrac12B(\mathbf m)\). If the plate lies along a habit line, that limit is zero; otherwise a remnant stress \(\sigma_{ss}\) stays trapped along the plate.
Orientation
Rotating the ellipse changes which direction carries the long interface. For a thin plate the energy follows \(\tfrac12B(\mathbf m)\) and is lowest at the habit orientation. A rounder inclusion averages over many normals, so its energy depends less on orientation.
Mathematics: interior field, jumps and limits
Step 1. The interior displacement gradient.
Let the ellipse have semi-axes a (along the angle φ) and b. For a direction n write \(Q(\mathbf n)_{ik}=C_{ijkl}n_jn_l\), \(\sigma^0=C:\epsilon^0\) and
The Fourier solution of equilibrium gives the total strain of any inclusion as \(\hat\epsilon(\mathbf k)=\operatorname{sym}\big(\mathbf g(\hat{\mathbf k})\otimes\hat{\mathbf k}\big)\,\hat\theta(\mathbf k)\), where \(\theta\) is 1 inside and 0 outside. For an ellipse the transform of θ depends on \(\mathbf k\) only through one elliptic radius, and the inverse transform at any interior point reduces to an average over directions (Mura, §11):
From D follow the interior strain, rotation, stress and energy:
w is Eshelby’s elastic energy per unit inclusion area (matrix energy included). The integral reproduces Mura’s closed-form isotropic Eshelby tensor to machine precision, and it works unchanged for the anisotropic cubic law.
Step 2. The jump at a boundary point.
At a point with outward normal \(\mathbf m\), Hadamard’s condition gives \([\![\epsilon]\!]=\operatorname{sym}(\mathbf g\otimes\mathbf m)\) for some \(\mathbf g\), and traction continuity \([\![\sigma]\!]\mathbf m=\mathbf 0\) with \([\![\sigma]\!]=C:([\![\epsilon]\!]-\epsilon^0)\) gives \(Q(\mathbf m)\mathbf g=\sigma^0\mathbf m\). So
The jump is the thin-plate field of Module I for the normal \(\mathbf m\): it is fixed locally, independent of shape. The table rotates each tensor into the \((\mathbf s,\mathbf m)\) frame with the sandwich rule of Module I, Step 1.
Step 3. The two limits.
Thin plate (\(b/a\to0\)): \(\mathbf n\) points along \(\pm\mathbf m_0\), the normal of the long faces, for almost every \(\vartheta\). Then \(D\to\mathbf g(\mathbf m_0)\otimes\mathbf m_0\), the field outside the faces vanishes, and \(w\to\tfrac12B(\mathbf m_0)\). Circle (\(b=a\)): \(\mathbf n\) is spread evenly; for isotropic stiffness and a dilatation \(\epsilon^0=\epsilon\,I\) (plane strain) the interior strain is \(\epsilon^c=\epsilon\,I/[2(1-\nu)]\).
Step 4. The field maps.
The maps solve the same equation on a periodic 256 × 256 cell: \(\hat\epsilon(\mathbf k)=\operatorname{sym}(\mathbf g(\hat{\mathbf k})\otimes\hat{\mathbf k})\hat\theta(\mathbf k)\) for \(\mathbf k\ne0\), mean total strain zero, and \(\sigma=C:(\epsilon-\theta\epsilon^0)\). The profiles plot this numerical field along \(\mathbf m\) through the chosen point, with the exact inside and outside values of Step 2 dotted.
Try it yourself
- Press “Circle”, set t = 1 (pure dilatation) and click round the boundary. The jump table changes with m, but the interior summary does not: the interior is uniform.
- Set t = −0.5, press “Thin plate” and “Align with habit line”. The energy drops to nearly zero and \(\sigma_{ss}\) on the long face vanishes. Now press “Worst orientation”.
- Drag the aspect-ratio slider from 1 to 0.01 and watch w approach the dashed thin-plate limit.
- Choose the map component \(\sigma_{xy}\) and move the interface point. In the profiles, the continuous components meet at the same one-sided interface value while the discontinuous components have distinct one-sided limits.
- Press “AZ = 3 stiffness”. The jump rule is unchanged, but the energy-orientation curve changes shape.
Computation
The interior field is the angular integral of Step 1 with 4096 midpoint samples (2048 for the sweeps), which resolves aspect ratios down to 0.01. Jumps use Step 2 exactly. The energy sweeps recompute Step 1 for 41 aspect ratios and 91 orientations.
The field maps use a 1024 × 1024 FFT grid by default (256 and 512 are selectable) on a periodic unit cell with an inclusion of semi-axis a = 0.16, a 4 × 4 supersampled indicator, and zero mean total strain. Periodic images and the finite area fraction shift the interior values by a few per cent relative to the infinite-matrix solution, and the sharp interface shows small Gibbs ripples. The separate sharp-interface profiles use the actual selected aspect ratio and a singularity-removing angular substitution, with 2048 midpoint samples per exterior point. The maps use b/a ≥ 0.08, the thinnest plate the grid resolves.
Elastic energy and directional dependence.
Change this module’s eigenstrain and elastic constants independently. The plots calculate the self elastic energy kernel \(B_{pp}(\mathbf n)\), using the same eigenstrain on both sides of the contraction. They do not plot a cross kernel or Young’s modulus. The polar angle denotes the Fourier wavevector direction: a modulation normal, not a plate tangent.
The eigenstrain is \(\boldsymbol\epsilon^{(p)}=\epsilon_p\operatorname{diag}(1,t)\). These inputs belong only to Module III. Equal components, \(t=1\), give in-plane dilatation. Unequal components can give a different angular symmetry, even with isotropic elastic constants. The self kernel has energy-density units for dimensional stiffnesses; this widget uses scaled model values.
Read this polar plot as an energy plot
The radius is proportional to \(B_{pp}(\mathbf n)\), with zero at the center. The concentric circles label kernel values. The drawing is scaled to fit the panel, so compare the labeled values when changing parameters. A minimum lies closer to the center. The previous offset scale placed zero away from the center and changed the apparent shape; the self kernel now uses a direct radial scale. Module IV retains its labeled signed scale because a cross kernel can be negative.
A surface-energy plot uses radius \(\gamma(\mathbf n)\), where \(\gamma\) is interfacial energy per area. This plot uses the elastic kernel, which depends on stiffness, eigenstrain and mechanical relaxation. Both are functions of direction, but they need not have the same angular shape. Neither raw polar curve is automatically the equilibrium particle shape: a surface-energy shape requires a Wulff construction, while an elastic particle shape requires the total elastic and interfacial energies.
Try the limiting case. Select equal dilatational components and isotropic elasticity. The kernel is independent of direction and the polar curve is a circle. Keeping isotropic stiffness but changing the strain ratio can restore angular variation. With cubic anisotropy and equal components, the curve has the corresponding fourfold symmetry. Ideal plate tangents are perpendicular to the plotted modulation normals.
Why did this happen?
Physics
Eigenstrain specifies what the phase wants to do. Stiffness specifies the cost of preventing it. Mechanical relaxation reduces that cost. Their combination gives B(n). A minimum identifies a low-penalty modulation normal; the corresponding ideal plate is perpendicular to it. Several equal minima are symmetry-related choices.
Mathematics
From my Appendix A, Eqs. (A.1–A.12), pp. 107–111. This is the elastic-energy derivation developed using Khachaturyan’s microelasticity framework. Spatial indices \(i,j,k,l\) are summed; phase indices \(p,q\) are summed only where shown. The laboratory uses two spatial dimensions.
1. Separate eigenstrain, mean strain and periodic relaxation.
For homogeneous stiffness and zero applied mean stress, \(E_{ij}=\langle\epsilon^0_{ij}\rangle\). The mean terms cancel in the six-integral expansion of Appendix A. Only the fluctuations \(\Delta\theta_p=\theta_p-\langle\theta_p\rangle\) remain.
2. Solve equilibrium for each nonzero Fourier mode.
The acoustic contraction leaves \(i,k\) free. It is not the contraction \(C_{ijkl}n_kn_l\), which leaves a different index pair free.
3. Substitute the relaxed displacement into the energy.
\(B_{pq}\) is the elastic energy kernel. Its first term is the unrelaxed eigenstrain contraction; the second is the energy reduction from mechanical relaxation. It becomes a total energy only after weighting by the field spectra and summing or integrating over modes. This is Eq. (A.12), with \(d=2\) here.
4. Differentiate the energy: Appendix B, Eqs. (B.1–B.12).
Here \(a\) is a fixed field label, not a summed spatial index. The two equal variations of the quadratic form cancel its factor of one half. We use \(B_{aq}=B_{qa}\), renaming dummy indices explicitly. For \(\theta_a=\eta_a^2\), the chain rule gives \(\delta F_e/\delta\eta_a=2\eta_a g_a\) in real space.
Boundary-condition distinction: Appendix A uses a relaxed mean strain. The finite-plate module fixes \(E_{ij}=0\) and therefore retains the additional mean energy \(\tfrac{V}{2}(E-\langle\epsilon^0\rangle):C:(E-\langle\epsilon^0\rangle)\). Its zero mode must not be silently dropped.
\(A_Z=2C_{44}/(C_{11}-C_{12})\). The 2-D law uses the in-plane cubic sub-block, with tensor shear \(\sigma_{xy}=2C_{44}e_{xy}\), no out-of-plane eigenstrain, and no plane-stress condensation.
Follow the full Appendix A and B derivations, including Julia and NumPy expressions →
Computation
For every sampled direction, form the 2 × 2 acoustic tensor Q, invert it and contract the stress vectors. Stability checks reject invalid constants. All minima come from evaluated curves. Flat curves are reported as directionally degenerate. Signed polar plots use an explicitly labeled offset radial scale; negative cross terms are never folded onto positive radii.
Sandeep’s thesis example: two precipitate phases in one matrix
This module introduces the three-phase system studied in Sandeep Sugathan’s thesis. The matrix is \(\alpha\); the two precipitate phases are \(\beta\) and \(\gamma\). Both precipitates have a dilatational misfit relative to \(\alpha\): their stress-free size differs from that of the matrix, without a preferred strain direction. In the X2 example, \(\beta\) expands by 1% and \(\gamma\) contracts by 1%.
This part of our work is described clearly in Sandeep Sugathan’s thesis. We have also written it up in our arXiv preprint, Effect of elastic anisotropy on phase separation in ternary alloys: A phase-field study, by Sandeep Sugathan and Saswata Bhattacharya. We are continuing this work by examining additional cases.
The elastic moduli are homogeneous and anisotropic. Homogeneous means that the same stiffness tensor \(C_{ijkl}\) is used in all three phases. Anisotropic means that the elastic response depends on direction. These are different properties: equal moduli across phases do not mean isotropic moduli.
The thesis reference is Chapter 5, Tables 5.1–5.2, Eq. (5.1) and Fig. 5.18. The parameters are \(G=2000\), \(\nu=1/3\), \(A_Z=3\), \(\epsilon_\beta=0.01\) and \(\epsilon_\gamma=-0.01\). Here \(G\) is the modulus parameter used in the thesis, \(\nu\) is Poisson’s ratio, \(A_Z\) is the Zener anisotropy ratio, and \(\epsilon_p\) is the stress-free linear misfit of precipitate phase \(p\) relative to \(\alpha\). Strains, \(\nu\) and \(A_Z\) are dimensionless. The stiffness and energy values below use the thesis’s model units; they are not labeled as pascals or joules without a dimensional conversion.
Set the elastic constants
Change each stiffness independently here. Alternatively, change the shear parameter \(\mu\) and Zener anisotropy ratio \(A_Z\), using a fixed Poisson parameter \(\nu=1/3\). The shear-parameter conversion uses a fixed \(\nu=1/3\). Direct stiffness inputs remain independent; an arbitrary triple need not correspond to that fixed Poisson parameter. These constants belong only to Module IV.
The stiffness constants and \(\mu\) use the same energy-density scale. This demonstration uses scaled values; no conversion to pascals is assumed. For an anisotropic solid, \(\mu\) and \(\nu\) are the parameters in the conversion below, rather than a shear modulus and Poisson ratio for every direction. When \(A_Z=1\), they reduce to the usual isotropic values.
Use \(\mu>0\), \(-1<\nu<1/2\), and \(A_Z>0\). Independent stiffness inputs must satisfy \(C_{44}>0\), \(C_{11}-C_{12}>0\), and \(C_{11}+2C_{12}>0\). Invalid inputs produce an explanation and retain the last valid plots.
At \(A_Z=3\), a dilatational self kernel favors Fourier directions \(\langle10\rangle\). Opposite misfits reverse the sign of the cross kernel, whose minimum is then along \(\langle11\rangle\). The pair-energy graph separately calculates the preferred separation direction by summing all Fourier modes.
Explore the three kernels
Select Sandeep thesis example to load the anisotropic reference case. Select Isotropic comparison to retain the misfits and set \(A_Z=1\). The stiffness controls above affect only Module IV. Misfits and strain ratios update the plots as you type. Keeping both ratios at one gives dilatational misfit; changing a ratio introduces non-dilatational strain and leaves the thesis example’s assumptions.
Read the signed values. Each polar panel displays one kernel: \(B_{\beta\beta}\), \(B_{\gamma\gamma}\), or \(B_{\beta\gamma}\). All three use the same scale. The bold circle marks zero; values inside it are negative and values outside it are positive. Moving inward means a lower signed value. Green dots mark the minimum. The angle is the Fourier normal direction \(\mathbf n\), measured from [10].
For equal and opposite misfits, the two self kernels are identical and the cross kernel is their negative. Separate panels let you see all three without hiding coincident curves. These are kernel plots; they do not show particle shapes.
Why did this happen?
When is the β–γ cross coefficient important?
The cross coefficient measures elastic coupling between the two phase fields. It is not the interaction energy of two particles. For the homogeneous elastic constants and dilatational misfits used here, all three kernels contain the same directional function:
Here \(\epsilon_\beta\) and \(\epsilon_\gamma\) are dimensionless linear misfits relative to the matrix. The unit vector \(\mathbf n\) gives the Fourier wavevector direction. The function \(K\) is the elastic kernel for unit dilatational misfit; it depends on the elastic constants and direction. The kernels have the dimensions of energy per volume when dimensional stiffnesses are used. This widget displays them in the model energy-density units defined above.
Thus the cross coefficient is comparable to both self coefficients when the misfit magnitudes are comparable. For Sandeep’s equal-and-opposite example, \(\epsilon_\beta=0.01\) and \(\epsilon_\gamma=-0.01\), so \(B_{\beta\gamma}=-B_{\beta\beta}=-B_{\gamma\gamma}\) at every direction. There is no anisotropy threshold at which the cross coefficient suddenly appears.
A large cross coefficient does not by itself establish a large particle interaction. The interaction energy also depends on the amounts, shapes, separation and arrangement of both phase fields. Contributions from different Fourier modes can reinforce or cancel. Elastic anisotropy gives the kernel angular variation; a nonzero isotropic cross coefficient alone does not establish an arrangement-dependent interaction. The Bitter–Crum derivation below explains the isotropic limit.
Now calculate a particle-pair energy.
The thesis preset selects a β–γ pair with opposite misfits. For this case at \(A_Z=3\), [11] has lower interaction energy. A β–β or γ–γ pair has same-sign misfits and gives a different preference. The positive self kernels are higher along [11]; the signed β–γ cross kernel is lower there. The three polar panels use signed radial values; a smaller radius means a lower kernel value.
Why did this happen?
Physics
Changing one misfit sign reverses the cross coupling because it is bilinear in the two eigenstrains. Both self terms remain nonnegative. However, \(B_{\beta\gamma}(\mathbf n)\) is a reciprocal-space coefficient. To obtain real-space interaction energy at separation R, we must also include particle shape spectra and Fourier phases. Single-particle energy and pair-interaction energy are different quantities.
Mathematics
The \(\mathbf k=0\) coefficient in this clamped periodic calculation is \(\epsilon^\beta:C:\epsilon^\gamma\). The square contains identical circular Gaussian profiles h, with width 0.045 of the box. Each Fourier mode contributes; the kernel is not simply evaluated at the separation angle.
Computation
Use a 64 × 64 periodic square and exact Fourier translations of the sampled Gaussian. Evaluate a 1° separation-angle sweep. Periodic copies are included. The profiles have diffuse tails. The on-screen circles indicate their positions, not sharp boundaries.
Thesis reference: Sandeep Sugathan, Chapter 5, Tables 5.1–5.2 (p. 116), Eq. (5.1) and Fig. 5.18 (p. 133), discussion p. 132. The X2 example has G = 2000, ν = 1/3, AZ = 3 and misfits ±0.01. Eqs. (3.27–3.29) give \(C_{11}=7000\), \(C_{12}=5000\), \(C_{44}=3000\). Press “Sandeep thesis example” to see the comparison table: the self terms are 0.342857 along ⟨10⟩ and 0.8 along ⟨11⟩, and the cross term is their negative. The separate two-Gaussian experiment is a teaching calculation, not a reproduction of the thesis’s full phase-field microstructure.
When does the interaction control the arrangement?
For the two prescribed circular profiles in this widget, rotating their separation leaves each self energy unchanged. The angle-dependent interaction therefore selects the preferred separation direction. In the thesis example, opposite misfits with \(A_Z=3\) give the lowest pair interaction along [11]. The self kernels still favor Fourier normals along [10].
Use the pair calculation to compare actual arrangements, rather than reading a pair direction directly from the kernel polar plot. Define
Here \(E_{\mathrm{int}}\) is the interaction energy reported by this widget in its normalized model units. The directions [10] and [11] describe the line joining the particle centers; they are different from the Fourier normal \(\mathbf n\) in the kernel plots. A negative difference favors [11], a positive difference favors [10], and zero gives no preference between these two arrangements. This comparison keeps particle size, separation and elastic constants fixed. It does not exclude a lower energy at another angle.
In an evolving microstructure, the precipitates can also change shape and composition. The change in interaction energy must then be compared with changes in self elastic energy, interfacial energy and chemical free energy. Larger misfits increase elastic effects relative to the nonelastic contributions when other parameters are fixed. Both precipitate phases must have appreciable misfits and spatially coupled fields for their cross contribution to be important. Increasing anisotropy can strengthen a directional preference, but it does not give a universal threshold for interaction dominance.
For one Fourier mode, define \(S_\beta=B_{\beta\beta}|\hat h_\beta|^2\), \(S_\gamma=B_{\gamma\gamma}|\hat h_\gamma|^2\), and \(X=2B_{\beta\gamma}\operatorname{Re}(\hat h_\beta\hat h_\gamma^*)\). The mode energy is \((S_\beta+S_\gamma+X)/2\). The cross contribution is most significant when both self contributions are comparable and the Fourier phases give strong correlation. Mechanical stability gives \(|X|\leq S_\beta+S_\gamma\): the interaction can control the preferred arrangement without exceeding the sum of the self contributions in magnitude. The full phase-field calculation is needed to decide the final morphology.
Derivation and the Bitter–Crum limit
1. Write the transformation strain of each phase
Let \(h_\beta(\mathbf x)\) and \(h_\gamma(\mathbf x)\) be dimensionless fields describing where the precipitates are located. A field is one inside its phase and zero outside for a sharp interface; a smooth field can represent a diffuse interface. The position is \(\mathbf x\). Taking the matrix as the zero-misfit reference gives
Here \(\epsilon^0_{ij}\) is the prescribed transformation strain, and \(\delta_{ij}\) is one when \(i=j\) and zero otherwise. The indices identify Cartesian components. A three-dimensional dilatation has three equal diagonal components. This widget uses the two-dimensional in-plane elastic law described in Module III: \(\boldsymbol\epsilon^p=\epsilon_p\operatorname{diag}(1,1)\). It omits the out-of-plane transformation strain; it is not a full three-dimensional plane-strain reduction.
2. Impose mechanical equilibrium
The displacement components are \(u_i\). The total strain is \(\epsilon_{ij}=(\partial_i u_j+\partial_j u_i)/2\), where \(\partial_i\) denotes differentiation with respect to position component \(x_i\). Elastic strain is total strain minus transformation strain. Stress \(\sigma_{ij}\) and elastic energy \(E_{\mathrm{el}}\) are therefore
Repeated indices are summed. \(V\) is the body volume; the two-dimensional calculation uses area and energy per unit out-of-plane thickness. Stiffness and stress have units of pressure in a dimensional model, and the integral then has units of energy. Mechanical equilibrium assumes no body force. The displacement adjusts to reduce the energy while satisfying the boundary conditions.
3. Solve each nonzero Fourier mode
A Fourier mode has wavevector \(\mathbf k\), with units of inverse length in a dimensional model. Its unit normal is \(\mathbf n=\mathbf k/|\mathbf k|\). Define the stress associated with the prescribed strain, the acoustic tensor, and the inverse acoustic tensor by
\(Q\) describes the directional stiffness of the mode. \(\Omega\) has units of inverse stiffness. With the Fourier convention \(\widehat{\partial_j u_i}=\mathrm i k_j\hat u_i\), mechanical equilibrium gives
Here \(\mathrm i^2=-1\), a hat denotes a Fourier coefficient, and \(\hat\sigma^0_{ij}=\sum_p\sigma^p_{ij}\hat h_p\). Substituting this equilibrium displacement into the energy removes the displacement as an independent unknown.
4. Obtain the two self kernels and one cross kernel
The first term is the energy coefficient before displacement relaxation. The second term is the reduction produced by mechanical relaxation. The vector \(a^p_i\) is the traction associated with the prescribed strain on a plane normal to \(\mathbf n\). Each \(B_{pq}\) has units of energy density in a dimensional model.
For clarity, define \(\hat h_p(\mathbf k)=V^{-1}\int_V h_p(\mathbf x)e^{-\mathrm i\mathbf k\cdot\mathbf x}\,dV\) in a periodic body. These coefficients are dimensionless. The nonzero-mode energy density is
\(\operatorname{Re}\) means real part, and an asterisk means complex conjugate. The two self kernels describe each precipitate phase’s contribution. The cross kernel describes the coupling between their fields. Symmetry gives \(B_{\beta\gamma}=B_{\gamma\beta}\), so there are three independent kernels, although the system contains three phases. The matrix contributes through the common stiffness and the reference strain.
For the two-dimensional square-symmetric stiffness, write \(n_x=\cos\theta\) and \(n_y=\sin\theta\), where \(\theta\) is the wavevector-normal angle. The matrix actually inverted by the widget is
The thesis parameters give \(C_{11}=7000\), \(C_{12}=5000\), \(C_{44}=3000\) in model units. The equal and opposite misfits imply \(B_{\beta\beta}=B_{\gamma\gamma}\) and \(B_{\beta\gamma}=-B_{\beta\beta}\). For this in-plane calculation, the self coefficient is \(0.342857\) along [10] and \(0.8\) along [11]; the cross coefficient is their negative. These numbers check the kernel calculation, not the reproduction of a complete thesis microstructure.
5. Derive the Bitter–Crum theorem from the isotropic kernel
The isotropic example demonstrates the Bitter–Crum theorem. We can show why it holds by calculating \(B_{pq}\). Assume that the elastic moduli are the same in all phases, the elasticity is isotropic, and each phase has a purely dilatational transformation strain. Keep the phase fractions and misfits fixed. Use an infinite body or the compatible periodic boundary conditions described below.
For the two-dimensional in-plane law used here, let \(\lambda=C_{12}\) and \(\mu=C_{44}\) be the Lamé elastic constants. Both have units of pressure in a dimensional model. Isotropy requires \(C_{11}=\lambda+2\mu\), or \(A_Z=1\). Since \(\epsilon^p_{ij}=\epsilon_p\delta_{ij}\), the prescribed stress is \(\sigma^p_{ij}=2(\lambda+\mu)\epsilon_p\delta_{ij}\). Thus
Here \(I\) is the identity matrix, and \(\mathbf n\otimes\mathbf n\) has components \(n_i n_j\). Substitute these expressions into the kernel. The two terms give
The coefficient \(K_B\) has units of pressure. There is no direction \(\mathbf n\) left in this expression. Equal and opposite misfits give \(B_{\beta\gamma}=-B_{\beta\beta}\), but this negative coefficient alone does not imply attraction between precipitates. We must substitute it into the complete energy expression.
Define the scalar misfit field \(e(\mathbf x)=\epsilon_\beta h_\beta(\mathbf x)+\epsilon_\gamma h_\gamma(\mathbf x)\). It is dimensionless. Parseval’s identity states that the sum of the squared Fourier coefficients equals the spatial average of the squared field. Removing the zero mode therefore gives
Angle brackets denote the spatial average. For disjoint phase regions, \(h_p\) is one inside phase \(p\) and zero outside. Consequently \(h_p^2=h_p\) and \(h_\beta h_\gamma=0\). If \(f_p=\langle h_p\rangle\) is the dimensionless phase fraction, then
This is the Bitter–Crum result in this model. The energy depends on the phase fractions and misfits. It does not depend on precipitate shape, position, separation or arrangement. Moving disjoint precipitates while keeping their amounts fixed does not change the energy: there is no elastic interaction driving their relative motion. Opposite misfits do not change this conclusion.
For imposed mean strain \(\bar\epsilon\), add \(E_0/V=\tfrac12(\bar\epsilon-\langle\epsilon^0\rangle):C:(\bar\epsilon-\langle\epsilon^0\rangle)\). The colons denote tensor contraction. This term also depends only on the mean transformation strain when the imposed strain is fixed. At relaxed zero mean stress, it vanishes. Comparing arrangements at fixed phase amounts is essential; a cross term obtained by subtracting two separate clamped calculations is not, by itself, a separation-dependent interaction.
Select Isotropic comparison to check the constant kernel numerically. Its circle and the horizontal signed curves demonstrate the directional cancellation derived above; the energy reduction proves the shape-and-arrangement independence. With anisotropic moduli, the kernel generally retains its directional dependence, and this Bitter–Crum proof no longer applies. Eshelby’s uniform-interior-strain theorem for an ellipsoidal inclusion is a separate result.
Soft or hard: rotate a fixed plate.
For this prescribed thin plate, what orientation minimises the elastic energy? The crystal axes and eigenstrain stay fixed while the plate shape rotates. Choose β or γ; each retains its own stiffness and eigenstrain settings. This module contains one plate at a time.
The selected phase uses the misfit inputs in this module. Its matrix and precipitate constants are also independent of every other module.
Choose parameters and calculate a sweep.
Why did this happen?
Physics
Soft and hard refer to the elastic response relative to the matrix. The same eigenstrain is accommodated differently when stiffness varies in space. With isotropic phases and dilatational misfit, an infinite medium has no preferred material direction. Small angular differences in this periodic calculation require particular care: the square cell and finite grid also introduce directions. We fix the plate shape to isolate its orientation energy, rather than claim to predict its unconstrained equilibrium morphology.
Physical insight: A minimum is the smallest elastic penalty among the sampled orientations of this prescribed plate, under the stated boundary condition. It is not automatically a habit plane or a stress-free state.
Mathematics
Inhomogeneous equilibrium: Bhattacharyya, Heo, Chang and Chen, Communications in Computational Physics 11 (2012), 726–738, Eqs. (2.6–2.16). Tushar Jogi’s thesis, §§3.1.4–3.1.6 and Appendix B, supplies the related phase-field formulation and driving-force derivation. Indices are retained here so each contraction can be checked.
The spectral iteration, with every index visible.
The wave-number factor follows directly by transforming two spatial derivatives. If the Green tensor is defined using \(k_jk_l\), its inverse already contains \(|\mathbf k|^{-2}\). If it is defined using the unit direction \(n_jn_l\), that factor must appear explicitly. The zero mode is handled separately.
Mean strain is a boundary condition.
For prescribed mean stress, invert \(\langle C\rangle\) on symmetric tensors. In general \(\langle C\rangle^{-1}\ne\langle C^{-1}\rangle\). This is the meaning of the compliance-like quantity defined in CICP Eq. (2.16). The current plate experiment instead prescribes \(E=0\); it does not claim to solve a traction-free macroscopic cell.
Elastic driving force, with the interpolation distinguished from the kernel.
This is the first-order result of Tushar’s Appendix B, pp. 181–182, written in indicial form. The displacement variation vanishes at mechanical equilibrium under the chosen boundary conditions. For this laboratory, \(C=C^m+h(C^p-C^m)\) and \(\epsilon^0=h\epsilon^p\), so \(g_h=\tfrac12 e:(C^p-C^m):e-\sigma:\epsilon^p\). The scalar interpolation called \(B(\phi)\) in Tushar’s thesis is not the elastic energy kernel \(B_{pq}(\mathbf n)\).
CICP paper, DOI and publication details · Equation-by-equation verification notes
Computation
The production solver uses preconditioned conjugate gradients rather than the paper’s fixed-point iteration. Both solve the same heterogeneous equilibrium equation. An independent CICP-style spectral iteration is included in the validation suite for comparison. A spectral symmetric-gradient operator and its adjoint form A = D* C(x) D. Preconditioned conjugate gradients solve A u = D* C(x) ε⁰. A homogeneous acoustic inverse is only a preconditioner; the residual uses the full heterogeneous stiffness. Failed solves are rejected, not plotted as minima.
The unit square contains a smooth ellipse with semiaxes 0.22 and 0.045. Its diffuse area is held at π × 0.22 × 0.045 for every orientation. The tanh transition width is 0.18 in elliptic-radius units. Fourier derivatives on even grids set Nyquist derivatives to zero, consistently in both the solver and reference. These modes retain their unrelaxed energy. Refine the grid to check their influence and resolve the interface.
Circles on the plots mark sampled local minima; larger circles mark sampled global minima. Recalculate at 64² or 128² and a smaller angular step to assess convergence. The solver reports its residual and iterations. The “Homogeneous limit” copies matrix constants to the selected plate and compares against an independent Fourier-energy evaluation.
Equations, numerical checks and source code
I compare the calculated strains and elastic energies with known solutions and with results from separate implementations. I also repeat selected calculations at finer grid sizes and smaller angle steps to check how much the results change.
The saved numerical results and explanations of the checks are available below.
- Numerical checks and results: what was compared, the differences found and what the checks establish.
- Equations and assumptions: the quantities, units, equations and boundary conditions used in each calculation.
- References and notation: the books, papers and theses used, with explanations of the notation.
- How to use and run the laboratory: the controls, source files and commands for running the checks.
- Download the source code and check programs.
This laboratory is part of my teaching and research in mesoscale microstructure modeling. It connects eigenstrain, elastic equilibrium and elastic energy with microstructure, complementing our teaching notebooks and research codes for homogeneous and inhomogeneous elasticity, phase-field evolution and ferroelectric modeling.
Explore the mesoscale modeling resources and elastic stress effects examples. Source code, numerical tests and licensing details are provided in the documentation above.
Read, run and change the elasticity solvers.
Explore the corrected course notebook Elastic.ipynb, a direct iterative perturbative solver and an equivalent-eigenstrain solver. Each notebook derives the equations, exposes the iteration and lets you change the physical parameters.
Open the explained solver notebooks →
My elasticity research codes
My historical hexagonal-to-orthorhombic codes and Julia thesis model describe elastic effects on microstructural evolution. The historical repository includes the GSL version. The Julia model includes conserved noise in composition, separate nonconserved noise in all three variants, and my particle-moment utility for measuring changes in particle shape. Each repository describes its equations and tests.
I discuss phase separation inside a particle with the Cahn–Hilliard example, including the role of the field that distinguishes the particle from its surroundings.
From my thesis: a growing precipitate tilts.
Sources: Saswata Bhattacharyya, Evolution of Multivariant Microstructures with Anisotropic Misfit: A Phase Field Study, Appendix A, Eqs. (A.1–A.12), pp. 107–111; Appendix B, Eqs. (B.1–B.12), pp. 112–114; and the habit-line discussion in Chapter 4. These derivations follow A. G. Khachaturyan, Theory of Structural Transformations in Solids (Wiley, 1983). Other reference examples: Sandeep Sugathan, A Phase-Field Study of Elastic Stress Effects on Phase Separation in Ternary Alloy Systems, IIT Hyderabad, June 2019, Chapters 3 and 5 as cited in Module IV. Read the elastic-kernel and driving-force derivations.