Grand Potential Multi-Phase, Multi-Order Parameter Model
A phase-field model for an arbitrary number of phases, grains per phase, and chemical constituents has been implemented in MOOSE. This model is based on a functional of the grand potential density (rather than the Helmholtz free energy density as is more commonly used in phase-field modeling). The model was originally described in Aagesen et al. (2018).
The model has certain advantages and disadvantages relative to other phase-field models in the literature.
Advantages
Allows interfacial thickness and energy to be set independently, enabling coarser mesh, improved computational performance
Chemical free energy contribution is removed from interfacial energy
Similar to Kim-Kim-Suzuki (KKS) phase-field model in this respect, but do not need separate phase concentration variables, so performance is improved
Prevents spurious formation of any additional phase ("ghost phase") at two-phase interfaces
Works well with ideal solution model/low equilibrium concentration situations compared with KKS
Disadvantages
Currently limited to parabolic, ideal solution, dilute solution chemical free energies
Evolution equation for chemical potential less intuitive compared to evolution equation for composition
When using the evolution equation for chemical potential, discretization error results in small variations in mass (typically < 0.1 %). To prevent this from occurring, an alternative version of the model that strictly enforces mass conservation is available. However, this requires an additional concentration variable, increasing computational cost.
What is the grand potential?
The grand potential is the thermodynamic potential that is minimized at constant temperature, volume, and chemical potential. For phase , the grand potential density per unit volume is given by where is the Helmholtz free energy density of phase , is the number density of atoms of solute species , and is the chemical potential of solute species .
Grand potential functional
The model is derived from a functional of the total grand potential : where is a multi-well contribution to the grand potential density, is the gradient energy contribution, and is the chemical energy contribution. is a function of the grand potential densities of each of the phases in the system: where is a switching function that interpolates the chemical contribution based on what phase the microstructure is in at each point. is given by where is a constant with dimensions of energy per unit volume and is given by where is an order parameter representing grain of phase ( index phases and index grains). is a set of dimensionless parameters that adjust grain boundary energies and interfacial energies. This function has minima at the equilibrium values of the order parameters. is given by where is the gradient energy cofficient, which has dimensions of energy per unit length.
Evolution equations
The order parameters representing the grains of each phase each evolve by an Allen-Cahn equation, derived from the grand potential functional: where is the Allen-Cahn interfacial mobility.
In the original model formulation, the evolution of composition is transformed into an evolution equation for chemical potential. The general multi-species evolution equation is given in Aagesen et al. (2018). The simplest form of the evolution equation for solute species is given by where is the self-diffusivity of species and is the susceptibility, defined as depends on the form of the chemical free energy density used.
MOOSE implementation of the evolution equations
The kernels used to implement the terms in the evolution equations are shown below with underbraces. For the Allen-Cahn equation:
For the chemical potential evolution equation,
Creating an input file for this model can be simplified through the use of the MOOSE action GrandPotentialKernelAction
, which automates the process of adding the required kernels.
Parameterization: Interface thickness and energy
For an interface between grain of phase and grain of phase , the interfacial thickness and interfacial energy are set by the combination of parameters , , and . For , analytical relationships exist: A convenient strategy for parameterization is to pick for one of the interfaces, preferably the one with the median interfacial energy of all the types of interface. The analytical relationships above can be used to calculate and . Once calculated, , , and can be set, normally using a GenericConstantMaterial
. Once and are set, the interfacial energy for other types of interface can be set using the other parameters: where is a dimensionless function of for the other types of interfaces and can be determined based on the known , and for the other interfaces. The following polynomial approximation can be used to determine as a function of : The values for for the other interfaces can be specified using a GenericConstantMaterial
. Alternatively, rather than calculating and specifying these values by hand, the material GrandPotentialInterface
can be used.
Version with strict mass conservation
In the original formulation of the model, small variations in total solute concentration may occur due to discretization error of the chemical potential evolution equation. Typically such variation is less than 0.1% of total solute in the system, and decreases with decreasing time step. The careful choice of time step can be used to lower solute variation to a level low enough that it does not affect microstructural evolution. To eliminate this variation entirely, an alternative formulation has been developed that maintains the use of an evolution equation for composition as a field variable. However, the use of the additional field variable results in increased computational cost. In this formulation, the following equation i used:
where is the atomic volume of species and is related to using . The relationship between composition and chemical potential must also be specified, such as Eq. 22 in Aagesen et al. (2018) for parabolic free energy forms:
An example input file can be found here: (../moose/modules/phase_field/examples/multiphase/GrandPotential3Phase_masscons.i)
References
- Larry K. Aagesen, Yipeng Gao, Daniel Schwen, and Karim Ahmed.
Grand-potential-based phase-field model for multiple phases, grains, and chemical components.
Phys. Rev. E, 98:023309, Aug 2018.
URL: https://link.aps.org/doi/10.1103/PhysRevE.98.023309, doi:10.1103/PhysRevE.98.023309.[BibTeX]
@article{AagesenGP2018,
author = "Aagesen, Larry K. and Gao, Yipeng and Schwen, Daniel and Ahmed, Karim",
title = "Grand-potential-based phase-field model for multiple phases, grains, and chemical components",
journal = "Phys. Rev. E",
volume = "98",
issue = "2",
pages = "023309",
numpages = "16",
year = "2018",
month = "Aug",
publisher = "American Physical Society",
doi = "10.1103/PhysRevE.98.023309",
url = "https://link.aps.org/doi/10.1103/PhysRevE.98.023309"
}
(../moose/modules/phase_field/test/tests/GrandPotentialPFM/GrandPotentialMultiphase.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 20
ny = 20
xmin = -20
xmax = 20
ymin = -20
ymax = 20
[]
[GlobalParams]
op_num = 2
var_name_base = etab
[]
[Variables]
[./w]
[../]
[./etaa0]
[../]
[./etab0]
[../]
[./etab1]
[../]
[]
[AuxVariables]
[./bnds]
order = FIRST
family = LAGRANGE
[../]
[]
[ICs]
[./IC_etaa0]
type = FunctionIC
variable = etaa0
function = ic_func_etaa0
[../]
[./IC_etab0]
type = FunctionIC
variable = etab0
function = ic_func_etab0
[../]
[./IC_etab1]
type = FunctionIC
variable = etab1
function = ic_func_etab1
[../]
[./IC_w]
type = ConstantIC
value = -0.05
variable = w
[../]
[]
[Functions]
[./ic_func_etaa0]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2);0.5*(1.0-tanh((r-10.0)/sqrt(2.0)))'
[../]
[./ic_func_etab0]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2);0.5*(1.0+tanh((r-10)/sqrt(2.0)))*0.5*(1.0+tanh((y)/sqrt(2.0)))'
[../]
[./ic_func_etab1]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2);0.5*(1.0+tanh((r-10)/sqrt(2.0)))*0.5*(1.0-tanh((y)/sqrt(2.0)))'
[../]
[]
[BCs]
[]
[Kernels]
# Order parameter eta_alpha0
[./ACa0_bulk]
type = ACGrGrMulti
variable = etaa0
v = 'etab0 etab1'
gamma_names = 'gab gab'
[../]
[./ACa0_sw]
type = ACSwitching
variable = etaa0
Fj_names = 'omegaa omegab'
hj_names = 'ha hb'
coupled_variables = 'etab0 etab1 w'
[../]
[./ACa0_int]
type = ACInterface
variable = etaa0
kappa_name = kappa
[../]
[./ea0_dot]
type = TimeDerivative
variable = etaa0
[../]
# Order parameter eta_beta0
[./ACb0_bulk]
type = ACGrGrMulti
variable = etab0
v = 'etaa0 etab1'
gamma_names = 'gab gbb'
[../]
[./ACb0_sw]
type = ACSwitching
variable = etab0
Fj_names = 'omegaa omegab'
hj_names = 'ha hb'
coupled_variables = 'etaa0 etab1 w'
[../]
[./ACb0_int]
type = ACInterface
variable = etab0
kappa_name = kappa
[../]
[./eb0_dot]
type = TimeDerivative
variable = etab0
[../]
# Order parameter eta_beta1
[./ACb1_bulk]
type = ACGrGrMulti
variable = etab1
v = 'etaa0 etab0'
gamma_names = 'gab gbb'
[../]
[./ACb1_sw]
type = ACSwitching
variable = etab1
Fj_names = 'omegaa omegab'
hj_names = 'ha hb'
coupled_variables = 'etaa0 etab0 w'
[../]
[./ACb1_int]
type = ACInterface
variable = etab1
kappa_name = kappa
[../]
[./eb1_dot]
type = TimeDerivative
variable = etab1
[../]
#Chemical potential
[./w_dot]
type = SusceptibilityTimeDerivative
variable = w
f_name = chi
coupled_variables = '' # in this case chi (the susceptibility) is simply a constant
[../]
[./Diffusion]
type = MatDiffusion
variable = w
diffusivity = Dchi
args = ''
[../]
[./coupled_etaa0dot]
type = CoupledSwitchingTimeDerivative
variable = w
v = etaa0
Fj_names = 'rhoa rhob'
hj_names = 'ha hb'
coupled_variables = 'etaa0 etab0 etab1'
[../]
[./coupled_etab0dot]
type = CoupledSwitchingTimeDerivative
variable = w
v = etab0
Fj_names = 'rhoa rhob'
hj_names = 'ha hb'
coupled_variables = 'etaa0 etab0 etab1'
[../]
[./coupled_etab1dot]
type = CoupledSwitchingTimeDerivative
variable = w
v = etab1
Fj_names = 'rhoa rhob'
hj_names = 'ha hb'
coupled_variables = 'etaa0 etab0 etab1'
[../]
[]
[AuxKernels]
[./BndsCalc]
type = BndsCalcAux
variable = bnds
execute_on = timestep_end
[../]
[]
# enable_jit set to false in many materials to make this test start up faster.
# It is recommended to set enable_jit = true or just remove these lines for
# production runs with this model
[Materials]
[./ha]
type = SwitchingFunctionMultiPhaseMaterial
h_name = ha
all_etas = 'etaa0 etab0 etab1'
phase_etas = 'etaa0'
[../]
[./hb]
type = SwitchingFunctionMultiPhaseMaterial
h_name = hb
all_etas = 'etaa0 etab0 etab1'
phase_etas = 'etab0 etab1'
[../]
[./omegaa]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = omegaa
material_property_names = 'Vm ka caeq'
expression = '-0.5*w^2/Vm^2/ka-w/Vm*caeq'
derivative_order = 2
enable_jit = false
[../]
[./omegab]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = omegab
material_property_names = 'Vm kb cbeq'
expression = '-0.5*w^2/Vm^2/kb-w/Vm*cbeq'
derivative_order = 2
enable_jit = false
[../]
[./rhoa]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = rhoa
material_property_names = 'Vm ka caeq'
expression = 'w/Vm^2/ka + caeq/Vm'
derivative_order = 2
enable_jit = false
[../]
[./rhob]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = rhob
material_property_names = 'Vm kb cbeq'
expression = 'w/Vm^2/kb + cbeq/Vm'
derivative_order = 2
enable_jit = false
[../]
[./const]
type = GenericConstantMaterial
prop_names = 'kappa_c kappa L D chi Vm ka caeq kb cbeq gab gbb mu'
prop_values = '0 1 1.0 1.0 1.0 1.0 10.0 0.1 10.0 0.9 4.5 1.5 1.0'
[../]
[./Mobility]
type = DerivativeParsedMaterial
property_name = Dchi
material_property_names = 'D chi'
expression = 'D*chi'
derivative_order = 2
enable_jit = false
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = bdf2
petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_pc_type -pc_asm_overlap'
petsc_options_value = 'asm 31 lu 1'
l_tol = 1.0e-3
nl_rel_tol = 1.0e-8
nl_abs_tol = 1e-8
num_steps = 2
[./TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 0.1
[../]
[]
[Outputs]
exodus = true
[]
(../moose/modules/phase_field/examples/multiphase/GrandPotential3Phase_masscons.i)
# This is an example of implementation of the multi-phase, multi-order parameter
# grand potential based phase-field model described in Phys. Rev. E, 98, 023309
# (2019). It includes 3 phases with 1 grain of each phase.
# This is a revised version of the model that eliminates small variations in mass
# that have been observed with the original formulation. In this version, rather
# than evolving the chemical potential as a field variable, we evolve the composition
# field using a normal Cahn-Hilliard equation, then relate chemical potential to
# composition using Eq. (22) from the paper (this relationship is derived from the
# grand potential functional and is valid only for parabolic free energies).
[Mesh]
type = GeneratedMesh
dim = 2
nx = 60
ny = 60
xmin = -15
xmax = 15
ymin = -15
ymax = 15
[]
[Variables]
[w]
[]
[c]
[]
[etaa0]
[]
[etab0]
[]
[etad0]
[]
[]
[ICs]
[IC_etaa0]
type = BoundingBoxIC
variable = etaa0
x1 = -10
y1 = -10
x2 = 10
y2 = 10
inside = 1.0
outside = 0.0
[]
[IC_etad0]
type = BoundingBoxIC
variable = etad0
x1 = -10
y1 = -10
x2 = 10
y2 = 10
inside = 0.0
outside = 1.0
[]
[IC_c]
type = BoundingBoxIC
variable = c
x1 = -10
y1 = -10
x2 = 10
y2 = 10
inside = 0.1
outside = 0.5
[]
[IC_w]
type = FunctionIC
variable = w
function = ic_func_w
[]
[]
[Functions]
[ic_func_w]
type = ConstantFunction
value = 0
[]
[]
[Kernels]
# Order parameter eta_alpha0
[ACa0_bulk]
type = ACGrGrMulti
variable = etaa0
v = 'etab0 etad0'
gamma_names = 'gab gad'
[]
[ACa0_sw]
type = ACSwitching
variable = etaa0
Fj_names = 'omegaa omegab omegad'
hj_names = 'ha hb hd'
coupled_variables = 'etab0 etad0 w'
[]
[ACa0_int]
type = ACInterface
variable = etaa0
kappa_name = kappa
[]
[ea0_dot]
type = TimeDerivative
variable = etaa0
[]
# Order parameter eta_beta0
[ACb0_bulk]
type = ACGrGrMulti
variable = etab0
v = 'etaa0 etad0'
gamma_names = 'gab gbd'
[]
[ACb0_sw]
type = ACSwitching
variable = etab0
Fj_names = 'omegaa omegab omegad'
hj_names = 'ha hb hd'
coupled_variables = 'etaa0 etad0 w'
[]
[ACb0_int]
type = ACInterface
variable = etab0
kappa_name = kappa
[]
[eb0_dot]
type = TimeDerivative
variable = etab0
[]
# Order parameter eta_delta0
[ACd0_bulk]
type = ACGrGrMulti
variable = etad0
v = 'etaa0 etab0'
gamma_names = 'gad gbd'
[]
[ACd0_sw]
type = ACSwitching
variable = etad0
Fj_names = 'omegaa omegab omegad'
hj_names = 'ha hb hd'
coupled_variables = 'etaa0 etab0 w'
[]
[ACd0_int]
type = ACInterface
variable = etad0
kappa_name = kappa
[]
[ed0_dot]
type = TimeDerivative
variable = etad0
[]
#Concentration
[c_dot]
type = TimeDerivative
variable = c
[]
[Diffusion]
type = MatDiffusion
variable = c
v = w
diffusivity = DchiVm
args = ''
[]
#The following relate chemical potential to composition using Eq. (22)
[w_rxn]
type = MatReaction
variable = w
v = c
mob_name = -1
[]
[ca_rxn]
type = MatReaction
variable = w
mob_name = 'hoverk_a'
args = 'etaa0 etab0 etad0'
[]
[ca_bodyforce]
type = MaskedBodyForce
variable = w
mask = ha
coupled_variables = 'etaa0 etab0 etad0'
value = 0.1 #caeq
[]
[cb_rxn]
type = MatReaction
variable = w
mob_name = 'hoverk_b'
args = 'etaa0 etab0 etad0'
[]
[cb_bodyforce]
type = MaskedBodyForce
variable = w
mask = hb
coupled_variables = 'etaa0 etab0 etad0'
value = 0.9 #cbeq
[]
[cd_rxn]
type = MatReaction
variable = w
mob_name = 'hoverk_d'
args = 'etaa0 etab0 etad0'
[]
[cd_bodyforce]
type = MaskedBodyForce
variable = w
mask = hd
coupled_variables = 'etaa0 etab0 etad0'
value = 0.5 #cdeq
[]
[]
[Materials]
[ha_test]
type = SwitchingFunctionMultiPhaseMaterial
h_name = ha
all_etas = 'etaa0 etab0 etad0'
phase_etas = 'etaa0'
[]
[hb_test]
type = SwitchingFunctionMultiPhaseMaterial
h_name = hb
all_etas = 'etaa0 etab0 etad0'
phase_etas = 'etab0'
[]
[hd_test]
type = SwitchingFunctionMultiPhaseMaterial
h_name = hd
all_etas = 'etaa0 etab0 etad0'
phase_etas = 'etad0'
[]
[omegaa]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = omegaa
material_property_names = 'Vm ka caeq'
expression = '-0.5*w^2/Vm^2/ka-w/Vm*caeq'
derivative_order = 2
[]
[omegab]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = omegab
material_property_names = 'Vm kb cbeq'
expression = '-0.5*w^2/Vm^2/kb-w/Vm*cbeq'
derivative_order = 2
[]
[omegad]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = omegad
material_property_names = 'Vm kd cdeq'
expression = '-0.5*w^2/Vm^2/kd-w/Vm*cdeq'
derivative_order = 2
[]
[rhoa]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = rhoa
material_property_names = 'Vm ka caeq'
expression = 'w/Vm^2/ka + caeq/Vm'
derivative_order = 2
[]
[rhob]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = rhob
material_property_names = 'Vm kb cbeq'
expression = 'w/Vm^2/kb + cbeq/Vm'
derivative_order = 2
[]
[rhod]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = rhod
material_property_names = 'Vm kd cdeq'
expression = 'w/Vm^2/kd + cdeq/Vm'
derivative_order = 2
[]
[const]
type = GenericConstantMaterial
prop_names = 'kappa_c kappa L D Vm ka caeq kb cbeq kd cdeq gab gad gbd mu tgrad_corr_mult'
prop_values = '0 1 1.0 1.0 1.0 10.0 0.1 10.0 0.9 10.0 0.5 1.5 1.5 1.5 1.0 0.0'
[]
[Mobility]
type = DerivativeParsedMaterial
property_name = DchiVm
material_property_names = 'D chi Vm' #Factor of Vm is needed to evolve c instead of rho
expression = 'D*chi*Vm'
derivative_order = 2
[]
[chi]
type = DerivativeParsedMaterial
property_name = chi
material_property_names = 'Vm ha(etaa0,etab0,etad0) ka hb(etaa0,etab0,etad0) kb hd(etaa0,etab0,etad0) kd'
expression = '(ha/ka + hb/kb + hd/kd) / Vm^2'
coupled_variables = 'etaa0 etab0 etad0'
derivative_order = 2
[]
[hoverk_a]
type = DerivativeParsedMaterial
material_property_names = 'ha(etaa0,etab0,etad0) Vm ka'
property_name = hoverk_a
expression = 'ha / Vm / ka'
[]
[hoverk_b]
type = DerivativeParsedMaterial
material_property_names = 'hb(etaa0,etab0,etad0) Vm kb'
property_name = hoverk_b
expression = 'hb / Vm / kb'
[]
[hoverk_d]
type = DerivativeParsedMaterial
material_property_names = 'hd(etaa0,etab0,etad0) Vm kd'
property_name = hoverk_d
expression = 'hd / Vm / kd'
[]
[]
[Postprocessors]
[c_total]
type = ElementIntegralVariablePostprocessor
variable = c
[]
[]
[Executioner]
type = Transient
nl_max_its = 15
scheme = bdf2
solve_type = NEWTON
petsc_options_iname = -pc_type
petsc_options_value = asm
l_max_its = 15
l_tol = 1.0e-3
nl_rel_tol = 1.0e-8
start_time = 0.0
num_steps = 20
nl_abs_tol = 1e-10
dt = 1.0
[]
[Outputs]
csv = true
exodus = true
[]