Representation of Equilibrium Shapes, γ-plots and Stiffness

J. Friedli
EPFL-STI-IMX-LSMX
Station 12
CH-1015 Lausanne
jonathan.friedli@epfl.ch

Introduction

In solidification of metallic alloys, dendrites have a tendency to develop along directions corresponding to convexities (maximum curvature regions) of the equilibrium shape of crystal. Surface tension (solid-liquid interfacial energy) anisotropy seems to play a key role in the selection of those growth directions.  In fact a variation of this anisotropy leads to a shift in growth directions. Such a variation can for example be caused by a change in composition, as in Al-Zn alloys where dendrites progressively transit from a <100> to a <110> growth directions with an increasing zinc content [1].

To efficiently predict the solidification microstructure and thus the resulting mechanical properties, a computational model like Phase Field requires quantification of the surface tension in 3D. The determination of its anisotropy, conditioned by the crystal structure, is done by measuring the equilibrium shapes of liquid droplets enclosed in a solid. Wulffs [2] theory states that if a crystalline body is in its equilibrium shape, there exists a point whose perpendicular distances from the surfaces of the body are in proportion to their free energies. With this information it is then possible to obtain the so called γ-plot, i.e. the value of surface tension as function of crystallographic orientation.

Mathematica will be used to plot Equilibrium Shape (r-plot), γ-plots and stiffness as function of the crystallographic orientation, for different sets of ((Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]))

Method

The influence of the interfacial energy on the equilibrium of the solid-liquid interface can be quantified with the generalized Gibbs-Thomson relation which link undercooling, surface tension and curvature. According to Herrings [3] relation the curvature undercooling can be written:

(Subscript[ΔT, r])((OverVector[n]))=(1/(Subscript[ΔS, f]))((Subscript[κ, 1])((Subscript[γ, sl])((OverVector[n]))+(((∂^2) (Subscript[γ, sl])((OverVector[n])))/(∂(Subscript[θ, 1]^2))))+(Subscript[κ, 2])((Subscript[γ, sl])((OverVector[n]))+(((∂^2) (Subscript[γ, sl])((OverVector[n])))/(∂(Subscript[θ, 2]^2)))))

where, (OverVector[n]) is the normal vector to the interface, (Subscript[ΔS, f]) the entropy of melting and κ the principal curvatures of the interface. The terms in parentheses are the two principal interfacial stiffness components which can be normalized and generalized to one term for a locally spherical interface:

(Subscript[Φ, r])((OverVector[n]))=(1/(Subscript[γ, sl]^0))((Subscript[γ, sl])((OverVector[n]))+(1/2)((((∂^2) (Subscript[γ, sl])((OverVector[n])))/(∂(Subscript[θ, 1]^2)))+(((∂^2) (Subscript[γ, sl])((OverVector[n])))/(∂(Subscript[θ, 2]^2)))))

The driving force for solidification will be highest where the curvature is lowest, i.e. at the minima of the interface stiffness. Thus one needs to know (Subscript[γ, sl]). Mathematically the orientation dependence of (Subscript[γ, sl]) in three dimensions is described as an expansion of spherical harmonics satisfying the cubic symmetry (for the example of FCC aluminum on which this project focuses). As suggested by Fehlner and Vosko [4] and with the more recent convention by Asta and al. [5], the anisotropy of the surface tension can be expressed as:

(Subscript[γ, sl])((OverVector[n]))=(Subscript[γ, sl]^0)[1+(Subscript[ϵ, 1])(Q-(3/5))+(Subscript[ϵ, 2])(3Q+66S-(17/7))+(Subscript[ϵ, 3])(65(Q^2)-94Q-208S-33)+...]

where :

Q=(Subscript[n, x]^4)+(Subscript[n, y]^4)+(Subscript[n, z]^4)

S=(Subscript[n, x]^2)(Subscript[n, y]^2)(Subscript[n, z]^2)

(Subscript[n, x]), (Subscript[n, y]) and (Subscript[n, z]) are the Cartesian components of the normal unit vector to the solid-liquid interface expressed in the crystallographic coordinate system. For this work, the expansion to the third term of the expression of (Subscript[γ, sl]) is sufficient.

Finally we use the ξ-plot introduced by Cahn and Hofmann [6], which is geometrically similar to the equilibrium shape, to represent the solid-liquid interface. Therefore, we introduce the ξ-vector defined as:

ξ=grad(r (Subscript[γ, sl])((OverVector[n])))

we thus get that :

ξ((OverVector[n]))=(Subscript[γ, sl]) (OverVector[n]) + (Subscript[(((∂ (Subscript[γ, sl]))/(∂θ ))), max])(Subscript[(OverVector[t]), max])=(Subscript[γ, sl]) (OverVector[n])+((∂ (Subscript[γ, sl]))/(∂θ ))(Subscript[(OverVector[t]), θ])+(1/(sin(θ)))((∂ (Subscript[γ, sl]))/(∂ϕ ))(Subscript[(OverVector[t]), ϕ])

with (Subscript[(OverVector[t]), θ]) and (Subscript[(OverVector[t]), ϕ]) the tangential unitary vector along θ and ϕ.

(Subscript[γ, sl])((OverVector[n])), (Subscript[Φ, r])((OverVector[n])) and ξ((OverVector[n])) can now be used to represent the interfacial energy, the stiffness and the equilibrium shape for a given set of parameters (Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]).

Results

The results are separated in two kinds of representations. First the 3D plots, where the normalized value plotted in a color scale is superimposed on the rendered surface to aid visualization and second the 2D cuts in some important crystallographic sections for the cubic system.

3D representation

Given the symmetry of the cubic crystal the whole information, meaning anisotropy, stiffness and equilibrium shape is contained in only 1/48 of the sphere, i.e. a portion between the <100>,<110> and <111> directions. For representation purpose however, we will use 1/8 of the sphere, allowing both a lower computation time and an instinctive representation.

A important element for the calculation is that in spherical coordinates, the angle ϕ doesn't run along a great circle. This means that for calculation of Equation (2) the derivative with respect to the second principal curvature (i.e. not θ) must be corrected for the angle ϕ to be used. In numerical derivation, dϕ', the angle increment along the second great circle is given by ArcTan[Tan[dϕ]Sin[θ]]. For the analytical derivative one does have:

(((∂^2) (Subscript[γ, sl])((OverVector[n])))/(∂(Subscript[θ, 2]^2)))=(1/((Sin[θ])^2))(((∂^2) (Subscript[γ, sl])((OverVector[n])))/(∂(ϕ^2)))+(1/(Tan[θ]))((∂ (Subscript[γ, sl])((OverVector[n])))/(∂θ ))

thus Equation (2) becomes :

(Subscript[Φ, r])((OverVector[n]))=(1/(Subscript[γ, sl]^0))((Subscript[γ, sl])((OverVector[n]))+(1/2)((((∂^2) (Subscript[γ, sl])((OverVector[n])))/(∂(θ^2)))+(1/((Sin[θ])^2))(((∂^2) (Subscript[γ, sl])((OverVector[n])))/(∂(ϕ^2)))+(1/(Tan[θ]))((∂ (Subscript[γ, sl])((OverVector[n])))/(∂θ ))))

To make this notebook a powerful representation tool, one introduces the new functions from Mathematica 6, Manipulate, SphericalPlot3D and ColorFunction. This allows to evaluate the cell leading to the representation only once and then change the desired parameters, here (Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]), with rulers. The default values of these parameters have been chosen as calculated by Niederberger and al.  [7] and match their results as can be seen by comparing Figure 1. and Figure 2.. They also match the earlier results obtained by Napolitano and Liu [8] who used only two terms in the spherical harmonics expansion ((Subscript[ϵ, 1])=0.0181 and (Subscript[ϵ, 2])= -0.0112 in Figure 3.) .

Off[General::Stop,General::Spell,Spell1,Power::infy,∞::indet,FindMinimum::"lstol",FindMaximum::"lstol"]

nx:=Cos[ϕ] Sin[θ];
ny:=Sin[ϕ] Sin[θ];
nz:=Cos[θ];
Q:=(nx^4)+(ny^4)+(nz^4)
S:=(nx^2)(ny^2) (nz^2)

γ[ϵ1_, ϵ2_,ϵ3_]:=1.+(ϵ1(Q-(3/5))+ϵ2(3 Q + 66 S -(17/7))+ϵ3 (65 (Q^2)-94 Q - 208 S + 33));

g[ϵ1_, ϵ2_,ϵ3_]:=1+10*(ϵ1(Q-(3/5))+ϵ2(3 Q + 66 S -(17/7))+ϵ3 (65 (Q^2)-94 Q - 208 S + 33));

d1γθ=D[γ[ϵ1, ϵ2,ϵ3], {θ,1}];

d2γθ=D[γ[ϵ1, ϵ2,ϵ3], {θ,2}];

d1γϕ=D[γ[ϵ1, ϵ2,ϵ3], {ϕ,1}];

d2γϕ=D[γ[ϵ1, ϵ2,ϵ3], {ϕ,2}];

Φ[ϵ1_, ϵ2_,ϵ3_]=γ[ϵ1, ϵ2,ϵ3]+(1/2)(d2γθ+d2γϕ/((Sin[θ])^2)+d1γθ/Tan[θ]);

esx=γ[ϵ1, ϵ2,ϵ3] nx+d1γθ Cos[ϕ] Cos [θ]-(1/(Sin[θ]))d1γϕ Sin[ϕ];
esy=γ[ϵ1, ϵ2,ϵ3] ny+d1γθ Sin[ϕ] Cos [θ]+(1/(Sin[θ]))d1γϕ Cos[ϕ];
esz=γ[ϵ1, ϵ2,ϵ3] nz-d1γθ Sin [θ];

Es[ϵ1_, ϵ2_,ϵ3_]=(Sqrt[(esx^2)+(esy^2)+(esz^2)]);

SetOptions[SphericalPlot3D,ViewPoint-> {5,5,5},
ImageSize->{200},BoundaryStyle->Directive[Black],Mesh->False,Boxed->False, Ticks-> False,PlotRange-> {0,1.25},AxesEdge->{{0,1.5},{0,1.5},{0,1.5},AxesOrigin->{{0},{0},{0}}}];

SetOptions[DensityPlot,PlotRangePadding->0.0005,ImageSize->{190},AspectRatio-> 1/5,FrameTicks->{{None,None},{Automatic,None}},ColorFunctionScaling->True,LabelStyle->Directive[Bold,FontSize->9]];

Manipulate[Grid[{{
SphericalPlot3D[Evaluate@γ[e1, e2, e3], {θ,0,π/2},{ϕ,0,π/2},PlotPoints->ControlActive[5,20],ColorFunction->Function[{x, y,z,θ,ϕ,γ},ColorData["Rainbow"][γ]]],SphericalPlot3D[Evaluate@Φ[e1, e2, e3], {θ,0,π/2},{ϕ,0,π/2},PlotPoints->ControlActive[5,20],ColorFunction->Function[{x, y,z,θ,ϕ,Φ},ColorData["Rainbow"][Φ]]],SphericalPlot3D[Evaluate@Es[e1, e2, e3], {θ,0,π/2},{ϕ,0,π/2},PlotPoints->ControlActive[5,20],ColorFunction->Function[{x, y,z,θ,ϕ,Es},ColorData["Rainbow"][Es]]]},{DensityPlot[x,{x,FindMinimum[γ[e1,e2,e3],{θ,π/7},{ϕ,π/8}][[1]],FindMaximum[γ[e1,e2,e3],{θ,π/32},{ϕ,π}][[1]]},{y,0,0.002},ColorFunction->ColorData["Rainbow"]],DensityPlot[x,{x,FindMinimum[Φ[e1,e2,e3],{θ,π/7},{ϕ,π/8}][[1]],FindMaximum[Φ[e1,e2,e3],{θ,0.955},{ϕ,π/4}][[1]]},{y,0,0.002},ColorFunction->ColorData["Rainbow"]],DensityPlot[x,{x,FindMinimum[Es[e1,e2,e3],{θ,π/7},{ϕ,π/8}][[1]],FindMaximum[Es[e1,e2,e3],{θ,π/32},{ϕ,π/4}][[1]]},{y,0,0.002},ColorFunction->ColorData["Rainbow"]]}}],
{{e1,0.025},0,0.5},{{e2,-0.0047},-0.02,0.02},{{e3,-0.000194},-0.001,0.001},FrameLabel->Style["Interfacial energy γ-plot                  Stiffness Φ-plot                   Equilibrium Shape ξ-plot", Medium]]

Figure 1. Plots of (Subscript[γ, sl])((OverVector[n])), (Subscript[Φ, r])((OverVector[n])) and ξ((OverVector[n])) versus crystallographic orientation with adjustable (Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]) parameters. Default values are (Subscript[ϵ, 1])=0.025, (Subscript[ϵ, 2])=-0.0047 and (Subscript[ϵ, 3])=-0.000194. The axes represent the <100> directions.

Figure 2. Plots of (Subscript[γ, sl])((OverVector[n]))/(Subscript[γ, sl]^0) and (Subscript[Φ, r])((OverVector[n])) versus crystallographic orientation for the values of anisotropy (Subscript[ϵ, 1])=0.025, (Subscript[ϵ, 2])=-0.0047 and (Subscript[ϵ, 3])=-0.000194 from Niederberger and al. [7].

Figure 3. Plots of (Subscript[γ, sl])((OverVector[n]))/(Subscript[γ, sl]^0) and (Subscript[Φ, r])((OverVector[n])) versus crystallographic orientation for the values of anisotropy (Subscript[ϵ, 1])=0.0181, (Subscript[ϵ, 2])=-0.0112 and (Subscript[ϵ, 3])=0 from Napolitano and Liu [8].

Niederberger and al., Napolitano and Liu have both observed what Figure 2. and 3. suggest, that dendrites grow at the minima of interface stiffness, i.e. in 8 specific directions (near <320>) in the (001) plane for Al-Zn-Si alloys [7] and in <110> in an Al-Sn alloy [8].

2D representation

Given the four-fold in plane symmetry of the crystal, one quadrant is sufficient to display all the information. Since anisotropy in metallic systems like aluminum based alloys is so small, the γ-plot can hardly be differentiated from a sphere, i.e. a circle in 2D. Thus for better visualization, the anisotropy parameters (Subscript[ϵ, 1]),(Subscript[ϵ, 2]),(Subscript[ϵ, 3]) have been multiplied by a factor 10 in Figure 4.. The stiffness plot being enough anisotropic, its effect on the Φ-plot didn't have to be exaggerated.

The Mathematica function Manipulate is once more used to make this notebook as interactive as possible.

The default parameters of Figure 4. are once more those calculated by Niederberger and al. [7] and the plots are in agreement with both their work and the paper by Napolitano and Liu [8].

SetOptions[PolarPlot,ImageSize->{190,200},PlotRange->{0,1.3},PlotStyle->{Directive[Thick,Red],Directive[Thick,Blue],Black,Black,Black},Ticks->{{0.0,0.5,1.0,1.25},None},TicksStyle->Directive[10]];
SetOptions[ParametricPlot,ImageSize->{190,200},PlotRange->{0,1.3},PlotStyle->{Directive[Thick,Red],Directive[Thick,Blue],Black,Black,Black},Ticks->{{0.0,0.5,1.0,1.25},None},TicksStyle->Directive[10]];

Θ[ϕ_]:=ArcTan[((-1)/((Sqrt[2])Cos[ϕ-π/4]))]/;7π/4>ϕ∧ϕ>3π/4
Θ[ϕ_]:=π+ArcTan[((-1)/((Sqrt[2])Cos[ϕ-π/4]))]/;ϕ<3π/4 ∨ϕ>7π/4
Θ[ϕ_]:=π/2/;ϕ==3π/4∨ϕ==7π/4

P111={{((-1)/(Sqrt[2])),((-1)/(Sqrt[6])),(1/(Sqrt[3]))},{(1/(Sqrt[2])),((-1)/(Sqrt[6])),(1/(Sqrt[3]))},{0,(Sqrt[(2/3)]),(1/(Sqrt[3]))}};

Mrot111=Inverse[P111];

gnew[ϵ1_,ϵ2_,ϵ3_]:=Mrot111.((g[ϵ1,ϵ2,ϵ3]/.{θ->Θ[ϕ]}){Sin[Θ[ϕ]]Cos[ϕ],Sin[Θ[ϕ]]Sin[ϕ],Cos[Θ[ϕ]]})//Simplify

fnew[ϵ1_,ϵ2_,ϵ3_]:=Mrot111.((Φ[ϵ1,ϵ2,ϵ3]/.{θ->Θ[ϕ]}){Sin[Θ[ϕ]]Cos[ϕ],Sin[Θ[ϕ]]Sin[ϕ],Cos[Θ[ϕ]]})//Simplify

Manipulate[Grid[{{
PolarPlot[{g[e1,e2,e3]/.θ->π/2,Φ[e1,e2,e3]/.θ->π/2,0.5,1,1.25},{ϕ,0,4π/2},AxesLabel->{Style["[100]",Bold],Style["[010]",Bold]}],PolarPlot[{g[e1,e2,e3]/.{ϕ->π/4,θ->π/2-q},Φ[e1,e2,e3]/.{ϕ->π/4,θ->π/2-q},0.5,1,1.25},{q,0,π/2},AxesLabel->{Style["[(OverBar[1])10]",Bold],Style["[001]",Bold]}],
ParametricPlot[{Drop[gnew[e1,e2,e3],-1],Drop[fnew[e1,e2,e3],-1],{0.5 Cos[ϕ],0.5 Sin[ϕ]},{ 1Cos[ϕ],1 Sin[ϕ]},{1.25 Cos[ϕ],1.25 Sin[ϕ]}},{ϕ,0,2π},AxesLabel->{Style["[(OverBar[1])10]",Bold],Style["[(OverBar[1])(OverBar[1])2]",Bold]}]},{Style["(100) Section",Medium,Bold],Style["(110) Section",Medium,Bold],Style["(111) Section",Medium,Bold]}}],
{{e1,0.025},0,0.5},{{e2,-0.0047},-0.02,0.02},{{e3,-0.000194},-0.001,0.001}]

Figure 4. Two dimensional polar plots of the interfacial energy (red, exaggerated ten times) and stiffness (blue) in three important cyrstallographic sections of the cubic system. Default values are (Subscript[ϵ, 1])=0.025, (Subscript[ϵ, 2])=-0.0047 and (Subscript[ϵ, 3])=-0.000194.

Conclusion and future work

The tool developed here will help to ascertain the assumption that dendrites grow at the minima of interface stiffness by correlating results of Gonzales and Rappaz [1] showing a growth direction change with zinc concentration in Al-Zn alloys and the measurements of the solid-liquid interfacial energy for various composition that will be done in this Ph. D. thesis.

In future development this Mathematica notebook will be extended to calculate the "real" stiffness, i.e. without the spherical assumption for (Subscript[γ, sl]), calculating the effective local principal curvatures. Further analysis will also emphasis on which parameter range leads to which "plane of growth" for dendrites. Finally a section for the visualization of the results in an inverse pole figure mode will probably be added.

References

[1] F.Gonzales and M.Rappaz, Metall. Mater. Trans. A 37 (2006) 2797

[2] G. Wulff, Z. Krist. 34 (1901) 449

[3] C. Herring, in The Physics of Powder Metallurgy, edited by W.E. Kingston (McGraw-Hill, NewYork, 1951)

[4] W.R. Fehlner and S.H. Vosko, Can. J. Phys. 54 (1976) 2159

[5] M. Asta, J.J. Hoyt and A. Karma, Phys. Rev. B 66, 100101 (2002)

[6] J.W. Cahn and D.W. Hoffman, Acta. Met. 22 (1974) 1205-1214

[7] C. Niederberger, J. Michler and A. Jacot, Phys. Rev. E 74, 021604 (2006)

[8] R.E. Napolitano and S. Liu, Phys. Rev. B 70, 214103 (2004)