Homogenization Model
Example - Fe-Cr-Ni system
The homogenization model can simulate multiphase diffusion without having to resort to more complex methods such as phase field modeling. The model relies on the assumption that every volume element is in local equilibrium. Then fluxes are determined by the mobility and chemical potential gradient.
$$ J_k = -\Gamma_k^* \frac{\partial \mu_k^{eq}}{\partial z} $$
$$ \Gamma_k^\phi = M_k^\phi x_k^\phi $$
$$ \Gamma_k^* = f(\Gamma_k^\alpha, \Gamma_k^\beta, …) $$
$\Gamma_k^*$ is an average mobility term that assumes certain geometry in the system. The following averaging functions are available in kawin:
- Upper Wiener - assumes phases are continuous layers parallel to flux
- Lower Wiener - assumes phases are continuous layers orthogonal to flux
- Upper Hashin-Shtrikman - assumes a matrix of the phase with the fastest mobility with spheres of all other phases
- Lower Hashin-Shtrikman - assumes a matrix of the phase with the slowest mobility with spheres of all other phases
- Labyrinth - assumes phases as precipitates
Note that the Hashin-Shtrikman bounds are much narrower than the Wiener bounds.
The fluxes are calculated in a lattice fixed frame of reference. To convert to a volume fixed frame, the flux is then defined by:
$$ J_k^v = J_k - x_k \sum{J_j} $$
In this example a Fe-25.7Cr-6.5Ni / Fe-42.3Cr-27.6Ni diffusion couple will be simulated using the lower and upper Hashin-Shtrikman bounds. Both sides of the diffusion couple are $\alpha+\gamma$.
The first step is the load the thermodynamic database.
1 import matplotlib.pyplot as plt
2 from kawin.thermo import GeneralThermodynamics
3 from kawin.diffusion import HomogenizationModel
4 from kawin.solver import explicitEulerIterator
5
6 elements = ['FE', 'CR', 'NI']
7 phases = ['FCC_A1', 'BCC_A2']
8
9 therm = GeneralThermodynamics('FeCrNi.tdb', elements, phases)
10 from kawin.diffusion import HomogenizationParameters
11 from kawin.diffusion.mesh import Cartesian1D, ProfileBuilder, StepProfile1D
12
13 profile = ProfileBuilder()
14 profile.addBuildStep(StepProfile1D(0, [0.257, 0.065], [0.423, 0.276]), ['CR', 'NI'])
15 mesh = Cartesian1D(['CR', 'NI'], [-5e-4, 5e-4], 200)
16 mesh.setResponseProfile(profile)
17
18 temperature = 1100+273.15
19
20 mobLower = HomogenizationParameters(homogenizationFunction='hashin lower', eps=0.01)
21 mobUpper = HomogenizationParameters(homogenizationFunction='hashin upper', eps=0.01)
Defining the homogenization model is similar to defining the single phase diffusion model where the bounds of the domain, the number of volume elements, the defined elements and the defined phases are needed.
As with the single phase diffusion model, inputting the composition profile and parameters are also the same. The only difference is that two extra parameters will be defined for the homogenization model:
Smoothing factor ($\varepsilon$) - this factor allows for the composition to smooth out when the chemical potential gradient is zero but the composition gradient is non-zero (in n-phase regions where n is the number of components). This can be viewed as an ideal contribution where the composition smoothes out to maximize entropy. By default, it is set to 0.05, but here, we will set it to 0.01.
Mobility function - this defined which of the above mentioned mobility functions to use. We will start with the lower Hashin-Shtrikman bounds.
Solving the model is also similar to the single phase diffusion model.
22 ml = HomogenizationModel(mesh, elements, phases,
23 thermodynamics=therm,
24 temperature=temperature,
25 homogenizationParameters=mobLower, record=False)
26 ml.solve(100*3600, iterator=explicitEulerIterator, verbose=True, vIt=500)
Iteration Sim Time (h) Run time (s)
0 0.0e+00 0.0
240 1.0e+02 20.8
The next model will be the exact same except the mobility function will be switched to the upper Hashin-Shtrikman bounds.
27 mu = HomogenizationModel(mesh, elements, phases,
28 thermodynamics=therm,
29 temperature=temperature,
30 homogenizationParameters=mobUpper, record=False)
31 mu.solve(100*3600, iterator=explicitEulerIterator, verbose=True, vIt=500)
Iteration Sim Time (h) Run time (s)
0 0.0e+00 0.0
500 1.4e+01 35.8
1000 3.0e+01 68.0
1500 4.5e+01 90.3
2000 6.0e+01 106.2
2500 7.4e+01 120.0
3000 8.8e+01 132.6
3405 1.0e+02 143.0
To compare the two mobility functions, the Cr composition, Ni composition and $\alpha$ phase fraction profile will be plotted. By default, the plotting functions will plot all components or phases; however, an individual component or phase can be defined to have it be the only thing that is plotted.
Here, we can see that the upper Hashin-Shtrikman bounds gives a smoother Cr and Ni profile. Additionally, the lower Hashin-Shtrikman bounds shows a pure $\gamma$ layer near the interface of around 4-6 $\mu m$.
32 from kawin.diffusion.Plot import plot1D, plot1DPhases
33
34 fig, ax = plt.subplots(1,3, figsize=(15,4))
35
36 plot1D(ml, 'CR', zScale=1e-3, ax=ax[0], label='lower')
37 plot1D(mu, 'CR', zScale=1e-3, ax=ax[0], label='upper')
38 ax[0].set_xlabel('Distance (mm)')
39 ax[0].set_ylabel('Composition Cr (%at)')
40 ax[0].set_ylim([0.2, 0.45])
41 ax[0].legend()
42
43 plot1D(ml, 'NI', zScale=1e-3, ax=ax[1], label='lower')
44 plot1D(mu, 'NI', zScale=1e-3, ax=ax[1], label='upper')
45 ax[1].set_xlabel('Distance (mm)')
46 ax[1].set_ylabel('Composition Ni (%at)')
47 ax[1].set_ylim([0, 0.35])
48
49 plot1DPhases(ml, 'BCC_A2', zScale=1e-3, ax=ax[2], label='lower')
50 plot1DPhases(mu, 'BCC_A2', zScale=1e-3, ax=ax[2], label='upper')
51 ax[2].set_xlabel('Distance (mm)')
52 ax[2].set_ylabel(r'Fraction $\alpha$')
53 ax[2].set_ylim([0, 0.8])
54 plt.tight_layout()
55 plt.show()
56 from pycalphad import ternplot, variables as v
57 from pycalphad.plot import triangular
58 from pycalphad.plot.utils import phase_legend
59
60 fig = plt.figure(figsize=(6,6))
61 ax = fig.add_subplot(projection='triangular')
62
63 conds = {v.T: 1100+273.15, v.P:101325, v.X('CR'): (0,1,0.015), v.X('NI'): (0,1,0.015)}
64 ternplot(therm.db, ['FE', 'CR', 'NI', 'VA'], phases, conds, x=v.X('CR'), y=v.X('NI'), ax = ax)
65
66 # Compositions will contain all elements (dependent, independent)
67 compl = ml.getCompositions()
68 ln1, = ax.plot(compl[:,1], compl[:,2], label='lower')
69
70 compu = mu.getCompositions()
71 ln2, = ax.plot(compu[:,1], compu[:,2], label='upper')
72
73 #The pycalphad ternplot function will automatically add a legend for the phases,
74 #but the legend has to be added again to add labels for the diffusion paths
75 handles, _ = phase_legend(phases)
76 ax.legend(handles = handles + [ln1, ln2])
77
78 plt.show()
Periodic boundaries
79 from kawin.diffusion.mesh import PeriodicBoundary1D
80
81 profile = ProfileBuilder()
82 profile.addBuildStep(StepProfile1D(0.5e-3, [0.257, 0.065], [0.423, 0.276]), ['CR', 'NI'])
83
84 mesh = Cartesian1D(['CR', 'NI'], [0, 1e-3], 200)
85 mesh.setResponseProfile(profile, boundaryConditions=PeriodicBoundary1D())
86
87 m = HomogenizationModel(mesh, elements, phases,
88 thermodynamics=therm, temperature=temperature,
89 homogenizationParameters=mobLower, record=False)
90 m.solve(100*3600, iterator=explicitEulerIterator, verbose=True, vIt=500)
91
92 fig, ax = plt.subplots(1,3, figsize=(15,4))
93
94 plot1D(m, 'CR', zScale=1e-3, ax=ax[0])
95 ax[0].set_xlabel('Distance (mm)')
96 ax[0].set_ylabel('Composition Cr (%at)')
97 ax[0].set_ylim([0.2, 0.45])
98
99 plot1D(m, 'NI', zScale=1e-3, ax=ax[1])
100 ax[1].set_xlabel('Distance (mm)')
101 ax[1].set_ylabel('Composition Ni (%at)')
102 ax[1].set_ylim([0, 0.35])
103
104 plot1DPhases(m, 'BCC_A2', zScale=1e-3, ax=ax[2])
105 ax[2].set_xlabel('Distance (mm)')
106 ax[2].set_ylabel(r'Fraction $\alpha$')
107 ax[2].set_ylim([0, 0.8])
108 plt.tight_layout()
109 plt.show()
Iteration Sim Time (h) Run time (s)
0 0.0e+00 0.0
240 1.0e+02 18.5
References
- H. Larsson and L. Hoglund, “Multiphase diffusion simulations in 1D using the DICTRA homogenization model” Calphad 33 (2009) p. 495
- H. Larsson and A. Engstrom, “A homogenization approach to diffusion simulations applied to $\alpha+\gamma$ Fe-Cr-Ni diffusion couples” Acta Materialia 54 (2006) p. 2431