The precipitation simulation of Ni4Ti3 phase

[solo2025]

We tried to reproduce the precipitation results of Ni4ti3 in Section 3.32 of the paper (https://doi.org/10.1016/j.actamat.2023.118988) to test the effect of adding elastic energy on precipitation. The parameters are completely based on the reference given in Figure 13 of the paper, but the calculation does not converge. Can you provide a script that can be run for this case? Attached below is our script.

from kawin.thermo import BinaryThermodynamics
from kawin.precipitation import PrecipitateModel,VolumeParameter, ShapeFactor
from kawin.precipitation.non_ideal.ElasticFactors import StrainEnergy
from kawin.precipitation.non_ideal.ShapeFactors import ShapeFactor

import numpy as np
import matplotlib.pyplot as plt
import matplotlib ``
matplotlib.use('TkAgg')


# initialize model
model = PrecipitateModel(phases=['TI3NI4'], elements=['TI'])

# set up thermodynamics
therm = BinaryThermodynamics(
    'niti-2013calphad.tdb', 
    ['NI', 'TI'], 
    ['BCC_B2', 'TI3NI4'],
    interfacialCompMethod='equilibrium',
    drivingForceMethod='approximate'
    )
# therm.setMobilityCorrection('all', 100)

# Override guess composition(s) to reduce number of calculations
# when finding interfacial composition
therm.setGuessComposition(0.56)


# Model parameters
x_init = 0.508
gamma = 0.053
T = 450 + 273.15
Dni  = lambda x, T: 1.8E-8 * np.exp(-155000/(8.314*T))
vaBCC, nBCC = 0.0268114E-27, 2
vaTI3NI4, nTI3NI4 = 0.184615E-27, 14

#Initialize model
model.setInitialComposition(x_init)
model.setTemperature(T)
model.setInterfacialEnergy(gamma)
model.setDiffusivity(Dni)
model.setVolumeAlpha(vaBCC,VolumeParameter.MOLAR_VOLUME, nBCC)
model.setVolumeBeta(vaTI3NI4, VolumeParameter.MOLAR_VOLUME, nTI3NI4)
model.setThermodynamics(therm)
model.setNucleationSite('bulk')
model.setNucleationDensity(bulkN0=1e30)


se = StrainEnergy()

# strain energy parameters
B2e = np.asarray([175,45,35]) * 1e9
eigenstrain = [-0.00417, -0.00417, -0.0257]
rotate = [
    [-4/np.sqrt(42), 5/np.sqrt(42), -1/np.sqrt(42)],
    [-2/np.sqrt(14), -1/np.sqrt(14), 3/np.sqrt(14)], 
    [1/np.sqrt(3), 1/np.sqrt(3), 1/np.sqrt(3)]]

# se.setup()
se.setEigenstrain(eigenstrain)
se.setElasticConstants(B2e [0],B2e [1],B2e [2])
se.setRotationMatrix(rotate)
se.setEllipsoidal()

model.setStrainEnergy(se, calculateAspectRatio=True)
model.setPrecipitateShape(ShapeFactor.PLATE)

# model.solve(1e1, verbose=True, vIt=2000)
model.solve(1e2, verbose=True, vIt=10)


import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(12, 8))

model.plot(axes[0,0], 'Precipitate Density', bounds=[1e-2, 1e4], timeUnits='min')
axes[0,0].set_ylim([1e10, 1e28])
axes[0,0].set_yscale('log')

model.plot(axes[0,1], 'Average Radius', bounds=[1e-2, 1e4], timeUnits='min', label='Radius')
axes[0,1].set_ylim([0, 7e-9])

ax1 = axes[0,1].twinx()
model.plot(ax1, 'Aspect Ratio', bounds=[1e-2, 1e4], timeUnits='min', label='Aspect Ratio', linestyle=':')
ax1.set_ylim([1,4])

plt.tight_layout()
plt.show()

[nury12n]

I had to make some small changes to your script to get it working. You can see if these changes work on your end.

Main things were:

• When you define the Thermodynamics class, it assumes that the first element is the solvent and the remaining elements are solutes. For the Ni-Ti system here, TI is the solvent and NI is the solute, so I switched the order of the elements.
  • This also relates to the guess composition since Ti3Ni4 has a composition of 57.1 at% NI while the BCC_B2 phase at equilibrium is near 50 at% NI, so the guess composition is set between these two values. (so if you had the elements as ['NI', 'TI'], then it would try to compute equilibrium at 57.1 at% TI (42.9 at% NI), where Ti3Ni4 would be unstable, so an interfacial composition or driving force would not be found).
• The volume of the two phases are in atomic volume, so the volume parameter is set to VolumeParameter.ATOMIC_VOLUME
• When adding the thermodynamics class to the model, an extra parameter addDiffusivity is set to False. Unfortunately, this is not clearly stated in the documentation, but when you add the thermodynamics to the model, it will attempt to add diffusivity or mobility models that the thermodynamics class found from the TDB file. In this case, we manually set the diffusivity, so we want the model to not check whether it's defined in the thermodynamics class
• Minor changes to the plotting functions at the end

from kawin.thermo import BinaryThermodynamics
from kawin.precipitation import PrecipitateModel,VolumeParameter, ShapeFactor
from kawin.precipitation.non_ideal.ElasticFactors import StrainEnergy
from kawin.precipitation.non_ideal.ShapeFactors import ShapeFactor

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('TkAgg')


# initialize model
model = PrecipitateModel(phases=['TI3NI4'], elements=['NI'])

# set up thermodynamics
therm = BinaryThermodynamics(
    'niti-2013calphad.tdb', 
    ['TI', 'NI'], 
    ['BCC_B2', 'TI3NI4'],
    interfacialCompMethod='equilibrium',
    drivingForceMethod='approximate'
    )
# therm.setMobilityCorrection('all', 100)

# Override guess composition(s) to reduce number of calculations
# when finding interfacial composition
therm.setGuessComposition(0.56)


# Model parameters
x_init = 0.508
gamma = 0.053
T = 450 + 273.15
Dni  = lambda x, T: 1.8E-8 * np.exp(-155000/(8.314*T))
vaBCC, nBCC = 0.0268114E-27, 2
vaTI3NI4, nTI3NI4 = 0.184615E-27, 14

#Initialize model
model.setInitialComposition(x_init)
model.setTemperature(T)
model.setInterfacialEnergy(gamma)
model.setDiffusivity(Dni)
model.setVolumeAlpha(vaBCC,VolumeParameter.ATOMIC_VOLUME, nBCC)
model.setVolumeBeta(vaTI3NI4, VolumeParameter.ATOMIC_VOLUME, nTI3NI4)
model.setThermodynamics(therm, addDiffusivity=False)
model.setNucleationSite('bulk')
model.setNucleationDensity(bulkN0=1e30)


se = StrainEnergy()

# strain energy parameters
B2e = np.asarray([175,45,35]) * 1e9
eigenstrain = [-0.00417, -0.00417, -0.0257]
rotate = [
    [-4/np.sqrt(42), 5/np.sqrt(42), -1/np.sqrt(42)],
    [-2/np.sqrt(14), -1/np.sqrt(14), 3/np.sqrt(14)], 
    [1/np.sqrt(3), 1/np.sqrt(3), 1/np.sqrt(3)]]

# se.setup()
se.setEigenstrain(eigenstrain)
se.setElasticConstants(B2e [0],B2e [1],B2e [2])
se.setRotationMatrix(rotate)
se.setEllipsoidal()

model.setStrainEnergy(se, calculateAspectRatio=True)
model.setPrecipitateShape(ShapeFactor.PLATE)

# model.solve(1e1, verbose=True, vIt=2000)
#model.solve(1e2, verbose=True, vIt=10)
model.solve(1e8, verbose=True, vIt=1000)


import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(12, 8))

model.plot(axes[0], 'Precipitate Density', bounds=[1e-2, 1e7], timeUnits='min')
axes[0].set_ylim([1e10, 1e28])
axes[0].set_yscale('log')

model.plot(axes[1], 'Average Radius', bounds=[1e-2, 1e7], timeUnits='min', label='Radius')
axes[1].set_ylim([0, 7e-8])

ax1 = axes[1].twinx()
model.plot(ax1, 'Aspect Ratio', bounds=[1e-2, 1e7], timeUnits='min', label='Aspect Ratio', linestyle=':')
ax1.set_ylim([1,4])

plt.tight_layout()
plt.show()