- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
- constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
- declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
- etaOrder parameter variable
C++ Type:std::vector<VariableName>
Controllable:No
Description:Order parameter variable
- function_namehactual name for f(eta), i.e. 'h' or 'g'
Default:h
C++ Type:std::string
Controllable:No
Description:actual name for f(eta), i.e. 'h' or 'g'
- h_orderSIMPLEPolynomial order of the switching function h(eta)
Default:SIMPLE
C++ Type:MooseEnum
Controllable:No
Description:Polynomial order of the switching function h(eta)
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
SwitchingFunctionMaterial
The SwitchingFunctionMaterial has not been documented. The content listed below should be used as a starting point for documenting the class, which includes the typical automatic documentation associated with a MooseObject; however, what is contained is ultimately determined by what is necessary to make the documentation clear for users.
Helper material to provide and its derivative in one of two polynomial forms. SIMPLE: HIGH:
Overview
Example Input File Syntax
Input Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
- outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
Outputs Parameters
Input Files
- (modules/combined/examples/phase_field-mechanics/kks_mechanics_VTS.i)
- (modules/combined/examples/phase_field-mechanics/Pattern1.i)
- (modules/combined/examples/phase_field-mechanics/Nonconserved.i)
- (modules/phase_field/test/tests/MultiPhase/orderparameterfunctionmaterial.i)
- (modules/phase_field/examples/kim-kim-suzuki/kks_example_dirichlet.i)
- (modules/combined/test/tests/multiphase_mechanics/twophasestress.i)
- (modules/phase_field/examples/slkks/CrFe.i)
- (modules/phase_field/test/tests/KKS_system/auxkernel.i)
- (modules/combined/test/tests/multiphase_mechanics/multiphasestress.i)
- (modules/phase_field/test/tests/KKS_system/kks_example_nested.i)
- (modules/combined/test/tests/surface_tension_KKS/surface_tension_KKS.i)
- (modules/phase_field/test/tests/KKS_system/two_phase.i)
- (modules/phase_field/test/tests/MultiPhase/derivativetwophasematerial.i)
- (modules/phase_field/test/tests/MultiPhase/acmultiinterface_aux.i)
- (modules/phase_field/test/tests/MultiPhase/lagrangemult.i)
- (modules/phase_field/test/tests/KKS_system/kks_phase_concentration.i)
- (modules/combined/examples/publications/rapid_dev/fig7b.i)
- (modules/combined/test/tests/surface_tension_KKS/surface_tension_VDWgas.i)
- (modules/combined/examples/publications/rapid_dev/fig7a.i)
- (modules/phase_field/test/tests/KKS_system/kks_example.i)
- (modules/phase_field/test/tests/MultiPhase/asymmetriccrosstermbarrierfunction.i)
- (modules/phase_field/test/tests/MultiPhase/acmultiinterface.i)
- (modules/phase_field/test/tests/GrandPotentialPFM/GrandPotentialPFM.i)
- (modules/phase_field/examples/multiphase/DerivativeMultiPhaseMaterial.i)
- (modules/phase_field/test/tests/KKS_system/nonlinear.i)
- (modules/phase_field/examples/kim-kim-suzuki/kks_example_noflux.i)
- (modules/combined/examples/publications/rapid_dev/fig8.i)
- (modules/phase_field/test/tests/KKS_system/kks_example_split.i)
- (modules/phase_field/examples/kim-kim-suzuki/kks_example_ternary.i)
- (modules/combined/examples/phase_field-mechanics/kks_mechanics_KHS.i)
- (modules/phase_field/test/tests/KKS_system/kks_example_offset.i)
- (modules/phase_field/test/tests/MultiPhase/penalty.i)
- (modules/combined/test/tests/linear_elasticity/extra_stress.i)
- (modules/phase_field/test/tests/slkks/full_solve.i)
- (modules/phase_field/test/tests/KKS_system/lagrange_multiplier.i)
(modules/combined/examples/phase_field-mechanics/kks_mechanics_VTS.i)
# KKS phase-field model coupled with elasticity using the Voigt-Taylor scheme as
# described in L.K. Aagesen et al., Computational Materials Science, 140, 10-21 (2017)
# Original run #170329e
[Mesh]
type = GeneratedMesh
dim = 3
nx = 640
ny = 1
nz = 1
xmin = -10
xmax = 10
ymin = 0
ymax = 0.03125
zmin = 0
zmax = 0.03125
elem_type = HEX8
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# solute concentration
[./c]
order = FIRST
family = LAGRANGE
[../]
# chemical potential
[./w]
order = FIRST
family = LAGRANGE
[../]
# solute phase concentration (matrix)
[./cm]
order = FIRST
family = LAGRANGE
[../]
# solute phase concentration (precipitate)
[./cp]
order = FIRST
family = LAGRANGE
[../]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[./disp_z]
order = FIRST
family = LAGRANGE
[../]
[]
[ICs]
[./eta_ic]
variable = eta
type = FunctionIC
function = ic_func_eta
block = 0
[../]
[./c_ic]
variable = c
type = FunctionIC
function = ic_func_c
block = 0
[../]
[./w_ic]
variable = w
type = ConstantIC
value = 0.00991
block = 0
[../]
[./cm_ic]
variable = cm
type = ConstantIC
value = 0.131
block = 0
[../]
[./cp_ic]
variable = cp
type = ConstantIC
value = 0.236
block = 0
[../]
[]
[Functions]
[./ic_func_eta]
type = ParsedFunction
expression = '0.5*(1.0+tanh((x)/delta_eta/sqrt(2.0)))'
symbol_names = 'delta_eta'
symbol_values = '0.8034'
[../]
[./ic_func_c]
type = ParsedFunction
expression = '0.2388*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^3*(6*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^2-15*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))+10)+0.1338*(1-(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^3*(6*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^2-15*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))+10))'
symbol_names = 'delta'
symbol_values = '0.8034'
[../]
[./psi_eq_int]
type = ParsedFunction
expression = 'volume*psi_alpha'
symbol_names = 'volume psi_alpha'
symbol_values = 'volume psi_alpha'
[../]
[./gamma]
type = ParsedFunction
expression = '(psi_int - psi_eq_int) / dy / dz'
symbol_names = 'psi_int psi_eq_int dy dz'
symbol_values = 'psi_int psi_eq_int 0.03125 0.03125'
[../]
[]
[AuxVariables]
[./sigma11]
order = CONSTANT
family = MONOMIAL
[../]
[./sigma22]
order = CONSTANT
family = MONOMIAL
[../]
[./sigma33]
order = CONSTANT
family = MONOMIAL
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
[../]
[./e33]
order = CONSTANT
family = MONOMIAL
[../]
[./e_el11]
order = CONSTANT
family = MONOMIAL
[../]
[./e_el12]
order = CONSTANT
family = MONOMIAL
[../]
[./e_el22]
order = CONSTANT
family = MONOMIAL
[../]
[./f_el]
order = CONSTANT
family = MONOMIAL
[../]
[./eigen_strain00]
order = CONSTANT
family = MONOMIAL
[../]
[./Fglobal]
order = CONSTANT
family = MONOMIAL
[../]
[./psi]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./matl_sigma11]
type = RankTwoAux
rank_two_tensor = stress
index_i = 0
index_j = 0
variable = sigma11
[../]
[./matl_sigma22]
type = RankTwoAux
rank_two_tensor = stress
index_i = 1
index_j = 1
variable = sigma22
[../]
[./matl_sigma33]
type = RankTwoAux
rank_two_tensor = stress
index_i = 2
index_j = 2
variable = sigma33
[../]
[./matl_e11]
type = RankTwoAux
rank_two_tensor = total_strain
index_i = 0
index_j = 0
variable = e11
[../]
[./matl_e12]
type = RankTwoAux
rank_two_tensor = total_strain
index_i = 0
index_j = 1
variable = e12
[../]
[./matl_e22]
type = RankTwoAux
rank_two_tensor = total_strain
index_i = 1
index_j = 1
variable = e22
[../]
[./matl_e33]
type = RankTwoAux
rank_two_tensor = total_strain
index_i = 2
index_j = 2
variable = e33
[../]
[./f_el]
type = MaterialRealAux
variable = f_el
property = f_el_mat
execute_on = timestep_end
[../]
[./GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fm
fb_name = fp
w = 0.0264
kappa_names = kappa
interfacial_vars = eta
[../]
[./psi_potential]
variable = psi
type = ParsedAux
coupled_variables = 'Fglobal w c f_el sigma11 e11'
expression = 'Fglobal - w*c + f_el - sigma11*e11'
[../]
[]
[BCs]
[./left_x]
type = DirichletBC
variable = disp_x
boundary = left
value = 0
[../]
[./right_x]
type = DirichletBC
variable = disp_x
boundary = right
value = 0
[../]
[./front_y]
type = DirichletBC
variable = disp_y
boundary = front
value = 0
[../]
[./back_y]
type = DirichletBC
variable = disp_y
boundary = back
value = 0
[../]
[./top_z]
type = DirichletBC
variable = disp_z
boundary = top
value = 0
[../]
[./bottom_z]
type = DirichletBC
variable = disp_z
boundary = bottom
value = 0
[../]
[]
[Materials]
# Chemical free energy of the matrix
[./fm]
type = DerivativeParsedMaterial
property_name = fm
coupled_variables = 'cm'
expression = '6.55*(cm-0.13)^2'
[../]
# Elastic energy of the matrix
[./elastic_free_energy_m]
type = ElasticEnergyMaterial
base_name = matrix
f_name = fe_m
args = ' '
outputs = exodus
[../]
# Total free energy of the matrix
[./Total_energy_matrix]
type = DerivativeSumMaterial
property_name = f_total_matrix
sum_materials = 'fm fe_m'
coupled_variables = 'cm'
[../]
# Free energy of the precipitate phase
[./fp]
type = DerivativeParsedMaterial
property_name = fp
coupled_variables = 'cp'
expression = '6.55*(cp-0.235)^2'
[../]
# Elastic energy of the precipitate
[./elastic_free_energy_p]
type = ElasticEnergyMaterial
base_name = ppt
f_name = fe_p
args = ' '
outputs = exodus
[../]
# Total free energy of the precipitate
[./Total_energy_ppt]
type = DerivativeSumMaterial
property_name = f_total_ppt
sum_materials = 'fp fe_p'
coupled_variables = 'cp'
[../]
# Total elastic energy
[./Total_elastic_energy]
type = DerivativeTwoPhaseMaterial
eta = eta
f_name = f_el_mat
fa_name = fe_m
fb_name = fe_p
outputs = exodus
W = 0
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L kappa misfit'
prop_values = '0.7 0.7 0.01704 0.00377'
[../]
#Mechanical properties
[./Stiffness_matrix]
type = ComputeElasticityTensor
C_ijkl = '103.3 74.25 74.25 103.3 74.25 103.3 46.75 46.75 46.75'
base_name = matrix
fill_method = symmetric9
[../]
[./Stiffness_ppt]
type = ComputeElasticityTensor
C_ijkl = '100.7 71.45 71.45 100.7 71.45 100.7 50.10 50.10 50.10'
base_name = ppt
fill_method = symmetric9
[../]
[./stress_matrix]
type = ComputeLinearElasticStress
base_name = matrix
[../]
[./stress_ppt]
type = ComputeLinearElasticStress
base_name = ppt
[../]
[./strain_matrix]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
base_name = matrix
[../]
[./strain_ppt]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
base_name = ppt
eigenstrain_names = 'eigenstrain_ppt'
[../]
[./eigen_strain]
type = ComputeEigenstrain
base_name = ppt
eigen_base = '1 1 1 0 0 0'
prefactor = misfit
eigenstrain_name = 'eigenstrain_ppt'
[../]
[./global_stress]
type = TwoPhaseStressMaterial
base_A = matrix
base_B = ppt
[../]
[./global_strain]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
[../]
[]
[Kernels]
[./TensorMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
# enforce c = (1-h(eta))*cm + h(eta)*cp
[./PhaseConc]
type = KKSPhaseConcentration
ca = cm
variable = cp
c = c
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotVacancies]
type = KKSPhaseChemicalPotential
variable = cm
cb = cp
fa_name = f_total_matrix
fb_name = f_total_ppt
[../]
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = KKSSplitCHCRes
variable = c
ca = cm
fa_name = f_total_matrix
w = w
[../]
[./dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[./ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = f_total_matrix
fb_name = f_total_ppt
w = 0.0264
args = 'cp cm'
[../]
[./ACBulkC]
type = KKSACBulkC
variable = eta
ca = cm
cb = cp
fa_name = f_total_matrix
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'asm ilu nonzero'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-8
nl_abs_tol = 1.0e-11
num_steps = 200
[./TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 0.5
[../]
[]
[VectorPostprocessors]
#[./eta]
# type = LineValueSampler
# start_point = '-10 0 0'
# end_point = '10 0 0'
# variable = eta
# num_points = 321
# sort_by = id
#[../]
#[./eta_position]
# type = FindValueOnLineSample
# vectorpostprocessor = eta
# variable_name = eta
# search_value = 0.5
#[../]
# [./f_el]
# type = LineMaterialRealSampler
# start = '-20 0 0'
# end = '20 0 0'
# sort_by = id
# property = f_el
# [../]
# [./f_el_a]
# type = LineMaterialRealSampler
# start = '-20 0 0'
# end = '20 0 0'
# sort_by = id
# property = fe_m
# [../]
# [./f_el_b]
# type = LineMaterialRealSampler
# start = '-20 0 0'
# end = '20 0 0'
# sort_by = id
# property = fe_p
# [../]
# [./h_out]
# type = LineMaterialRealSampler
# start = '-20 0 0'
# end = '20 0 0'
# sort_by = id
# property = h
# [../]
# [./fm_out]
# type = LineMaterialRealSampler
# start = '-20 0 0'
# end = '20 0 0'
# sort_by = id
# property = fm
# [../]
[]
[Postprocessors]
[./f_el_int]
type = ElementIntegralMaterialProperty
mat_prop = f_el_mat
[../]
[./c_alpha]
type = SideAverageValue
boundary = left
variable = c
[../]
[./c_beta]
type = SideAverageValue
boundary = right
variable = c
[../]
[./e11_alpha]
type = SideAverageValue
boundary = left
variable = e11
[../]
[./e11_beta]
type = SideAverageValue
boundary = right
variable = e11
[../]
[./s11_alpha]
type = SideAverageValue
boundary = left
variable = sigma11
[../]
[./s22_alpha]
type = SideAverageValue
boundary = left
variable = sigma22
[../]
[./s33_alpha]
type = SideAverageValue
boundary = left
variable = sigma33
[../]
[./s11_beta]
type = SideAverageValue
boundary = right
variable = sigma11
[../]
[./s22_beta]
type = SideAverageValue
boundary = right
variable = sigma22
[../]
[./s33_beta]
type = SideAverageValue
boundary = right
variable = sigma33
[../]
[./f_el_alpha]
type = SideAverageValue
boundary = left
variable = f_el
[../]
[./f_el_beta]
type = SideAverageValue
boundary = right
variable = f_el
[../]
[./f_c_alpha]
type = SideAverageValue
boundary = left
variable = Fglobal
[../]
[./f_c_beta]
type = SideAverageValue
boundary = right
variable = Fglobal
[../]
[./chem_pot_alpha]
type = SideAverageValue
boundary = left
variable = w
[../]
[./chem_pot_beta]
type = SideAverageValue
boundary = right
variable = w
[../]
[./psi_alpha]
type = SideAverageValue
boundary = left
variable = psi
[../]
[./psi_beta]
type = SideAverageValue
boundary = right
variable = psi
[../]
[./total_energy]
type = ElementIntegralVariablePostprocessor
variable = Fglobal
[../]
# Get simulation cell size from postprocessor
[./volume]
type = ElementIntegralMaterialProperty
mat_prop = 1
[../]
[./psi_eq_int]
type = FunctionValuePostprocessor
function = psi_eq_int
[../]
[./psi_int]
type = ElementIntegralVariablePostprocessor
variable = psi
[../]
[./gamma]
type = FunctionValuePostprocessor
function = gamma
[../]
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Outputs]
[./exodus]
type = Exodus
time_step_interval = 20
[../]
[./csv]
type = CSV
execute_on = 'final'
[../]
#[./console]
# type = Console
# output_file = true
# [../]
[]
(modules/combined/examples/phase_field-mechanics/Pattern1.i)
#
# Pattern example 1
#
# Phase changes driven by a combination mechanical (elastic) and chemical
# driving forces. In this three phase system a matrix phase, an oversized and
# an undersized precipitate phase compete. The chemical free energy favors a
# phase separation into either precipitate phase. A mix of both precipitate
# emerges to balance lattice expansion and contraction.
#
# This example demonstrates the use of
# * ACMultiInterface
# * SwitchingFunctionConstraintEta and SwitchingFunctionConstraintLagrange
# * DerivativeParsedMaterial
# * ElasticEnergyMaterial
# * DerivativeMultiPhaseMaterial
# * MultiPhaseStressMaterial
# which are the components to se up a phase field model with an arbitrary number
# of phases
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 80
ny = 80
nz = 0
xmin = -20
xmax = 20
ymin = -20
ymax = 20
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[GlobalParams]
# CahnHilliard needs the third derivatives
derivative_order = 3
enable_jit = true
displacements = 'disp_x disp_y'
[]
# AuxVars to compute the free energy density for outputting
[AuxVariables]
[./local_energy]
order = CONSTANT
family = MONOMIAL
[../]
[./cross_energy]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./local_free_energy]
type = TotalFreeEnergy
variable = local_energy
interfacial_vars = 'c'
kappa_names = 'kappa_c'
additional_free_energy = cross_energy
[../]
[./cross_terms]
type = CrossTermGradientFreeEnergy
variable = cross_energy
interfacial_vars = 'eta1 eta2 eta3'
kappa_names = 'kappa11 kappa12 kappa13
kappa21 kappa22 kappa23
kappa31 kappa32 kappa33'
[../]
[]
[Variables]
# Solute concentration variable
[./c]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = 0
max = 0.8
seed = 1235
[../]
[../]
# Order parameter for the Matrix
[./eta1]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[../]
# Order parameters for the 2 different inclusion orientations
[./eta2]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[../]
[./eta3]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[../]
# Mesh displacement
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
# Lagrange-multiplier
[./lambda]
order = FIRST
family = LAGRANGE
initial_condition = 1.0
[../]
[]
[Kernels]
# Set up stress divergence kernels
[./TensorMechanics]
[../]
# Cahn-Hilliard kernels
[./c_res]
type = CahnHilliard
variable = c
f_name = F
args = 'eta1 eta2 eta3'
[../]
[./time]
type = TimeDerivative
variable = c
[../]
# Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 1
[./deta1dt]
type = TimeDerivative
variable = eta1
[../]
[./ACBulk1]
type = AllenCahn
variable = eta1
args = 'eta2 eta3 c'
mob_name = L1
f_name = F
[../]
[./ACInterface1]
type = ACMultiInterface
variable = eta1
etas = 'eta1 eta2 eta3'
mob_name = L1
kappa_names = 'kappa11 kappa12 kappa13'
[../]
[./lagrange1]
type = SwitchingFunctionConstraintEta
variable = eta1
h_name = h1
lambda = lambda
[../]
# Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 2
[./deta2dt]
type = TimeDerivative
variable = eta2
[../]
[./ACBulk2]
type = AllenCahn
variable = eta2
args = 'eta1 eta3 c'
mob_name = L2
f_name = F
[../]
[./ACInterface2]
type = ACMultiInterface
variable = eta2
etas = 'eta1 eta2 eta3'
mob_name = L2
kappa_names = 'kappa21 kappa22 kappa23'
[../]
[./lagrange2]
type = SwitchingFunctionConstraintEta
variable = eta2
h_name = h2
lambda = lambda
[../]
# Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 3
[./deta3dt]
type = TimeDerivative
variable = eta3
[../]
[./ACBulk3]
type = AllenCahn
variable = eta3
args = 'eta1 eta2 c'
mob_name = L3
f_name = F
[../]
[./ACInterface3]
type = ACMultiInterface
variable = eta3
etas = 'eta1 eta2 eta3'
mob_name = L3
kappa_names = 'kappa31 kappa32 kappa33'
[../]
[./lagrange3]
type = SwitchingFunctionConstraintEta
variable = eta3
h_name = h3
lambda = lambda
[../]
# Lagrange-multiplier constraint kernel for lambda
[./lagrange]
type = SwitchingFunctionConstraintLagrange
variable = lambda
etas = 'eta1 eta2 eta3'
h_names = 'h1 h2 h3'
epsilon = 1e-6
[../]
[]
[Materials]
# declare a few constants, such as mobilities (L,M) and interface gradient prefactors (kappa*)
[./consts]
type = GenericConstantMaterial
prop_names = 'M kappa_c L1 L2 L3 kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
prop_values = '0.2 0 1 1 1 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00 '
[../]
# We use this to output the level of constraint enforcement
# ideally it should be 0 everywhere, if the constraint is fully enforced
[./etasummat]
type = ParsedMaterial
property_name = etasum
coupled_variables = 'eta1 eta2 eta3'
material_property_names = 'h1 h2 h3'
expression = 'h1+h2+h3-1'
outputs = exodus
[../]
# This parsed material creates a single property for visualization purposes.
# It will be 0 for phase 1, -1 for phase 2, and 1 for phase 3
[./phasemap]
type = ParsedMaterial
property_name = phase
coupled_variables = 'eta2 eta3'
expression = 'if(eta3>0.5,1,0)-if(eta2>0.5,1,0)'
outputs = exodus
[../]
# matrix phase
[./elasticity_tensor_1]
type = ComputeElasticityTensor
base_name = phase1
C_ijkl = '3 3'
fill_method = symmetric_isotropic
[../]
[./strain_1]
type = ComputeSmallStrain
base_name = phase1
displacements = 'disp_x disp_y'
[../]
[./stress_1]
type = ComputeLinearElasticStress
base_name = phase1
[../]
# oversized phase
[./elasticity_tensor_2]
type = ComputeElasticityTensor
base_name = phase2
C_ijkl = '7 7'
fill_method = symmetric_isotropic
[../]
[./strain_2]
type = ComputeSmallStrain
base_name = phase2
displacements = 'disp_x disp_y'
eigenstrain_names = eigenstrain
[../]
[./stress_2]
type = ComputeLinearElasticStress
base_name = phase2
[../]
[./eigenstrain_2]
type = ComputeEigenstrain
base_name = phase2
eigen_base = '0.02'
eigenstrain_name = eigenstrain
[../]
# undersized phase
[./elasticity_tensor_3]
type = ComputeElasticityTensor
base_name = phase3
C_ijkl = '7 7'
fill_method = symmetric_isotropic
[../]
[./strain_3]
type = ComputeSmallStrain
base_name = phase3
displacements = 'disp_x disp_y'
eigenstrain_names = eigenstrain
[../]
[./stress_3]
type = ComputeLinearElasticStress
base_name = phase3
[../]
[./eigenstrain_3]
type = ComputeEigenstrain
base_name = phase3
eigen_base = '-0.05'
eigenstrain_name = eigenstrain
[../]
# switching functions
[./switching1]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
h_order = SIMPLE
[../]
[./switching2]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
h_order = SIMPLE
[../]
[./switching3]
type = SwitchingFunctionMaterial
function_name = h3
eta = eta3
h_order = SIMPLE
[../]
[./barrier]
type = MultiBarrierFunctionMaterial
etas = 'eta1 eta2 eta3'
[../]
# chemical free energies
[./chemical_free_energy_1]
type = DerivativeParsedMaterial
property_name = Fc1
expression = '4*c^2'
coupled_variables = 'c'
derivative_order = 2
[../]
[./chemical_free_energy_2]
type = DerivativeParsedMaterial
property_name = Fc2
expression = '(c-0.9)^2-0.4'
coupled_variables = 'c'
derivative_order = 2
[../]
[./chemical_free_energy_3]
type = DerivativeParsedMaterial
property_name = Fc3
expression = '(c-0.9)^2-0.5'
coupled_variables = 'c'
derivative_order = 2
[../]
# elastic free energies
[./elastic_free_energy_1]
type = ElasticEnergyMaterial
base_name = phase1
f_name = Fe1
derivative_order = 2
args = 'c' # should be empty
[../]
[./elastic_free_energy_2]
type = ElasticEnergyMaterial
base_name = phase2
f_name = Fe2
derivative_order = 2
args = 'c' # should be empty
[../]
[./elastic_free_energy_3]
type = ElasticEnergyMaterial
base_name = phase3
f_name = Fe3
derivative_order = 2
args = 'c' # should be empty
[../]
# phase free energies (chemical + elastic)
[./phase_free_energy_1]
type = DerivativeSumMaterial
property_name = F1
sum_materials = 'Fc1 Fe1'
coupled_variables = 'c'
derivative_order = 2
[../]
[./phase_free_energy_2]
type = DerivativeSumMaterial
property_name = F2
sum_materials = 'Fc2 Fe2'
coupled_variables = 'c'
derivative_order = 2
[../]
[./phase_free_energy_3]
type = DerivativeSumMaterial
property_name = F3
sum_materials = 'Fc3 Fe3'
coupled_variables = 'c'
derivative_order = 2
[../]
# global free energy
[./free_energy]
type = DerivativeMultiPhaseMaterial
f_name = F
fi_names = 'F1 F2 F3'
hi_names = 'h1 h2 h3'
etas = 'eta1 eta2 eta3'
coupled_variables = 'c'
W = 3
[../]
# Generate the global stress from the phase stresses
[./global_stress]
type = MultiPhaseStressMaterial
phase_base = 'phase1 phase2 phase3'
h = 'h1 h2 h3'
[../]
[]
[BCs]
# the boundary conditions on the displacement enforce periodicity
# at zero total shear and constant volume
[./bottom_y]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0
[../]
[./top_y]
type = DirichletBC
variable = disp_y
boundary = 'top'
value = 0
[../]
[./left_x]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0
[../]
[./right_x]
type = DirichletBC
variable = disp_x
boundary = 'right'
value = 0
[../]
[./Periodic]
[./disp_x]
auto_direction = 'y'
[../]
[./disp_y]
auto_direction = 'x'
[../]
# all other phase field variables are fully periodic
[./c]
auto_direction = 'x y'
[../]
[./eta1]
auto_direction = 'x y'
[../]
[./eta2]
auto_direction = 'x y'
[../]
[./eta3]
auto_direction = 'x y'
[../]
[./lambda]
auto_direction = 'x y'
[../]
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
# We monitor the total free energy and the total solute concentration (should be constant)
[Postprocessors]
[./total_free_energy]
type = ElementIntegralVariablePostprocessor
variable = local_energy
[../]
[./total_solute]
type = ElementIntegralVariablePostprocessor
variable = c
[../]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type'
petsc_options_value = 'asm ilu'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-8
nl_abs_tol = 1.0e-10
start_time = 0.0
num_steps = 200
[./TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 0.1
[../]
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[./table]
type = CSV
delimiter = ' '
[../]
[]
[Debug]
# show_var_residual_norms = true
[]
(modules/combined/examples/phase_field-mechanics/Nonconserved.i)
#
# Example 2
# Phase change driven by a mechanical (elastic) driving force.
# An oversized phase inclusion grows under a uniaxial tensile stress.
# Check the file below for comments and suggestions for parameter modifications.
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 40
ny = 40
nz = 0
xmin = 0
xmax = 50
ymin = 0
ymax = 50
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[Variables]
[./eta]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 0
y1 = 0
radius = 30.0
invalue = 1.0
outvalue = 0.0
int_width = 10.0
[../]
[../]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[]
[Kernels]
[./TensorMechanics]
displacements = 'disp_x disp_y'
[../]
[./eta_bulk]
type = AllenCahn
variable = eta
f_name = F
[../]
[./eta_interface]
type = ACInterface
variable = eta
kappa_name = 1
[../]
[./time]
type = TimeDerivative
variable = eta
[../]
[]
#
# Try visualizing the stress tensor components as done in Conserved.i
#
[Materials]
[./consts]
type = GenericConstantMaterial
block = 0
prop_names = 'L'
prop_values = '1'
[../]
# matrix phase
[./stiffness_a]
type = ComputeElasticityTensor
base_name = phasea
block = 0
# lambda, mu values
C_ijkl = '7 7'
# Stiffness tensor is created from lambda=7, mu=7 for symmetric_isotropic fill method
fill_method = symmetric_isotropic
# See RankFourTensor.h for details on fill methods
[../]
[./strain_a]
type = ComputeSmallStrain
block = 0
displacements = 'disp_x disp_y'
base_name = phasea
[../]
[./stress_a]
type = ComputeLinearElasticStress
block = 0
base_name = phasea
[../]
[./elastic_free_energy_a]
type = ElasticEnergyMaterial
base_name = phasea
f_name = Fea
block = 0
args = ''
[../]
# oversized precipitate phase (simulated using thermal expansion)
[./stiffness_b]
type = ComputeElasticityTensor
base_name = phaseb
block = 0
# Stiffness tensor lambda, mu values
# Note that the two phases could have different stiffnesses.
# Try reducing the precipitate stiffness (to '1 1') rather than making it oversized
C_ijkl = '7 7'
fill_method = symmetric_isotropic
[../]
[./strain_b]
type = ComputeSmallStrain
block = 0
displacements = 'disp_x disp_y'
base_name = phaseb
eigenstrain_names = eigenstrain
[../]
[./eigenstrain_b]
type = ComputeEigenstrain
base_name = phaseb
eigen_base = '0.1 0.1 0.1'
eigenstrain_name = eigenstrain
[../]
[./stress_b]
type = ComputeLinearElasticStress
block = 0
base_name = phaseb
[../]
[./elastic_free_energy_b]
type = ElasticEnergyMaterial
base_name = phaseb
f_name = Feb
block = 0
args = ''
[../]
# Generate the global free energy from the phase free energies
[./switching]
type = SwitchingFunctionMaterial
block = 0
eta = eta
h_order = SIMPLE
[../]
[./barrier]
type = BarrierFunctionMaterial
block = 0
eta = eta
g_order = SIMPLE
[../]
[./free_energy]
type = DerivativeTwoPhaseMaterial
block = 0
f_name = F
fa_name = Fea
fb_name = Feb
eta = eta
args = ''
W = 0.1
derivative_order = 2
[../]
# Generate the global stress from the phase stresses
[./global_stress]
type = TwoPhaseStressMaterial
block = 0
base_A = phasea
base_B = phaseb
[../]
[]
[BCs]
[./bottom_y]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0
[../]
[./top_y]
type = DirichletBC
variable = disp_y
boundary = 'top'
value = 5
[../]
[./left_x]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0
[../]
[]
[Preconditioning]
# active = ' '
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = bdf2
# this gives best performance on 4 cores
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type '
petsc_options_value = 'asm lu'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-8
nl_abs_tol = 1.0e-10
start_time = 0.0
num_steps = 200
[./TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 0.2
[../]
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(modules/phase_field/test/tests/MultiPhase/orderparameterfunctionmaterial.i)
#
# This test validates the helper materials that generate material properties for
# the h(eta) switching function and the g(eta) double well function
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 40
ny = 5
nz = 0
xmin = 0
xmax = 1
ymin = 0
ymax = 1
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[BCs]
[./left1]
type = DirichletBC
variable = eta1
boundary = 'left'
value = 0
[../]
[./right1]
type = DirichletBC
variable = eta1
boundary = 'right'
value = 1
[../]
[./left2]
type = DirichletBC
variable = eta2
boundary = 'left'
value = 0
[../]
[./right2]
type = DirichletBC
variable = eta2
boundary = 'right'
value = 1
[../]
[]
[Variables]
# order parameter 1
[./eta1]
order = FIRST
family = LAGRANGE
[../]
# order parameter 2
[./eta2]
order = FIRST
family = LAGRANGE
[../]
[]
[Materials]
[./h_eta1]
type = SwitchingFunctionMaterial
h_order = SIMPLE
eta = eta1
function_name = h1
outputs = exodus
[../]
[./h_eta2]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta2
function_name = h2
outputs = exodus
[../]
[./g_eta1]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta1
function_name = g1
outputs = exodus
[../]
[./g_eta2]
type = BarrierFunctionMaterial
g_order = LOW
eta = eta2
function_name = g2
outputs = exodus
[../]
[]
[Kernels]
[./eta1diff]
type = Diffusion
variable = eta1
[../]
[./eta2diff]
type = Diffusion
variable = eta2
[../]
[]
[Executioner]
type = Steady
solve_type = 'PJFNK'
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(modules/phase_field/examples/kim-kim-suzuki/kks_example_dirichlet.i)
#
# KKS simple example in the split form
#
[Mesh]
type = GeneratedMesh
dim = 2
elem_type = QUAD4
nx = 50
ny = 2
nz = 0
xmin = 0
xmax = 20
ymin = 0
ymax = 0.4
zmin = 0
zmax = 0
[]
[AuxVariables]
[./Fglobal]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# hydrogen concentration
[./c]
order = FIRST
family = LAGRANGE
[../]
# chemical potential
[./w]
order = FIRST
family = LAGRANGE
[../]
# Liquid phase solute concentration
[./cl]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[../]
# Solid phase solute concentration
[./cs]
order = FIRST
family = LAGRANGE
initial_condition = 0.9
[../]
[]
[Functions]
[./ic_func_eta]
type = ParsedFunction
expression = 0.5*(1.0-tanh((x)/sqrt(2.0)))
[../]
[./ic_func_c]
type = ParsedFunction
expression = '0.9*(0.5*(1.0-tanh(x/sqrt(2.0))))^3*(6*(0.5*(1.0-tanh(x/sqrt(2.0))))^2-15*(0.5*(1.0-tanh(x/sqrt(2.0))))+10)+0.1*(1-(0.5*(1.0-tanh(x/sqrt(2.0))))^3*(6*(0.5*(1.0-tanh(x/sqrt(2.0))))^2-15*(0.5*(1.0-tanh(x/sqrt(2.0))))+10))'
[../]
[]
[ICs]
[./eta]
variable = eta
type = FunctionIC
function = ic_func_eta
[../]
[./c]
variable = c
type = FunctionIC
function = ic_func_c
[../]
[]
[BCs]
[./left_c]
type = DirichletBC
variable = 'c'
boundary = 'left'
value = 0.5
[../]
[./left_eta]
type = DirichletBC
variable = 'eta'
boundary = 'left'
value = 0.5
[../]
[]
[Materials]
# Free energy of the liquid
[./fl]
type = DerivativeParsedMaterial
property_name = fl
coupled_variables = 'cl'
expression = '(0.1-cl)^2'
[../]
# Free energy of the solid
[./fs]
type = DerivativeParsedMaterial
property_name = fs
coupled_variables = 'cs'
expression = '(0.9-cs)^2'
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L eps_sq'
prop_values = '0.7 0.7 1.0 '
[../]
[]
[Kernels]
# enforce c = (1-h(eta))*cl + h(eta)*cs
[./PhaseConc]
type = KKSPhaseConcentration
ca = cl
variable = cs
c = c
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotSolute]
type = KKSPhaseChemicalPotential
variable = cl
cb = cs
fa_name = fl
fb_name = fs
[../]
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = KKSSplitCHCRes
variable = c
ca = cl
fa_name = fl
w = w
[../]
[./dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[./ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = fl
fb_name = fs
w = 1.0
coupled_variables = 'cl cs'
[../]
[./ACBulkC]
type = KKSACBulkC
variable = eta
ca = cl
cb = cs
fa_name = fl
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = eps_sq
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[AuxKernels]
[./GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fl
fb_name = fs
w = 1.0
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'asm ilu nonzero'
l_max_its = 100
nl_max_its = 100
nl_abs_tol = 1e-10
end_time = 800
dt = 4.0
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Postprocessors]
[./dofs]
type = NumDOFs
[../]
[./integral]
type = ElementL2Error
variable = eta
function = ic_func_eta
[../]
[]
[Outputs]
exodus = true
console = true
gnuplot = true
[]
(modules/combined/test/tests/multiphase_mechanics/twophasestress.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 20
ny = 20
xmin = 0
xmax = 2
ymin = 0
ymax = 2
elem_type = QUAD4
[]
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[]
[AuxVariables]
[./eta]
[./InitialCondition]
type = FunctionIC
function = 'x/2'
[../]
[../]
[./e11_aux]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./matl_e11]
type = RankTwoAux
rank_two_tensor = stress
index_i = 0
index_j = 0
variable = e11_aux
[../]
[]
[Kernels]
[./TensorMechanics]
[../]
[]
[Materials]
[./elasticity_tensor_A]
type = ComputeElasticityTensor
base_name = A
fill_method = symmetric9
C_ijkl = '1e6 1e5 1e5 1e6 0 1e6 .4e6 .2e6 .5e6'
[../]
[./strain_A]
type = ComputeSmallStrain
base_name = A
eigenstrain_names = eigenstrain
[../]
[./stress_A]
type = ComputeLinearElasticStress
base_name = A
[../]
[./eigenstrain_A]
type = ComputeEigenstrain
base_name = A
eigen_base = '0.1 0.05 0 0 0 0.01'
prefactor = -1
eigenstrain_name = eigenstrain
[../]
[./elasticity_tensor_B]
type = ComputeElasticityTensor
base_name = B
fill_method = symmetric9
C_ijkl = '1e6 0 0 1e6 0 1e6 .5e6 .5e6 .5e6'
[../]
[./strain_B]
type = ComputeSmallStrain
base_name = B
eigenstrain_names = 'B_eigenstrain'
[../]
[./stress_B]
type = ComputeLinearElasticStress
base_name = B
[../]
[./eigenstrain_B]
type = ComputeEigenstrain
base_name = B
eigen_base = '0.1 0.05 0 0 0 0.01'
prefactor = -1
eigenstrain_name = 'B_eigenstrain'
[../]
[./switching]
type = SwitchingFunctionMaterial
eta = eta
[../]
[./combined]
type = TwoPhaseStressMaterial
base_A = A
base_B = B
[../]
[]
[BCs]
[./bottom_y]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0
[../]
[./left_x]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Steady
solve_type = 'NEWTON'
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(modules/phase_field/examples/slkks/CrFe.i)
#
# SLKKS two phase example for the BCC and SIGMA phases. The sigma phase contains
# multiple sublattices. Free energy from
# Jacob, Aurelie, Erwin Povoden-Karadeniz, and Ernst Kozeschnik. "Revised thermodynamic
# description of the Fe-Cr system based on an improved sublattice model of the sigma phase."
# Calphad 60 (2018): 16-28.
#
# In this simulation we consider diffusion (Cahn-Hilliard) and phase transformation.
#
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 160
ny = 1
nz = 0
xmin = -25
xmax = 25
ymin = -2.5
ymax = 2.5
elem_type = QUAD4
[]
[]
[AuxVariables]
[Fglobal]
order = CONSTANT
family = MONOMIAL
[]
[]
[Functions]
[sigma_cr0]
type = PiecewiseLinear
data_file = CrFe_sigma_out_var_0001.csv
format = columns
x_index_in_file = 5
y_index_in_file = 2
xy_in_file_only = false
[]
[sigma_cr1]
type = PiecewiseLinear
data_file = CrFe_sigma_out_var_0001.csv
format = columns
x_index_in_file = 5
y_index_in_file = 3
xy_in_file_only = false
[]
[sigma_cr2]
type = PiecewiseLinear
data_file = CrFe_sigma_out_var_0001.csv
format = columns
x_index_in_file = 5
y_index_in_file = 4
xy_in_file_only = false
[]
[]
[Variables]
# order parameters
[eta1]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
[eta2]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
# solute concentration
[cCr]
order = FIRST
family = LAGRANGE
[InitialCondition]
type = FunctionIC
function = '(x+25)/50*0.5+0.1'
[]
[]
# sublattice concentrations
[BCC_CR]
initial_condition = 0.45
[]
[SIGMA_0CR]
[InitialCondition]
type = CoupledValueFunctionIC
function = sigma_cr0
v = cCr
variable = SIGMA_0CR
[]
[]
[SIGMA_1CR]
[InitialCondition]
type = CoupledValueFunctionIC
function = sigma_cr1
v = cCr
variable = SIGMA_1CR
[]
[]
[SIGMA_2CR]
[InitialCondition]
type = CoupledValueFunctionIC
function = sigma_cr2
v = cCr
variable = SIGMA_2CR
[]
[]
# Lagrange multiplier
[lambda]
[]
[]
[Materials]
# CALPHAD free energies
[F_BCC_A2]
type = DerivativeParsedMaterial
property_name = F_BCC_A2
outputs = exodus
output_properties = F_BCC_A2
expression = 'BCC_FE:=1-BCC_CR; G := 8.3145*T*(1.0*if(BCC_CR > 1.0e-15,BCC_CR*log(BCC_CR),0) + '
'1.0*if(BCC_FE > 1.0e-15,BCC_FE*plog(BCC_FE,eps),0) + 3.0*if(BCC_VA > '
'1.0e-15,BCC_VA*log(BCC_VA),0))/(BCC_CR + BCC_FE) + 8.3145*T*if(T < '
'548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - 932.5*BCC_CR*BCC_FE*BCC_VA + '
'311.5*BCC_CR*BCC_VA - '
'1043.0*BCC_FE*BCC_VA,-8.13674105561218e-49*T^15/(0.525599232981783*BCC_CR*BCC_FE*BCC_'
'VA*(BCC_CR - BCC_FE) - 0.894055608820709*BCC_CR*BCC_FE*BCC_VA + '
'0.298657718120805*BCC_CR*BCC_VA - BCC_FE*BCC_VA + 9.58772770853308e-13)^15 - '
'4.65558036243985e-30*T^9/(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^9 - '
'1.3485349181899e-10*T^3/(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^3 + 1 - '
'0.905299382744392*(548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'932.5*BCC_CR*BCC_FE*BCC_VA + 311.5*BCC_CR*BCC_VA - 1043.0*BCC_FE*BCC_VA + '
'1.0e-9)/T,if(T < -548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'932.5*BCC_CR*BCC_FE*BCC_VA - 311.5*BCC_CR*BCC_VA + '
'1043.0*BCC_FE*BCC_VA,-8.13674105561218e-49*T^15/(-0.525599232981783*BCC_CR*BCC_FE*BCC'
'_VA*(BCC_CR - BCC_FE) + 0.894055608820709*BCC_CR*BCC_FE*BCC_VA - '
'0.298657718120805*BCC_CR*BCC_VA + BCC_FE*BCC_VA + 9.58772770853308e-13)^15 - '
'4.65558036243985e-30*T^9/(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) '
'+ 0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^9 - '
'1.3485349181899e-10*T^3/(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^3 + 1 - '
'0.905299382744392*(-548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'932.5*BCC_CR*BCC_FE*BCC_VA - 311.5*BCC_CR*BCC_VA + 1043.0*BCC_FE*BCC_VA + '
'1.0e-9)/T,if(T > -548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'932.5*BCC_CR*BCC_FE*BCC_VA - 311.5*BCC_CR*BCC_VA + 1043.0*BCC_FE*BCC_VA & '
'548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - 932.5*BCC_CR*BCC_FE*BCC_VA + '
'311.5*BCC_CR*BCC_VA - 1043.0*BCC_FE*BCC_VA < '
'0,-79209031311018.7*(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^5/T^5 - '
'3.83095660520737e+42*(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^15/T^15 - '
'1.22565886734485e+72*(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^25/T^25,if(T > '
'548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - 932.5*BCC_CR*BCC_FE*BCC_VA + '
'311.5*BCC_CR*BCC_VA - 1043.0*BCC_FE*BCC_VA & 548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - '
'BCC_FE) - 932.5*BCC_CR*BCC_FE*BCC_VA + 311.5*BCC_CR*BCC_VA - 1043.0*BCC_FE*BCC_VA > '
'0,-79209031311018.7*(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^5/T^5 - '
'3.83095660520737e+42*(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^15/T^15 - '
'1.22565886734485e+72*(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^25/T^25,0))))*log((2.15*BCC_CR*BCC_FE*BCC_VA - '
'0.008*BCC_CR*BCC_VA + 2.22*BCC_FE*BCC_VA)*if(2.15*BCC_CR*BCC_FE*BCC_VA - '
'0.008*BCC_CR*BCC_VA + 2.22*BCC_FE*BCC_VA <= 0,-1.0,1.0) + 1)/(BCC_CR + BCC_FE) + '
'1.0*(BCC_CR*BCC_VA*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + '
'BCC_FE*BCC_VA*if(T >= 298.15 & T < 1811.0,77358.5*1/T - 23.5143*T*log(T) + 124.134*T '
'- 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= 1811.0 & T < '
'6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581,0)))/(BCC_CR '
'+ BCC_FE) + 1.0*(BCC_CR*BCC_FE*BCC_VA*(500.0 - 1.5*T)*(BCC_CR - BCC_FE) + '
'BCC_CR*BCC_FE*BCC_VA*(24600.0 - 14.98*T) + BCC_CR*BCC_FE*BCC_VA*(9.15*T - '
'14000.0)*(BCC_CR - BCC_FE)^2)/(BCC_CR + BCC_FE); G/100000'
coupled_variables = 'BCC_CR'
constant_names = 'BCC_VA T eps'
constant_expressions = '1 1000 0.01'
[]
[F_SIGMA]
type = DerivativeParsedMaterial
property_name = F_SIGMA
outputs = exodus
output_properties = F_SIGMA
expression = 'SIGMA_0FE := 1-SIGMA_0CR; SIGMA_1FE := 1-SIGMA_1CR; SIGMA_2FE := 1-SIGMA_2CR; G := '
'8.3145*T*(10.0*if(SIGMA_0CR > 1.0e-15,SIGMA_0CR*plog(SIGMA_0CR,eps),0) + '
'10.0*if(SIGMA_0FE > 1.0e-15,SIGMA_0FE*plog(SIGMA_0FE,eps),0) + 4.0*if(SIGMA_1CR > '
'1.0e-15,SIGMA_1CR*plog(SIGMA_1CR,eps),0) + 4.0*if(SIGMA_1FE > '
'1.0e-15,SIGMA_1FE*plog(SIGMA_1FE,eps),0) + 16.0*if(SIGMA_2CR > '
'1.0e-15,SIGMA_2CR*plog(SIGMA_2CR,eps),0) + 16.0*if(SIGMA_2FE > '
'1.0e-15,SIGMA_2FE*plog(SIGMA_2FE,eps),0))/(10.0*SIGMA_0CR + 10.0*SIGMA_0FE + '
'4.0*SIGMA_1CR + 4.0*SIGMA_1FE + 16.0*SIGMA_2CR + 16.0*SIGMA_2FE) + '
'(SIGMA_0FE*SIGMA_1CR*SIGMA_2CR*SIGMA_2FE*(-70.0*T - 170400.0) + '
'SIGMA_0FE*SIGMA_1FE*SIGMA_2CR*SIGMA_2FE*(-10.0*T - 330839.0))/(10.0*SIGMA_0CR + '
'10.0*SIGMA_0FE + 4.0*SIGMA_1CR + 4.0*SIGMA_1FE + 16.0*SIGMA_2CR + 16.0*SIGMA_2FE) + '
'(SIGMA_0CR*SIGMA_1CR*SIGMA_2CR*(30.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - '
'26.908*T*log(T) + 157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= '
'2180.0 & T < 6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) '
'+ 132000.0) + SIGMA_0CR*SIGMA_1CR*SIGMA_2FE*(-110.0*T + 16.0*if(T >= 298.15 & T < '
'1811.0,77358.5*1/T - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - '
'5.89269e-8*T^3.0 + 1225.7,if(T >= 1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - '
'46.0*T*log(T) + 299.31255*T - 25383.581,0)) + 14.0*if(T >= 298.15 & T < '
'2180.0,139250.0*1/T - 26.908*T*log(T) + 157.48*T + 0.00189435*T^2.0 - '
'1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < 6000.0,-2.88526e+32*T^(-9.0) - '
'50.0*T*log(T) + 344.18*T - 34869.344,0)) + 123500.0) + '
'SIGMA_0CR*SIGMA_1FE*SIGMA_2CR*(4.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 26.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 140486.0) '
'+ SIGMA_0CR*SIGMA_1FE*SIGMA_2FE*(20.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 10.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 148800.0) '
'+ SIGMA_0FE*SIGMA_1CR*SIGMA_2CR*(10.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 20.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 56200.0) + '
'SIGMA_0FE*SIGMA_1CR*SIGMA_2FE*(26.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 4.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 152700.0) '
'+ SIGMA_0FE*SIGMA_1FE*SIGMA_2CR*(14.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 16.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 46200.0) + '
'SIGMA_0FE*SIGMA_1FE*SIGMA_2FE*(30.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 173333.0))/(10.0*SIGMA_0CR + 10.0*SIGMA_0FE + 4.0*SIGMA_1CR + '
'4.0*SIGMA_1FE + 16.0*SIGMA_2CR + 16.0*SIGMA_2FE); G/100000'
coupled_variables = 'SIGMA_0CR SIGMA_1CR SIGMA_2CR'
constant_names = 'T eps'
constant_expressions = '1000 0.01'
[]
# h(eta)
[h1]
type = SwitchingFunctionMaterial
function_name = h1
h_order = HIGH
eta = eta1
[]
[h2]
type = SwitchingFunctionMaterial
function_name = h2
h_order = HIGH
eta = eta2
[]
# g(eta)
[g1]
type = BarrierFunctionMaterial
function_name = g1
g_order = SIMPLE
eta = eta1
[]
[g2]
type = BarrierFunctionMaterial
function_name = g2
g_order = SIMPLE
eta = eta2
[]
# constant properties
[constants]
type = GenericConstantMaterial
prop_names = 'D L kappa'
prop_values = '10 1 0.1 '
[]
# Coefficients for diffusion equation
[Dh1]
type = DerivativeParsedMaterial
material_property_names = 'D h1(eta1)'
expression = D*h1
property_name = Dh1
coupled_variables = eta1
derivative_order = 1
[]
[Dh2a]
type = DerivativeParsedMaterial
material_property_names = 'D h2(eta2)'
expression = D*h2*10/30
property_name = Dh2a
coupled_variables = eta2
derivative_order = 1
[]
[Dh2b]
type = DerivativeParsedMaterial
material_property_names = 'D h2(eta2)'
expression = D*h2*4/30
property_name = Dh2b
coupled_variables = eta2
derivative_order = 1
[]
[Dh2c]
type = DerivativeParsedMaterial
material_property_names = 'D h2(eta2)'
expression = D*h2*16/30
property_name = Dh2c
coupled_variables = eta2
derivative_order = 1
[]
[]
[Kernels]
#Kernels for diffusion equation
[diff_time]
type = TimeDerivative
variable = cCr
[]
[diff_c1]
type = MatDiffusion
variable = cCr
diffusivity = Dh1
v = BCC_CR
coupled_variables = eta1
[]
[diff_c2a]
type = MatDiffusion
variable = cCr
diffusivity = Dh2a
v = SIGMA_0CR
coupled_variables = eta2
[]
[diff_c2b]
type = MatDiffusion
variable = cCr
diffusivity = Dh2b
v = SIGMA_1CR
coupled_variables = eta2
[]
[diff_c2c]
type = MatDiffusion
variable = cCr
diffusivity = Dh2c
v = SIGMA_2CR
coupled_variables = eta2
[]
# enforce pointwise equality of chemical potentials
[chempot1a2a]
# The BCC phase has only one sublattice
# we tie it to the first sublattice with site fraction 10/(10+4+16) in the sigma phase
type = KKSPhaseChemicalPotential
variable = BCC_CR
cb = SIGMA_0CR
kb = '${fparse 10/30}'
fa_name = F_BCC_A2
fb_name = F_SIGMA
args_b = 'SIGMA_1CR SIGMA_2CR'
[]
[chempot2a2b]
# This kernel ties the first two sublattices in the sigma phase together
type = SLKKSChemicalPotential
variable = SIGMA_0CR
a = 10
cs = SIGMA_1CR
as = 4
F = F_SIGMA
coupled_variables = 'SIGMA_2CR'
[]
[chempot2b2c]
# This kernel ties the remaining two sublattices in the sigma phase together
type = SLKKSChemicalPotential
variable = SIGMA_1CR
a = 4
cs = SIGMA_2CR
as = 16
F = F_SIGMA
coupled_variables = 'SIGMA_0CR'
[]
[phaseconcentration]
# This kernel ties the sum of the sublattice concentrations to the global concentration cCr
type = SLKKSMultiPhaseConcentration
variable = SIGMA_2CR
c = cCr
ns = '1 3'
as = '1 10 4 16'
cs = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR'
h_names = 'h1 h2'
eta = 'eta1 eta2'
[]
# Kernels for Allen-Cahn equation for eta1
[deta1dt]
type = TimeDerivative
variable = eta1
[]
[ACBulkF1]
type = KKSMultiACBulkF
variable = eta1
Fj_names = 'F_BCC_A2 F_SIGMA'
hj_names = 'h1 h2'
gi_name = g1
eta_i = eta1
wi = 0.1
coupled_variables = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR eta2'
[]
[ACBulkC1]
type = SLKKSMultiACBulkC
variable = eta1
F = F_BCC_A2
c = BCC_CR
ns = '1 3'
as = '1 10 4 16'
cs = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR'
h_names = 'h1 h2'
eta = 'eta1 eta2'
[]
[ACInterface1]
type = ACInterface
variable = eta1
kappa_name = kappa
[]
[lagrange1]
type = SwitchingFunctionConstraintEta
variable = eta1
h_name = h1
lambda = lambda
coupled_variables = 'eta2'
[]
# Kernels for Allen-Cahn equation for eta1
[deta2dt]
type = TimeDerivative
variable = eta2
[]
[ACBulkF2]
type = KKSMultiACBulkF
variable = eta2
Fj_names = 'F_BCC_A2 F_SIGMA'
hj_names = 'h1 h2'
gi_name = g2
eta_i = eta2
wi = 0.1
coupled_variables = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR eta1'
[]
[ACBulkC2]
type = SLKKSMultiACBulkC
variable = eta2
F = F_BCC_A2
c = BCC_CR
ns = '1 3'
as = '1 10 4 16'
cs = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR'
h_names = 'h1 h2'
eta = 'eta1 eta2'
[]
[ACInterface2]
type = ACInterface
variable = eta2
kappa_name = kappa
[]
[lagrange2]
type = SwitchingFunctionConstraintEta
variable = eta2
h_name = h2
lambda = lambda
coupled_variables = 'eta1'
[]
# Lagrange-multiplier constraint kernel for lambda
[lagrange]
type = SwitchingFunctionConstraintLagrange
variable = lambda
h_names = 'h1 h2'
etas = 'eta1 eta2'
epsilon = 1e-6
[]
[]
[AuxKernels]
[GlobalFreeEnergy]
type = KKSMultiFreeEnergy
variable = Fglobal
Fj_names = 'F_BCC_A2 F_SIGMA'
hj_names = 'h1 h2'
gj_names = 'g1 g2'
interfacial_vars = 'eta1 eta2'
kappa_names = 'kappa kappa'
w = 0.1
[]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
line_search = none
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -ksp_gmres_restart'
petsc_options_value = 'asm lu nonzero 30'
l_max_its = 100
nl_max_its = 20
nl_abs_tol = 1e-10
end_time = 10000
[TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 12
iteration_window = 2
growth_factor = 1.5
cutback_factor = 0.7
dt = 0.1
[]
[]
[VectorPostprocessors]
[var]
type = LineValueSampler
start_point = '-25 0 0'
end_point = '25 0 0'
variable = 'cCr eta1 eta2 SIGMA_0CR SIGMA_1CR SIGMA_2CR'
num_points = 151
sort_by = id
execute_on = 'initial timestep_end'
[]
[mat]
type = LineMaterialRealSampler
start = '-25 0 0'
end = '25 0 0'
property = 'F_BCC_A2 F_SIGMA'
sort_by = id
execute_on = 'initial timestep_end'
[]
[]
[Postprocessors]
[F]
type = ElementIntegralVariablePostprocessor
variable = Fglobal
execute_on = 'initial timestep_end'
[]
[cmin]
type = NodalExtremeValue
value_type = min
variable = cCr
execute_on = 'initial timestep_end'
[]
[cmax]
type = NodalExtremeValue
value_type = max
variable = cCr
execute_on = 'initial timestep_end'
[]
[ctotal]
type = ElementIntegralVariablePostprocessor
variable = cCr
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
exodus = true
print_linear_residuals = false
csv = true
perf_graph = true
[]
(modules/phase_field/test/tests/KKS_system/auxkernel.i)
#
# This test checks if the two phase and lagrange multiplier solutions can be replicated
# with a two order parameter approach, where the second order parameter eta2 is an
# auxiliary variable that is set as eta2 := 1 - eta1
# The solution is reproduced, but convergence is suboptimal, as important Jacobian
# terms for eta1 (that should come indirectly from eta2) are missing.
#
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 5
[]
[AuxVariables]
[Fglobal]
order = CONSTANT
family = MONOMIAL
[]
# order parameter 2
[eta2]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
[]
#
# With this approach the derivative w.r.t. eta1 is lost in all terms depending on
# eta2 a potential fix would be to make eta2 a material property with derivatives.
# This would require a major rewrite of the phase field kernels, though.
#
[AuxKernels]
[eta2]
type = ParsedAux
variable = eta2
expression = '1-eta1'
coupled_variables = eta1
[]
[]
[Variables]
# concentration
[c]
order = FIRST
family = LAGRANGE
[InitialCondition]
type = FunctionIC
function = x/5
[]
[]
# order parameter 1
[eta1]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
# phase concentration 1
[c1]
order = FIRST
family = LAGRANGE
initial_condition = 0.9
[]
# phase concentration 2
[c2]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[]
[]
[Materials]
# simple toy free energies
[f1] # = fd
type = DerivativeParsedMaterial
property_name = F1
coupled_variables = 'c1'
expression = '(0.9-c1)^2'
[]
[f2] # = fm
type = DerivativeParsedMaterial
property_name = F2
coupled_variables = 'c2'
expression = '(0.1-c2)^2'
[]
# Switching functions for each phase
[h1_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta1
function_name = h1
[]
[h2_eta]
type = DerivativeParsedMaterial
material_property_names = 'h1(eta1)'
expression = '1-h1'
property_name = h2
coupled_variables = eta1
[]
# Coefficients for diffusion equation
[Dh1]
type = DerivativeParsedMaterial
material_property_names = 'D h1(eta1)'
expression = D*h1
property_name = Dh1
coupled_variables = eta1
[]
[Dh2]
type = DerivativeParsedMaterial
material_property_names = 'D h2(eta1)'
expression = 'D*h2'
property_name = Dh2
coupled_variables = eta1
[]
# Barrier functions for each phase
[g1]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta1
function_name = g1
[]
[g2]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta2
function_name = g2
[]
# constant properties
[constants]
type = GenericConstantMaterial
prop_names = 'D L kappa'
prop_values = '0.7 0.7 0.2'
[]
[]
[Kernels]
#Kernels for diffusion equation
[diff_time]
type = TimeDerivative
variable = c
[]
[diff_c1]
type = MatDiffusion
variable = c
diffusivity = Dh1
v = c1
args = eta1
[]
[diff_c2]
type = MatDiffusion
variable = c
diffusivity = Dh2
v = c2
args = eta1
[]
# Kernels for Allen-Cahn equation for eta1
[deta1dt]
type = TimeDerivative
variable = eta1
[]
[ACBulkF1]
type = KKSMultiACBulkF
variable = eta1
Fj_names = 'F1 F2'
hj_names = 'h1 h2'
gi_name = g1
eta_i = eta1
wi = 0.2
coupled_variables = 'c1 c2 eta2'
[]
[ACBulkC1]
type = KKSMultiACBulkC
variable = eta1
Fj_names = 'F1 F2'
hj_names = 'h1 h2'
cj_names = 'c1 c2'
eta_i = eta1
coupled_variables = 'eta2'
[]
[ACInterface1]
type = ACInterface
variable = eta1
kappa_name = kappa
[]
# Phase concentration constraints
[chempot12]
type = KKSPhaseChemicalPotential
variable = c1
cb = c2
fa_name = F1
fb_name = F2
[]
[phaseconcentration]
type = KKSMultiPhaseConcentration
variable = c2
cj = 'c1 c2'
hj_names = 'h1 h2'
etas = 'eta1 eta2'
c = c
[]
[]
[AuxKernels]
[Fglobal_total]
type = KKSMultiFreeEnergy
Fj_names = 'F1 F2 '
hj_names = 'h1 h2 '
gj_names = 'g1 g2 '
variable = Fglobal
w = 0.2
interfacial_vars = 'eta1 eta2 '
kappa_names = 'kappa kappa'
[]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu '
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-10
nl_abs_tol = 1.0e-11
end_time = 350
dt = 10
[]
[Preconditioning]
[full]
type = SMP
full = true
[]
[]
[VectorPostprocessors]
[c]
type = LineValueSampler
variable = c
start_point = '0 0 0'
end_point = '5 0 0'
num_points = 21
sort_by = x
[]
[]
[Outputs]
csv = true
execute_on = FINAL
[]
(modules/combined/test/tests/multiphase_mechanics/multiphasestress.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 20
ny = 20
xmin = 0
xmax = 2
ymin = 0
ymax = 2
elem_type = QUAD4
[]
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Variables]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[]
[AuxVariables]
[./eta1]
[./InitialCondition]
type = FunctionIC
function = 'x/2'
[../]
[../]
[./eta2]
[./InitialCondition]
type = FunctionIC
function = 'y/2'
[../]
[../]
[./eta3]
[./InitialCondition]
type = FunctionIC
function = '(2^0.5-(y-1)^2=(y-1)^2)/2'
[../]
[../]
[./e11_aux]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./matl_e11]
type = RankTwoAux
rank_two_tensor = stress
index_i = 0
index_j = 0
variable = e11_aux
[../]
[]
[Kernels]
[./TensorMechanics]
[../]
[]
[Materials]
[./elasticity_tensor_A]
type = ComputeElasticityTensor
base_name = A
fill_method = symmetric9
C_ijkl = '1e6 1e5 1e5 1e6 0 1e6 .4e6 .2e6 .5e6'
[../]
[./strain_A]
type = ComputeSmallStrain
base_name = A
eigenstrain_names = eigenstrain
[../]
[./stress_A]
type = ComputeLinearElasticStress
base_name = A
[../]
[./eigenstrain_A]
type = ComputeEigenstrain
base_name = A
eigen_base = '0.1 0.05 0 0 0 0.01'
prefactor = -1
eigenstrain_name = eigenstrain
[../]
[./elasticity_tensor_B]
type = ComputeElasticityTensor
base_name = B
fill_method = symmetric9
C_ijkl = '1e6 0 0 1e6 0 1e6 .5e6 .5e6 .5e6'
[../]
[./strain_B]
type = ComputeSmallStrain
base_name = B
eigenstrain_names = 'B_eigenstrain'
[../]
[./stress_B]
type = ComputeLinearElasticStress
base_name = B
[../]
[./eigenstrain_B]
type = ComputeEigenstrain
base_name = B
eigen_base = '0.1 0.05 0 0 0 0.01'
prefactor = -1
eigenstrain_name = 'B_eigenstrain'
[../]
[./elasticity_tensor_C]
type = ComputeElasticityTensor
base_name = C
fill_method = symmetric9
C_ijkl = '1.1e6 1e5 0 1e6 0 1e6 .5e6 .2e6 .5e6'
[../]
[./strain_C]
type = ComputeSmallStrain
base_name = C
eigenstrain_names = 'C_eigenstrain'
[../]
[./stress_C]
type = ComputeLinearElasticStress
base_name = C
[../]
[./eigenstrain_C]
type = ComputeEigenstrain
base_name = C
eigen_base = '0.1 0.05 0 0 0 0.01'
prefactor = -1
eigenstrain_name = 'C_eigenstrain'
[../]
[./switching_A]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
[../]
[./switching_B]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
[../]
[./switching_C]
type = SwitchingFunctionMaterial
function_name = h3
eta = eta3
[../]
[./combined]
type = MultiPhaseStressMaterial
phase_base = 'A B C'
h = 'h1 h2 h3'
[../]
[]
[BCs]
[./bottom_y]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0
[../]
[./left_x]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Steady
solve_type = 'NEWTON'
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(modules/phase_field/test/tests/KKS_system/kks_example_nested.i)
#
# KKS toy problem in the split form
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 15
ny = 15
nz = 0
xmin = -2.5
xmax = 2.5
ymin = -2.5
ymax = 2.5
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[AuxVariables]
[./Fglobal]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# hydrogen concentration
[./c]
order = FIRST
family = LAGRANGE
[../]
# chemical potential
[./w]
order = FIRST
family = LAGRANGE
[../]
[]
[ICs]
[./eta]
variable = eta
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 1.5
invalue = 0.2
outvalue = 0.1
int_width = 0.75
[../]
[./c]
variable = c
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 1.5
invalue = 0.6
outvalue = 0.4
int_width = 0.75
[../]
[]
[BCs]
[./Periodic]
[./all]
variable = 'eta w c'
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
# Free energy of the matrix
[./fm]
type = DerivativeParsedMaterial
f_name = fm
function = '(0.1-cm)^2'
material_property_names = 'cm'
additional_derivative_symbols = 'cm'
compute = false
[../]
# Free energy of the delta phase
[./fd]
type = DerivativeParsedMaterial
f_name = fd
function = '(0.9-cd)^2'
material_property_names = 'cd'
additional_derivative_symbols = 'cd'
compute = false
[../]
# Compute phase concentrations
[./PhaseConcentrationMaterial]
type = KKSPhaseConcentrationMaterial
global_cs = 'c'
ci_names = 'cm cd'
ci_IC = '0 0'
fa_name = fm
fb_name = fd
h_name = h
min_iterations = 1
max_iterations = 100
absolute_tolerance = 1e-9
relative_tolerance = 1e-9
nested_iterations = iter
outputs = exodus
[../]
# Compute chain rule terms
[./PhaseConcentrationDerivatives]
type = KKSPhaseConcentrationDerivatives
global_cs = 'c'
eta = eta
ci_names = 'cm cd'
fa_name = fm
fb_name = fd
h_name = h
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L kappa'
prop_values = '0.7 0.7 0.4 '
[../]
[]
[Kernels]
# full transient
active = 'CHBulk ACBulkF ACBulkC ACInterface dcdt detadt ckernel'
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = NestedKKSSplitCHCRes
variable = c
global_cs = 'c'
w = w
all_etas = eta
ca_names = 'cm cd'
fa_name = fm
args = 'eta w'
[../]
[./dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[./ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = NestedKKSACBulkF
variable = eta
global_cs = 'c'
ci_names = 'cm cd'
fa_name = fm
fb_name = fd
g_name = g
h_name = h
mob_name = L
w = 0.4
args = 'c'
[../]
[./ACBulkC]
type = NestedKKSACBulkC
variable = eta
global_cs = 'c'
ci_names = 'cm cd'
fa_name = fm
h_name = h
mob_name = L
args = 'c'
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[AuxKernels]
[./GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fm
fb_name = fd
w = 0.4
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pctype -sub_pc_type -sub_pc_factor_shift_type -pc_factor_shift_type'
petsc_options_value = ' asm lu nonzero nonzero'
l_max_its = 100
nl_max_its = 100
num_steps = 3
dt = 0.1
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Outputs]
file_base = kks_example_nested
exodus = true
[]
(modules/combined/test/tests/surface_tension_KKS/surface_tension_KKS.i)
#
# KKS coupled with elasticity. Physical parameters for matrix and precipitate phases
# are gamma and gamma-prime phases, respectively, in the Ni-Al system.
# Parameterization is as described in L.K. Aagesen et al., Computational Materials
# Science, 140, 10-21 (2017), with isotropic elastic properties in both phases
# and without eigenstrain.
#
[Mesh]
type = GeneratedMesh
dim = 1
nx = 200
xmax = 200
[]
[Problem]
coord_type = RSPHERICAL
[]
[GlobalParams]
displacements = 'disp_x'
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# solute concentration
[./c]
order = FIRST
family = LAGRANGE
[../]
# chemical potential
[./w]
order = FIRST
family = LAGRANGE
[../]
# solute phase concentration (matrix)
[./cm]
order = FIRST
family = LAGRANGE
initial_condition = 0.13
[../]
# solute phase concentration (precipitate)
[./cp]
order = FIRST
family = LAGRANGE
initial_condition = 0.235
[../]
[]
[AuxVariables]
[./energy_density]
family = MONOMIAL
[../]
[./extra_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./extra_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./extra_zz]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_xx]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./strain_zz]
order = CONSTANT
family = MONOMIAL
[../]
[]
[ICs]
[./eta_ic]
variable = eta
type = FunctionIC
function = ic_func_eta
[../]
[./c_ic]
variable = c
type = FunctionIC
function = ic_func_c
[../]
[]
[Functions]
[./ic_func_eta]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2+z^2);0.5*(1.0-tanh((r-r0)/delta_eta/sqrt(2.0)))'
symbol_names = 'delta_eta r0'
symbol_values = '6.431 100'
[../]
[./ic_func_c]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2+z^2);eta_an:=0.5*(1.0-tanh((r-r0)/delta/sqrt(2.0)));0.235*eta_an^3*(6*eta_an^2-15*eta_an+10)+0.13*(1-eta_an^3*(6*eta_an^2-15*eta_an+10))'
symbol_names = 'delta r0'
symbol_values = '6.431 100'
[../]
[]
[Modules/TensorMechanics/Master]
[./all]
add_variables = true
generate_output = 'hydrostatic_stress stress_xx stress_yy stress_zz'
[../]
[]
[Kernels]
# enforce c = (1-h(eta))*cm + h(eta)*cp
[./PhaseConc]
type = KKSPhaseConcentration
ca = cm
variable = cp
c = c
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotVacancies]
type = KKSPhaseChemicalPotential
variable = cm
cb = cp
fa_name = f_total_matrix
fb_name = f_total_ppt
[../]
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = KKSSplitCHCRes
variable = c
ca = cm
fa_name = f_total_matrix
w = w
[../]
[./dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[./ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = f_total_matrix
fb_name = f_total_ppt
w = 0.0033
args = 'cp cm'
[../]
[./ACBulkC]
type = KKSACBulkC
variable = eta
ca = cm
cb = cp
fa_name = f_total_matrix
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[AuxKernels]
[./extra_xx]
type = RankTwoAux
rank_two_tensor = extra_stress
index_i = 0
index_j = 0
variable = extra_xx
[../]
[./extra_yy]
type = RankTwoAux
rank_two_tensor = extra_stress
index_i = 1
index_j = 1
variable = extra_yy
[../]
[./extra_zz]
type = RankTwoAux
rank_two_tensor = extra_stress
index_i = 2
index_j = 2
variable = extra_zz
[../]
[./strain_xx]
type = RankTwoAux
rank_two_tensor = mechanical_strain
index_i = 0
index_j = 0
variable = strain_xx
[../]
[./strain_yy]
type = RankTwoAux
rank_two_tensor = mechanical_strain
index_i = 1
index_j = 1
variable = strain_yy
[../]
[./strain_zz]
type = RankTwoAux
rank_two_tensor = mechanical_strain
index_i = 2
index_j = 2
variable = strain_zz
[../]
[]
[Materials]
# Chemical free energy of the matrix
[./fm]
type = DerivativeParsedMaterial
property_name = fm
coupled_variables = 'cm'
expression = '6.55*(cm-0.13)^2'
[../]
# Elastic energy of the matrix
[./elastic_free_energy_m]
type = ElasticEnergyMaterial
base_name = matrix
f_name = fe_m
args = ' '
[../]
# Total free energy of the matrix
[./Total_energy_matrix]
type = DerivativeSumMaterial
property_name = f_total_matrix
sum_materials = 'fm fe_m'
coupled_variables = 'cm'
[../]
# Free energy of the precipitate phase
[./fp]
type = DerivativeParsedMaterial
property_name = fp
coupled_variables = 'cp'
expression = '6.55*(cp-0.235)^2'
[../]
# Elastic energy of the precipitate
[./elastic_free_energy_p]
type = ElasticEnergyMaterial
base_name = ppt
f_name = fe_p
args = ' '
[../]
# Total free energy of the precipitate
[./Total_energy_ppt]
type = DerivativeSumMaterial
property_name = f_total_ppt
sum_materials = 'fp fe_p'
coupled_variables = 'cp'
[../]
# Total elastic energy
[./Total_elastic_energy]
type = DerivativeTwoPhaseMaterial
eta = eta
f_name = f_el_mat
fa_name = fe_m
fb_name = fe_p
outputs = exodus
W = 0
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
outputs = exodus
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L kappa'
prop_values = '0.7 0.7 0.1365'
[../]
#Mechanical properties
[./Stiffness_matrix]
type = ComputeElasticityTensor
C_ijkl = '74.25 14.525'
base_name = matrix
fill_method = symmetric_isotropic
[../]
[./Stiffness_ppt]
type = ComputeElasticityTensor
C_ijkl = '74.25 14.525'
base_name = ppt
fill_method = symmetric_isotropic
[../]
[./strain_matrix]
type = ComputeRSphericalSmallStrain
base_name = matrix
[../]
[./strain_ppt]
type = ComputeRSphericalSmallStrain
base_name = ppt
[../]
[./stress_matrix]
type = ComputeLinearElasticStress
base_name = matrix
[../]
[./stress_ppt]
type = ComputeLinearElasticStress
base_name = ppt
[../]
[./global_stress]
type = TwoPhaseStressMaterial
base_A = matrix
base_B = ppt
[../]
[./interface_stress]
type = ComputeSurfaceTensionKKS
v = eta
kappa_name = kappa
w = 0.0033
[../]
[]
[BCs]
[./left_r]
type = DirichletBC
variable = disp_x
boundary = left
value = 0
[../]
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'asm lu nonzero'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-9
nl_abs_tol = 1.0e-10
num_steps = 2
dt = 0.5
[]
[Outputs]
exodus = true
[./csv]
type = CSV
execute_on = 'final'
[../]
[]
(modules/phase_field/test/tests/KKS_system/two_phase.i)
#
# This test ensures that the equilibrium solution using the dedicated two phase
# formulation is identical to the two order parameters with a Lagrange multiplier
# constraint in lagrange_multiplier.i
#
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 5
[]
[AuxVariables]
[Fglobal]
order = CONSTANT
family = MONOMIAL
[]
[]
[Variables]
# order parameter
[eta]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
# hydrogen concentration
[c]
order = FIRST
family = LAGRANGE
[InitialCondition]
type = FunctionIC
function = x/5
[]
[]
# chemical potential
[w]
order = FIRST
family = LAGRANGE
[]
# hydrogen phase concentration (matrix)
[cm]
order = FIRST
family = LAGRANGE
initial_condition = 0.2
[]
# hydrogen phase concentration (delta phase)
[cd]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
[]
[Materials]
# Free energy of the matrix
[fm]
type = DerivativeParsedMaterial
property_name = fm
coupled_variables = 'cm'
expression = '(0.1-cm)^2'
[]
# Free energy of the delta phase
[fd]
type = DerivativeParsedMaterial
property_name = fd
coupled_variables = 'cd'
expression = '(0.9-cd)^2'
[]
# h(eta)
[h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[]
# g(eta)
[g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[]
# constant properties
[constants]
type = GenericConstantMaterial
prop_names = 'M L kappa'
prop_values = '0.7 0.7 0.4 '
[]
[]
[Kernels]
# full transient
active = 'PhaseConc ChemPotVacancies CHBulk ACBulkF ACBulkC ACInterface dcdt detadt ckernel'
# enforce c = (1-h(eta))*cm + h(eta)*cd
[PhaseConc]
type = KKSPhaseConcentration
ca = cm
variable = cd
c = c
eta = eta
[]
# enforce pointwise equality of chemical potentials
[ChemPotVacancies]
type = KKSPhaseChemicalPotential
variable = cm
cb = cd
fa_name = fm
fb_name = fd
[]
#
# Cahn-Hilliard Equation
#
[CHBulk]
type = KKSSplitCHCRes
variable = c
ca = cm
fa_name = fm
w = w
[]
[dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[]
[ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[]
#
# Allen-Cahn Equation
#
[ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = fm
fb_name = fd
coupled_variables = 'cm cd'
w = 0.4
[]
[ACBulkC]
type = KKSACBulkC
variable = eta
ca = cm
cb = cd
fa_name = fm
mob_name = L
[]
[ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa
mob_name = L
[]
[detadt]
type = TimeDerivative
variable = eta
[]
[]
[AuxKernels]
[GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fm
fb_name = fd
w = 0.4
[]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pctype -sub_pc_type -sub_pc_factor_shift_type -pc_factor_shift_type'
petsc_options_value = ' asm lu nonzero nonzero'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-10
nl_abs_tol = 1.0e-11
num_steps = 35
dt = 10
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[full]
type = SMP
full = true
[]
[]
[VectorPostprocessors]
[c]
type = LineValueSampler
variable = c
start_point = '0 0 0'
end_point = '5 0 0'
num_points = 21
sort_by = x
[]
[]
[Outputs]
csv = true
execute_on = FINAL
[]
(modules/phase_field/test/tests/MultiPhase/derivativetwophasematerial.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 14
ny = 10
nz = 0
xmin = 10
xmax = 40
ymin = 15
ymax = 35
elem_type = QUAD4
[]
[Variables]
[./c]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 25.0
y1 = 25.0
radius = 6.0
invalue = 0.9
outvalue = 0.1
int_width = 3.0
[../]
[../]
[./w]
order = FIRST
family = LAGRANGE
[../]
[./eta]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 30.0
y1 = 25.0
radius = 4.0
invalue = 0.9
outvalue = 0.1
int_width = 2.0
[../]
[../]
[]
[Kernels]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[./ACBulk]
type = AllenCahn
variable = eta
coupled_variables = c
f_name = F
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa_eta
[../]
[./c_res]
type = SplitCHParsed
variable = c
f_name = F
kappa_name = kappa_c
w = w
coupled_variables = 'eta'
[../]
[./w_res]
type = SplitCHWRes
variable = w
mob_name = M
[../]
[./time]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[]
[BCs]
[./Periodic]
[./All]
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
[./consts]
type = GenericConstantMaterial
prop_names = 'L kappa_eta'
prop_values = '1 1 '
[../]
[./consts2]
type = GenericConstantMaterial
prop_names = 'M kappa_c'
prop_values = '1 1'
[../]
[./switching]
type = SwitchingFunctionMaterial
eta = eta
h_order = SIMPLE
[../]
[./barrier]
type = BarrierFunctionMaterial
eta = eta
g_order = SIMPLE
[../]
[./free_energy_A]
type = DerivativeParsedMaterial
property_name = Fa
coupled_variables = 'c'
expression = '(c-0.1)^2*(c-1)^2 + c*0.01'
derivative_order = 2
enable_jit = true
[../]
[./free_energy_B]
type = DerivativeParsedMaterial
property_name = Fb
coupled_variables = 'c'
expression = 'c^2*(c-0.9)^2 + (1-c)*0.01'
derivative_order = 2
enable_jit = true
[../]
[./free_energy]
type = DerivativeTwoPhaseMaterial
property_name = F
fa_name = Fa
fb_name = Fb
coupled_variables = 'c'
eta = eta
derivative_order = 2
outputs = exodus
output_properties = 'F dF/dc dF/deta d^2F/dc^2 d^2F/dcdeta d^2F/deta^2'
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = 'bdf2'
solve_type = 'NEWTON'
l_max_its = 15
l_tol = 1.0e-4
nl_max_its = 10
nl_rel_tol = 1.0e-11
start_time = 0.0
num_steps = 1
dt = 0.1
[]
[Outputs]
exodus = true
[]
(modules/phase_field/test/tests/MultiPhase/acmultiinterface_aux.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 20
ny = 10
nz = 0
xmin = -10
xmax = 10
ymin = -5
ymax = 5
elem_type = QUAD4
[]
[AuxVariables]
[./eta1]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = -3.5
y1 = 0.0
radius = 4.0
invalue = 0.9
outvalue = 0.1
int_width = 2.0
[../]
[../]
[]
[Variables]
[./eta2]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 3.5
y1 = 0.0
radius = 4.0
invalue = 0.9
outvalue = 0.1
int_width = 2.0
[../]
[../]
[./eta3]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SpecifiedSmoothCircleIC
x_positions = '-4.0 4.0'
y_positions = ' 0.0 0.0'
z_positions = ' 0.0 0.0'
radii = '4.0 4.0'
invalue = 0.1
outvalue = 0.9
int_width = 2.0
[../]
[../]
[./lambda]
order = FIRST
family = LAGRANGE
initial_condition = 1.0
[../]
[]
[Kernels]
[./deta2dt]
type = TimeDerivative
variable = eta2
[../]
[./ACBulk2]
type = AllenCahn
variable = eta2
coupled_variables = 'eta1 eta3'
mob_name = L2
f_name = F
[../]
[./ACInterface2]
type = ACMultiInterface
variable = eta2
etas = 'eta1 eta2 eta3'
mob_name = L2
kappa_names = 'kappa21 kappa22 kappa23'
[../]
[./lagrange2]
type = SwitchingFunctionConstraintEta
variable = eta2
h_name = h2
lambda = lambda
[../]
[./deta3dt]
type = TimeDerivative
variable = eta3
[../]
[./ACBulk3]
type = AllenCahn
variable = eta3
coupled_variables = 'eta1 eta2'
mob_name = L3
f_name = F
[../]
[./ACInterface3]
type = ACMultiInterface
variable = eta3
etas = 'eta1 eta2 eta3'
mob_name = L3
kappa_names = 'kappa31 kappa32 kappa33'
[../]
[./lagrange3]
type = SwitchingFunctionConstraintEta
variable = eta3
h_name = h3
lambda = lambda
[../]
[./lagrange]
type = SwitchingFunctionConstraintLagrange
variable = lambda
etas = 'eta1 eta2 eta3'
h_names = 'h1 h2 h3'
epsilon = 0
[../]
[]
[BCs]
[./Periodic]
[./All]
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
[./consts]
type = GenericConstantMaterial
prop_names = 'Fx L1 L2 L3 kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
prop_values = '0 1 1 1 1 1 1 1 1 1 1 1 1 '
[../]
[./switching1]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
h_order = SIMPLE
[../]
[./switching2]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
h_order = SIMPLE
[../]
[./switching3]
type = SwitchingFunctionMaterial
function_name = h3
eta = eta3
h_order = SIMPLE
[../]
[./barrier]
type = MultiBarrierFunctionMaterial
etas = 'eta1 eta2 eta3'
[../]
[./free_energy]
type = DerivativeMultiPhaseMaterial
property_name = F
# we use a constant free energy (GeneriConstantmaterial property Fx)
fi_names = 'Fx Fx Fx'
hi_names = 'h1 h2 h3'
etas = 'eta1 eta2 eta3'
# the free energy is given by the MultiBarrierFunctionMaterial only
W = 1
derivative_order = 2
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = 'bdf2'
solve_type = 'PJFNK'
#petsc_options = '-snes_ksp -snes_ksp_ew'
#petsc_options = '-ksp_monitor_snes_lg-snes_ksp_ew'
#petsc_options_iname = '-ksp_gmres_restart'
#petsc_options_value = '1000 '
l_max_its = 15
l_tol = 1.0e-6
nl_max_its = 50
nl_rel_tol = 1.0e-8
nl_abs_tol = 1.0e-10
start_time = 0.0
num_steps = 2
dt = 0.2
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(modules/phase_field/test/tests/MultiPhase/lagrangemult.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 14
ny = 10
nz = 0
xmin = 10
xmax = 40
ymin = 15
ymax = 35
elem_type = QUAD4
[]
[Variables]
[./c]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 25.0
y1 = 25.0
radius = 6.0
invalue = 0.9
outvalue = 0.1
int_width = 3.0
[../]
[../]
[./w]
order = FIRST
family = LAGRANGE
[../]
[./eta1]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 30.0
y1 = 25.0
radius = 4.0
invalue = 0.9
outvalue = 0.1
int_width = 2.0
[../]
[../]
[./eta2]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[../]
[./lambda]
order = FIRST
family = LAGRANGE
initial_condition = 1.0
[../]
[]
[Kernels]
[./deta1dt]
type = TimeDerivative
variable = eta1
[../]
[./ACBulk1]
type = AllenCahn
variable = eta1
coupled_variables = 'c eta2'
f_name = F
[../]
[./ACInterface1]
type = ACInterface
variable = eta1
kappa_name = kappa_eta
[../]
[./lagrange1]
type = SwitchingFunctionConstraintEta
variable = eta1
h_name = h1
lambda = lambda
[../]
[./deta2dt]
type = TimeDerivative
variable = eta2
[../]
[./ACBulk2]
type = AllenCahn
variable = eta2
coupled_variables = 'c eta1'
f_name = F
[../]
[./ACInterface2]
type = ACInterface
variable = eta2
kappa_name = kappa_eta
[../]
[./lagrange2]
type = SwitchingFunctionConstraintEta
variable = eta2
h_name = h2
lambda = lambda
[../]
[./lagrange]
type = SwitchingFunctionConstraintLagrange
variable = lambda
etas = 'eta1 eta2'
h_names = 'h1 h2'
epsilon = 0
[../]
[./c_res]
type = SplitCHParsed
variable = c
f_name = F
kappa_name = kappa_c
w = w
coupled_variables = 'eta1 eta2'
[../]
[./w_res]
type = SplitCHWRes
variable = w
mob_name = M
[../]
[./time1]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[]
[BCs]
[./Periodic]
[./All]
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
[./consts]
type = GenericConstantMaterial
prop_names = 'L kappa_eta'
prop_values = '1 1 '
[../]
[./consts2]
type = GenericConstantMaterial
prop_names = 'M kappa_c'
prop_values = '1 1'
[../]
[./switching1]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
h_order = SIMPLE
outputs = exodus
[../]
[./switching2]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
h_order = SIMPLE
outputs = exodus
[../]
[./barrier]
type = MultiBarrierFunctionMaterial
etas = 'eta1 eta2'
[../]
[./free_energy_A]
type = DerivativeParsedMaterial
property_name = Fa
coupled_variables = 'c'
expression = '(c-0.1)^2'
derivative_order = 2
enable_jit = true
[../]
[./free_energy_B]
type = DerivativeParsedMaterial
property_name = Fb
coupled_variables = 'c'
expression = '(c-0.9)^2'
derivative_order = 2
enable_jit = true
[../]
[./free_energy]
type = DerivativeMultiPhaseMaterial
property_name = F
fi_names = 'Fa Fb'
hi_names = 'h1 h2'
etas = 'eta1 eta2'
coupled_variables = 'c'
derivative_order = 2
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = 'bdf2'
solve_type = 'PJFNK'
#petsc_options = '-snes_ksp -snes_ksp_ew'
#petsc_options = '-ksp_monitor_snes_lg-snes_ksp_ew'
#petsc_options_iname = '-ksp_gmres_restart'
#petsc_options_value = '1000 '
l_max_its = 15
l_tol = 1.0e-6
nl_max_its = 50
nl_rel_tol = 1.0e-8
nl_abs_tol = 1.0e-10
start_time = 0.0
num_steps = 1
dt = 0.01
dtmin = 0.01
[]
[Debug]
# show_var_residual_norms = true
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(modules/phase_field/test/tests/KKS_system/kks_phase_concentration.i)
#
# This test validates the phase concentration calculation for the KKS system
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 20
ny = 20
nz = 0
xmin = 0
xmax = 1
ymin = 0
ymax = 1
zmin = 0
zmax = 0
elem_type = QUAD4
[]
# We set c and eta...
[BCs]
# (and ca for debugging purposes)
[./left]
type = DirichletBC
variable = c
boundary = 'left'
value = 0.1
[../]
[./right]
type = DirichletBC
variable = c
boundary = 'right'
value = 0.9
[../]
[./top]
type = DirichletBC
variable = eta
boundary = 'top'
value = 0.1
[../]
[./bottom]
type = DirichletBC
variable = eta
boundary = 'bottom'
value = 0.9
[../]
[]
[Variables]
# concentration
[./c]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[../]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[../]
# phase concentration a
[./ca]
order = FIRST
family = LAGRANGE
initial_condition = 0.2
[../]
# phase concentration b
[./cb]
order = FIRST
family = LAGRANGE
initial_condition = 0.3
[../]
[]
[Materials]
# simple toy free energy
[./fa]
type = DerivativeParsedMaterial
property_name = Fa
coupled_variables = 'ca'
expression = 'ca^2'
[../]
[./fb]
type = DerivativeParsedMaterial
property_name = Fb
coupled_variables = 'cb'
expression = '(1-cb)^2'
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
outputs = exodus
[../]
[]
[Kernels]
active = 'cdiff etadiff phaseconcentration chempot'
##active = 'cbdiff cdiff etadiff chempot'
#active = 'cadiff cdiff etadiff phaseconcentration'
##active = 'cadiff cbdiff cdiff etadiff'
[./cadiff]
type = Diffusion
variable = ca
[../]
[./cbdiff]
type = Diffusion
variable = cb
[../]
[./cdiff]
type = Diffusion
variable = c
[../]
[./etadiff]
type = Diffusion
variable = eta
[../]
# ...and solve for ca and cb
[./phaseconcentration]
type = KKSPhaseConcentration
ca = ca
variable = cb
c = c
eta = eta
[../]
[./chempot]
type = KKSPhaseChemicalPotential
variable = ca
cb = cb
fa_name = Fa
fb_namee = Fb
[../]
[]
[Executioner]
type = Steady
solve_type = 'PJFNK'
#solve_type = 'NEWTON'
petsc_options_iname = '-pctype -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = ' asm lu nonzero'
[]
[Preconditioning]
active = 'full'
#active = 'mydebug'
#active = ''
[./full]
type = SMP
full = true
[../]
[./mydebug]
type = FDP
full = true
[../]
[]
[Outputs]
execute_on = 'timestep_end'
file_base = kks_phase_concentration
exodus = true
[]
(modules/combined/examples/publications/rapid_dev/fig7b.i)
#
# Fig. 7 input for 10.1016/j.commatsci.2017.02.017
# D. Schwen et al./Computational Materials Science 132 (2017) 36-45
# Dashed black curve (2)
# Eigenstrain is globally applied. Single global elastic free energies.
# Supply the RADIUS parameter (10-35) on the command line to generate data
# for all curves in the plot.
#
[Mesh]
type = GeneratedMesh
dim = 1
nx = 32
xmin = 0
xmax = 100
second_order = true
[]
[Problem]
coord_type = RSPHERICAL
[]
[GlobalParams]
displacements = 'disp_r'
[]
[Functions]
[./diff]
type = ParsedFunction
expression = '${RADIUS}-pos_c'
symbol_names = pos_c
symbol_values = pos_c
[../]
[]
# AuxVars to compute the free energy density for outputting
[AuxVariables]
[./local_energy]
order = CONSTANT
family = MONOMIAL
[../]
[./cross_energy]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./local_free_energy]
type = TotalFreeEnergy
variable = local_energy
interfacial_vars = 'c'
kappa_names = 'kappa_c'
execute_on = 'INITIAL TIMESTEP_END'
[../]
[]
[Variables]
# Solute concentration variable
[./c]
[./InitialCondition]
type = SmoothCircleIC
invalue = 1
outvalue = 0
x1 = 0
y1 = 0
radius = ${RADIUS}
int_width = 3
[../]
[../]
[./w]
[../]
# Phase order parameter
[./eta]
[./InitialCondition]
type = SmoothCircleIC
invalue = 1
outvalue = 0
x1 = 0
y1 = 0
radius = ${RADIUS}
int_width = 3
[../]
[../]
[./Fe_fit]
order = SECOND
[../]
[]
[Modules/TensorMechanics/Master/all]
add_variables = true
eigenstrain_names = eigenstrain
[]
[Kernels]
# Split Cahn-Hilliard kernels
[./c_res]
type = SplitCHParsed
variable = c
f_name = F
args = 'eta'
kappa_name = kappa_c
w = w
[../]
[./wres]
type = SplitCHWRes
variable = w
mob_name = M
[../]
[./time]
type = CoupledTimeDerivative
variable = w
v = c
[../]
# Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 1
[./detadt]
type = TimeDerivative
variable = eta
[../]
[./ACBulk1]
type = AllenCahn
variable = eta
args = 'c'
mob_name = L
f_name = F
[../]
[./ACInterface]
type = ACInterface
variable = eta
mob_name = L
kappa_name = kappa_eta
[../]
[./Fe]
type = MaterialPropertyValue
prop_name = Fe
variable = Fe_fit
[../]
[./autoadjust]
type = MaskedBodyForce
variable = w
function = diff
mask = mask
[../]
[]
[Materials]
# declare a few constants, such as mobilities (L,M) and interface gradient prefactors (kappa*)
[./consts]
type = GenericConstantMaterial
prop_names = 'M L kappa_c kappa_eta'
prop_values = '1.0 1.0 0.5 1'
[../]
# forcing function mask
[./mask]
type = ParsedMaterial
property_name = mask
expression = grad/dt
material_property_names = 'grad dt'
[../]
[./grad]
type = VariableGradientMaterial
variable = c
prop = grad
[../]
[./time]
type = TimeStepMaterial
[../]
# global mechanical properties
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1'
fill_method = symmetric_isotropic
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
# eigenstrain as a function of phase
[./eigenstrain]
type = ComputeVariableEigenstrain
eigen_base = '0.05 0.05 0.05 0 0 0'
prefactor = h
args = eta
eigenstrain_name = eigenstrain
[../]
# switching functions
[./switching]
type = SwitchingFunctionMaterial
function_name = h
eta = eta
h_order = SIMPLE
[../]
[./barrier]
type = BarrierFunctionMaterial
eta = eta
[../]
# chemical free energies
[./chemical_free_energy_1]
type = DerivativeParsedMaterial
property_name = Fc1
expression = 'c^2'
coupled_variables = 'c'
derivative_order = 2
[../]
[./chemical_free_energy_2]
type = DerivativeParsedMaterial
property_name = Fc2
expression = '(1-c)^2'
coupled_variables = 'c'
derivative_order = 2
[../]
# global chemical free energy
[./chemical_free_energy]
type = DerivativeTwoPhaseMaterial
f_name = Fc
fa_name = Fc1
fb_name = Fc2
eta = eta
args = 'c'
W = 4
[../]
# global elastic free energy
[./elastic_free_energy]
type = ElasticEnergyMaterial
f_name = Fe
args = 'eta'
output_properties = Fe
derivative_order = 2
[../]
# free energy
[./free_energy]
type = DerivativeSumMaterial
property_name = F
sum_materials = 'Fc Fe'
coupled_variables = 'c eta'
derivative_order = 2
[../]
[]
[BCs]
[./left_r]
type = DirichletBC
variable = disp_r
boundary = 'left'
value = 0
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
# We monitor the total free energy and the total solute concentration (should be constant)
[Postprocessors]
[./total_free_energy]
type = ElementIntegralVariablePostprocessor
variable = local_energy
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[./total_solute]
type = ElementIntegralVariablePostprocessor
variable = c
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[./pos_c]
type = FindValueOnLine
start_point = '0 0 0'
end_point = '100 0 0'
v = c
target = 0.582
tol = 1e-8
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[./pos_eta]
type = FindValueOnLine
start_point = '0 0 0'
end_point = '100 0 0'
v = eta
target = 0.5
tol = 1e-8
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[./c_min]
type = ElementExtremeValue
value_type = min
variable = c
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[]
[VectorPostprocessors]
[./line]
type = LineValueSampler
variable = 'Fe_fit c w'
start_point = '0 0 0'
end_point = '100 0 0'
num_points = 5000
sort_by = x
outputs = vpp
execute_on = 'INITIAL TIMESTEP_END'
[../]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type'
petsc_options_value = 'asm lu'
l_max_its = 30
nl_max_its = 15
l_tol = 1.0e-4
nl_rel_tol = 1.0e-8
nl_abs_tol = 2.0e-9
start_time = 0.0
end_time = 100000.0
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 8
iteration_window = 1
dt = 1
[../]
[./Adaptivity]
initial_adaptivity = 5
interval = 10
max_h_level = 5
refine_fraction = 0.9
coarsen_fraction = 0.1
[../]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
execute_on = 'INITIAL TIMESTEP_END'
[./table]
type = CSV
delimiter = ' '
file_base = radius_${RADIUS}/eigenstrain_pp
[../]
[./vpp]
type = CSV
delimiter = ' '
sync_times = '10 50 100 500 1000 5000 10000 50000 100000'
sync_only = true
time_data = true
file_base = radius_${RADIUS}/eigenstrain_vpp
[../]
[]
(modules/combined/test/tests/surface_tension_KKS/surface_tension_VDWgas.i)
# Test for ComputeExtraStressVDWGas
# Gas bubble with r = 15 nm in a solid matrix
# The gas pressure is counterbalanced by the surface tension of the solid-gas interface,
# which is included with ComputeSurfaceTensionKKS
[Mesh]
type = GeneratedMesh
dim = 1
nx = 300
xmin = 0
xmax = 30
[]
[Problem]
coord_type = RSPHERICAL
[]
[GlobalParams]
displacements = 'disp_x'
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# gas concentration
[./cg]
order = FIRST
family = LAGRANGE
[../]
# vacancy concentration
[./cv]
order = FIRST
family = LAGRANGE
[../]
# gas chemical potential
[./wg]
order = FIRST
family = LAGRANGE
[../]
# vacancy chemical potential
[./wv]
order = FIRST
family = LAGRANGE
[../]
# Matrix phase gas concentration
[./cgm]
order = FIRST
family = LAGRANGE
initial_condition = 1.01e-31
[../]
# Matrix phase vacancy concentration
[./cvm]
order = FIRST
family = LAGRANGE
initial_condition = 2.25e-11
[../]
# Bubble phase gas concentration
[./cgb]
order = FIRST
family = LAGRANGE
initial_condition = 0.2714
[../]
# Bubble phase vacancy concentration
[./cvb]
order = FIRST
family = LAGRANGE
initial_condition = 0.7286
[../]
[]
[ICs]
[./eta_ic]
variable = eta
type = FunctionIC
function = ic_func_eta
[../]
[./cv_ic]
variable = cv
type = FunctionIC
function = ic_func_cv
[../]
[./cg_ic]
variable = cg
type = FunctionIC
function = ic_func_cg
[../]
[]
[Functions]
[./ic_func_eta]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2+z^2);0.5*(1.0-tanh((r-r0)/delta_eta/sqrt(2.0)))'
symbol_names = 'delta_eta r0'
symbol_values = '0.321 15'
[../]
[./ic_func_cv]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2+z^2);eta_an:=0.5*(1.0-tanh((r-r0)/delta/sqrt(2.0)));cvbubinit*eta_an^3*(6*eta_an^2-15*eta_an+10)+cvmatrixinit*(1-eta_an^3*(6*eta_an^2-15*eta_an+10))'
symbol_names = 'delta r0 cvbubinit cvmatrixinit'
symbol_values = '0.321 15 0.7286 2.25e-11'
[../]
[./ic_func_cg]
type = ParsedFunction
expression = 'r:=sqrt(x^2+y^2+z^2);eta_an:=0.5*(1.0-tanh((r-r0)/delta/sqrt(2.0)));cgbubinit*eta_an^3*(6*eta_an^2-15*eta_an+10)+cgmatrixinit*(1-eta_an^3*(6*eta_an^2-15*eta_an+10))'
symbol_names = 'delta r0 cgbubinit cgmatrixinit'
symbol_values = '0.321 15 0.2714 1.01e-31'
[../]
[]
[Modules/TensorMechanics/Master]
[./all]
add_variables = true
generate_output = 'hydrostatic_stress stress_xx stress_yy stress_zz'
[../]
[]
[Kernels]
# enforce cg = (1-h(eta))*cgm + h(eta)*cgb
[./PhaseConc_g]
type = KKSPhaseConcentration
ca = cgm
variable = cgb
c = cg
eta = eta
[../]
# enforce cv = (1-h(eta))*cvm + h(eta)*cvb
[./PhaseConc_v]
type = KKSPhaseConcentration
ca = cvm
variable = cvb
c = cv
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotVacancies]
type = KKSPhaseChemicalPotential
variable = cvm
cb = cvb
fa_name = f_total_matrix
fb_name = f_total_bub
args_a = 'cgm'
args_b = 'cgb'
[../]
[./ChemPotGas]
type = KKSPhaseChemicalPotential
variable = cgm
cb = cgb
fa_name = f_total_matrix
fb_name = f_total_bub
args_a = 'cvm'
args_b = 'cvb'
[../]
#
# Cahn-Hilliard Equations
#
[./CHBulk_g]
type = KKSSplitCHCRes
variable = cg
ca = cgm
fa_name = f_total_matrix
w = wg
args_a = 'cvm'
[../]
[./CHBulk_v]
type = KKSSplitCHCRes
variable = cv
ca = cvm
fa_name = f_total_matrix
w = wv
args_a = 'cgm'
[../]
[./dcgdt]
type = CoupledTimeDerivative
variable = wg
v = cg
[../]
[./dcvdt]
type = CoupledTimeDerivative
variable = wv
v = cv
[../]
[./wgkernel]
type = SplitCHWRes
mob_name = M
variable = wg
[../]
[./wvkernel]
type = SplitCHWRes
mob_name = M
variable = wv
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = f_total_matrix
fb_name = f_total_bub
w = 0.356
args = 'cvm cvb cgm cgb'
[../]
[./ACBulkCv]
type = KKSACBulkC
variable = eta
ca = cvm
cb = cvb
fa_name = f_total_matrix
args = 'cgm'
[../]
[./ACBulkCg]
type = KKSACBulkC
variable = eta
ca = cgm
cb = cgb
fa_name = f_total_matrix
args = 'cvm'
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[Materials]
# Chemical free energy of the matrix
[./fm]
type = DerivativeParsedMaterial
property_name = fm
coupled_variables = 'cvm cgm'
material_property_names = 'kvmatrix kgmatrix cvmatrixeq cgmatrixeq'
expression = '0.5*kvmatrix*(cvm-cvmatrixeq)^2 + 0.5*kgmatrix*(cgm-cgmatrixeq)^2'
[../]
# Elastic energy of the matrix
[./elastic_free_energy_m]
type = ElasticEnergyMaterial
base_name = matrix
f_name = fe_m
args = ' '
[../]
# Total free energy of the matrix
[./Total_energy_matrix]
type = DerivativeSumMaterial
property_name = f_total_matrix
sum_materials = 'fm fe_m'
coupled_variables = 'cvm cgm'
[../]
# Free energy of the bubble phase
[./fb]
type = DerivativeParsedMaterial
property_name = fb
coupled_variables = 'cvb cgb'
material_property_names = 'kToverV nQ Va b f0 kpen kgbub kvbub cvbubeq cgbubeq'
expression = '0.5*kgbub*(cvb-cvbubeq)^2 + 0.5*kvbub*(cgb-cgbubeq)^2'
[../]
# Elastic energy of the bubble
[./elastic_free_energy_p]
type = ElasticEnergyMaterial
base_name = bub
f_name = fe_b
args = ' '
[../]
# Total free energy of the bubble
[./Total_energy_bub]
type = DerivativeSumMaterial
property_name = f_total_bub
sum_materials = 'fb fe_b'
# sum_materials = 'fb'
coupled_variables = 'cvb cgb'
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L kappa Va kvmatrix kgmatrix kgbub kvbub f0 kpen cvbubeq cgbubeq b T'
prop_values = '0.7 0.7 0.0368 0.03629 223.16 223.16 2.23 2.23 0.0224 1.0 0.6076 0.3924 0.085 800'
[../]
[./cvmatrixeq]
type = ParsedMaterial
property_name = cvmatrixeq
material_property_names = 'T'
constant_names = 'kB Efv'
constant_expressions = '8.6173324e-5 1.69'
expression = 'exp(-Efv/(kB*T))'
[../]
[./cgmatrixeq]
type = ParsedMaterial
property_name = cgmatrixeq
material_property_names = 'T'
constant_names = 'kB Efg'
constant_expressions = '8.6173324e-5 4.92'
expression = 'exp(-Efg/(kB*T))'
[../]
[./kToverV]
type = ParsedMaterial
property_name = kToverV
material_property_names = 'T Va'
constant_names = 'k C44dim' #k in J/K and dimensional C44 in J/m^3
constant_expressions = '1.38e-23 63e9'
expression = 'k*T*1e27/Va/C44dim'
[../]
[./nQ]
type = ParsedMaterial
property_name = nQ
material_property_names = 'T'
constant_names = 'k Pi M hbar' #k in J/K, M is Xe atomic mass in kg, hbar in J s
constant_expressions = '1.38e-23 3.14159 2.18e-25 1.05459e-34'
expression = '(M*k*T/2/Pi/hbar^2)^1.5 * 1e-27' #1e-27 converts from #/m^3 to #/nm^3
[../]
#Mechanical properties
[./Stiffness_matrix]
type = ComputeElasticityTensor
C_ijkl = '0.778 0.7935'
fill_method = symmetric_isotropic
base_name = matrix
[../]
[./Stiffness_bub]
type = ComputeElasticityTensor
C_ijkl = '0.0778 0.07935'
fill_method = symmetric_isotropic
base_name = bub
[../]
[./strain_matrix]
type = ComputeRSphericalSmallStrain
base_name = matrix
[../]
[./strain_bub]
type = ComputeRSphericalSmallStrain
base_name = bub
[../]
[./stress_matrix]
type = ComputeLinearElasticStress
base_name = matrix
[../]
[./stress_bub]
type = ComputeLinearElasticStress
base_name = bub
[../]
[./global_stress]
type = TwoPhaseStressMaterial
base_A = matrix
base_B = bub
[../]
[./surface_tension]
type = ComputeSurfaceTensionKKS
v = eta
kappa_name = kappa
w = 0.356
[../]
[./gas_pressure]
type = ComputeExtraStressVDWGas
T = T
b = b
cg = cgb
Va = Va
nondim_factor = 63e9
base_name = bub
outputs = exodus
[../]
[]
[BCs]
[./left_r]
type = DirichletBC
variable = disp_x
boundary = left
value = 0
[../]
[]
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'asm lu nonzero'
l_max_its = 30
nl_max_its = 15
l_tol = 1.0e-4
nl_rel_tol = 1.0e-10
nl_abs_tol = 1e-11
num_steps = 2
dt = 0.5
[]
[Outputs]
exodus = true
[]
(modules/combined/examples/publications/rapid_dev/fig7a.i)
#
# Fig. 7 input for 10.1016/j.commatsci.2017.02.017
# D. Schwen et al./Computational Materials Science 132 (2017) 36-45
# Solid gray curve (1)
# Eigenstrain and elastic energies ar computed per phase and then interpolated.
# Supply the RADIUS parameter (10-35) on the command line to generate data
# for all curves in the plot.
#
[Mesh]
type = GeneratedMesh
dim = 1
nx = 32
xmin = 0
xmax = 100
second_order = true
[]
[Problem]
coord_type = RSPHERICAL
[]
[GlobalParams]
displacements = 'disp_r'
[]
[Functions]
[./diff]
type = ParsedFunction
expression = '${RADIUS}-pos_c'
symbol_names = pos_c
symbol_values = pos_c
[../]
[]
# AuxVars to compute the free energy density for outputting
[AuxVariables]
[./local_energy]
order = CONSTANT
family = MONOMIAL
[../]
[./cross_energy]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./local_free_energy]
type = TotalFreeEnergy
variable = local_energy
interfacial_vars = 'c'
kappa_names = 'kappa_c'
execute_on = 'INITIAL TIMESTEP_END'
[../]
[]
[Variables]
# Solute concentration variable
[./c]
[./InitialCondition]
type = SmoothCircleIC
invalue = 1
outvalue = 0
x1 = 0
y1 = 0
radius = ${RADIUS}
int_width = 3
[../]
[../]
[./w]
[../]
# Phase order parameter
[./eta]
[./InitialCondition]
type = SmoothCircleIC
invalue = 1
outvalue = 0
x1 = 0
y1 = 0
radius = ${RADIUS}
int_width = 3
[../]
[../]
# Mesh displacement
[./disp_r]
order = SECOND
[../]
[./Fe_fit]
order = SECOND
[../]
[]
[Kernels]
# Set up stress divergence kernels
[./TensorMechanics]
[../]
# Split Cahn-Hilliard kernels
[./c_res]
type = SplitCHParsed
variable = c
f_name = F
args = 'eta'
kappa_name = kappa_c
w = w
[../]
[./wres]
type = SplitCHWRes
variable = w
mob_name = M
[../]
[./time]
type = CoupledTimeDerivative
variable = w
v = c
[../]
# Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 1
[./detadt]
type = TimeDerivative
variable = eta
[../]
[./ACBulk1]
type = AllenCahn
variable = eta
args = 'c'
mob_name = L
f_name = F
[../]
[./ACInterface]
type = ACInterface
variable = eta
mob_name = L
kappa_name = kappa_eta
[../]
[./Fe]
type = MaterialPropertyValue
prop_name = Fe
variable = Fe_fit
[../]
[./autoadjust]
type = MaskedBodyForce
variable = w
function = diff
mask = mask
[../]
[]
[Materials]
# declare a few constants, such as mobilities (L,M) and interface gradient prefactors (kappa*)
[./consts]
type = GenericConstantMaterial
prop_names = 'M L kappa_c kappa_eta'
prop_values = '1.0 1.0 0.5 1'
[../]
# forcing function mask
[./mask]
type = ParsedMaterial
property_name = mask
expression = grad/dt
material_property_names = 'grad dt'
[../]
[./grad]
type = VariableGradientMaterial
variable = c
prop = grad
[../]
[./time]
type = TimeStepMaterial
[../]
# global mechanical properties
[./elasticity_tensor_1]
type = ComputeElasticityTensor
C_ijkl = '1 1'
base_name = phase1
fill_method = symmetric_isotropic
[../]
[./elasticity_tensor_2]
type = ComputeElasticityTensor
C_ijkl = '1 1'
base_name = phase2
fill_method = symmetric_isotropic
[../]
[./strain_1]
type = ComputeRSphericalSmallStrain
base_name = phase1
[../]
[./strain_2]
type = ComputeRSphericalSmallStrain
base_name = phase2
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
base_name = phase1
[../]
[./stress_2]
type = ComputeLinearElasticStress
base_name = phase2
[../]
# eigenstrain per phase
[./eigenstrain2]
type = ComputeEigenstrain
eigen_base = '0.05 0.05 0.05 0 0 0'
base_name = phase2
eigenstrain_name = eigenstrain
[../]
# switching functions
[./switching]
type = SwitchingFunctionMaterial
function_name = h
eta = eta
h_order = SIMPLE
[../]
[./barrier]
type = BarrierFunctionMaterial
eta = eta
[../]
# chemical free energies
[./chemical_free_energy_1]
type = DerivativeParsedMaterial
property_name = Fc1
expression = 'c^2'
coupled_variables = 'c'
derivative_order = 2
[../]
[./chemical_free_energy_2]
type = DerivativeParsedMaterial
property_name = Fc2
expression = '(1-c)^2'
coupled_variables = 'c'
derivative_order = 2
[../]
# elastic free energies
[./elastic_free_energy_1]
type = ElasticEnergyMaterial
f_name = Fe1
args = ''
base_name = phase1
derivative_order = 2
[../]
[./elastic_free_energy_2]
type = ElasticEnergyMaterial
f_name = Fe2
args = ''
base_name = phase2
derivative_order = 2
[../]
# per phase free energies
[./free_energy_1]
type = DerivativeSumMaterial
property_name = F1
sum_materials = 'Fc1 Fe1'
coupled_variables = 'c'
derivative_order = 2
[../]
[./free_energy_2]
type = DerivativeSumMaterial
property_name = F2
sum_materials = 'Fc2 Fe2'
coupled_variables = 'c'
derivative_order = 2
[../]
# global chemical free energy
[./global_free_energy]
type = DerivativeTwoPhaseMaterial
f_name = F
fa_name = F1
fb_name = F2
eta = eta
args = 'c'
W = 4
[../]
# global stress
[./global_stress]
type = TwoPhaseStressMaterial
base_A = phase1
base_B = phase2
[../]
[./elastic_free_energy]
type = DerivativeTwoPhaseMaterial
f_name = Fe
fa_name = Fe1
fb_name = Fe2
eta = eta
args = 'c'
W = 0
[../]
[]
[BCs]
[./left_r]
type = DirichletBC
variable = disp_r
boundary = 'left'
value = 0
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
# We monitor the total free energy and the total solute concentration (should be constant)
[Postprocessors]
[./total_free_energy]
type = ElementIntegralVariablePostprocessor
variable = local_energy
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[./total_solute]
type = ElementIntegralVariablePostprocessor
variable = c
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[./pos_c]
type = FindValueOnLine
start_point = '0 0 0'
end_point = '100 0 0'
v = c
target = 0.582
tol = 1e-8
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[./pos_eta]
type = FindValueOnLine
start_point = '0 0 0'
end_point = '100 0 0'
v = eta
target = 0.5
tol = 1e-8
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[./c_min]
type = ElementExtremeValue
value_type = min
variable = c
execute_on = 'INITIAL TIMESTEP_END'
outputs = 'table console'
[../]
[]
[VectorPostprocessors]
[./line]
type = LineValueSampler
variable = 'Fe_fit c w'
start_point = '0 0 0'
end_point = '100 0 0'
num_points = 5000
sort_by = x
outputs = vpp
execute_on = 'INITIAL TIMESTEP_END'
[../]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type'
petsc_options_value = 'asm lu'
l_max_its = 30
nl_max_its = 15
l_tol = 1.0e-4
nl_rel_tol = 1.0e-8
nl_abs_tol = 2.0e-9
start_time = 0.0
end_time = 100000.0
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 7
iteration_window = 1
dt = 1
[../]
[./Adaptivity]
initial_adaptivity = 5
interval = 10
max_h_level = 5
refine_fraction = 0.9
coarsen_fraction = 0.1
[../]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
execute_on = 'INITIAL TIMESTEP_END'
[./table]
type = CSV
delimiter = ' '
file_base = radius_${RADIUS}/energy_pp
[../]
[./vpp]
type = CSV
delimiter = ' '
sync_times = '10 50 100 500 1000 5000 10000 50000 100000'
sync_only = true
time_data = true
file_base = radius_${RADIUS}/energy_vpp
[../]
[]
(modules/phase_field/test/tests/KKS_system/kks_example.i)
#
# KKS toy problem in the non-split form
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 5
ny = 5
nz = 0
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[Variables]
# order parameter
[./eta]
order = THIRD
family = HERMITE
[../]
# hydrogen concentration
[./c]
order = THIRD
family = HERMITE
[../]
# hydrogen phase concentration (matrix)
[./cm]
order = THIRD
family = HERMITE
initial_condition = 0.0
[../]
# hydrogen phase concentration (delta phase)
[./cd]
order = THIRD
family = HERMITE
initial_condition = 0.0
[../]
[]
[ICs]
[./eta]
variable = eta
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 0.2
invalue = 0.2
outvalue = 0.1
int_width = 0.05
[../]
[./c]
variable = c
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 0.2
invalue = 0.6
outvalue = 0.4
int_width = 0.05
[../]
[]
[BCs]
[./Periodic]
[./all]
variable = 'eta c cm cd'
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
# Free energy of the matrix
[./fm]
type = DerivativeParsedMaterial
property_name = fm
coupled_variables = 'cm'
expression = '(0.1-cm)^2'
outputs = oversampling
[../]
# Free energy of the delta phase
[./fd]
type = DerivativeParsedMaterial
property_name = fd
coupled_variables = 'cd'
expression = '(0.9-cd)^2'
outputs = oversampling
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
outputs = oversampling
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
outputs = oversampling
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'L '
prop_values = '0.7 '
[../]
[]
[Kernels]
# enforce c = (1-h(eta))*cm + h(eta)*cd
[./PhaseConc]
type = KKSPhaseConcentration
ca = cm
variable = cd
c = c
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotVacancies]
type = KKSPhaseChemicalPotential
variable = cm
cb = cd
fa_name = fm
fb_name = fd
[../]
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = KKSCHBulk
variable = c
ca = cm
cb = cd
fa_name = fm
fb_name = fd
mob_name = 0.7
[../]
[./dcdt]
type = TimeDerivative
variable = c
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = fm
fb_name = fd
coupled_variables = 'cm cd'
w = 0.4
[../]
[./ACBulkC]
type = KKSACBulkC
variable = eta
ca = cm
cb = cd
fa_name = fm
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = 0.4
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pctype -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = ' asm lu nonzero'
l_max_its = 100
nl_max_its = 100
nl_rel_tol = 1e-4
num_steps = 1
dt = 0.01
dtmin = 0.01
[]
[Preconditioning]
[./mydebug]
type = SMP
full = true
[../]
[]
[Outputs]
file_base = kks_example
[./oversampling]
type = Exodus
refinements = 3
[../]
[]
(modules/phase_field/test/tests/MultiPhase/asymmetriccrosstermbarrierfunction.i)
[Mesh]
type = GeneratedMesh
dim = 1
nx = 200
xmin = 0
xmax = 9
[]
[Functions]
[./func1]
type = ParsedFunction
expression = 'il:=x-7; ir:=2-x; if(x<1, 1,
if(x<2, 0.5-0.5*cos(ir*pi),
if(x<7, 0,
if(x<8, 0.5-0.5*cos(il*pi),
1))))'
[../]
[./func2]
type = ParsedFunction
expression = 'il:=x-1; ir:=5-x; if(x<1, 0,
if(x<2, 0.5-0.5*cos(il*pi),
if(x<4, 1,
if(x<5, 0.5-0.5*cos(ir*pi),
0))))'
[../]
[./func3]
type = ParsedFunction
expression = 'il:=x-4; ir:=8-x; if(x<4, 0,
if(x<5, 0.5-0.5*cos(il*pi),
if(x<7, 1,
if(x<8, 0.5-0.5*cos(ir*pi),
0))))'
[../]
[]
[AuxVariables]
[./eta1]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
function = func1
[../]
[../]
[./eta2]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
function = func2
[../]
[../]
[./eta3]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
function = func3
[../]
[../]
[]
[Materials]
[./symmetriccrosstermbarrier_low]
type = AsymmetricCrossTermBarrierFunctionMaterial
etas = 'eta1 eta2 eta3'
hi_names = 'h1 h2 h3'
W_ij = '0 1 2.2
1 0 3.1
2.2 3.1 0'
function_name = gsl
g_order = LOW
outputs = exodus
[../]
[./asymmetriccrosstermbarrier_low]
type = AsymmetricCrossTermBarrierFunctionMaterial
etas = 'eta1 eta2 eta3'
hi_names = 'h1 h2 h3'
W_ij = ' 0 1.2 5.2
0.8 0 2.1
-0.8 4.1 0'
function_name = gal
g_order = LOW
outputs = exodus
[../]
[./asymmetriccrosstermbarrie_simple]
type = AsymmetricCrossTermBarrierFunctionMaterial
etas = 'eta1 eta2 eta3'
hi_names = 'h1 h2 h3'
W_ij = '0 1.2 3.2
0.8 0 2.1
1.2 4.1 0'
function_name = gas
g_order = SIMPLE
outputs = exodus
[../]
[./switch1]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
[../]
[./switch2]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
[../]
[./switch3]
type = SwitchingFunctionMaterial
function_name = h3
eta = eta3
[../]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
num_steps = 1
[]
[Problem]
solve = false
kernel_coverage_check = false
[]
[Outputs]
exodus = true
execute_on = final
[]
(modules/phase_field/test/tests/MultiPhase/acmultiinterface.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 20
ny = 10
nz = 0
xmin = -10
xmax = 10
ymin = -5
ymax = 5
elem_type = QUAD4
[]
[Variables]
[./eta1]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = -3.5
y1 = 0.0
radius = 4.0
invalue = 0.9
outvalue = 0.1
int_width = 2.0
[../]
[../]
[./eta2]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 3.5
y1 = 0.0
radius = 4.0
invalue = 0.9
outvalue = 0.1
int_width = 2.0
[../]
[../]
[./eta3]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SpecifiedSmoothCircleIC
x_positions = '-4.0 4.0'
y_positions = ' 0.0 0.0'
z_positions = ' 0.0 0.0'
radii = '4.0 4.0'
invalue = 0.1
outvalue = 0.9
int_width = 2.0
[../]
[../]
[./lambda]
order = FIRST
family = LAGRANGE
initial_condition = 1.0
[../]
[]
[Kernels]
[./deta1dt]
type = TimeDerivative
variable = eta1
[../]
[./ACBulk1]
type = AllenCahn
variable = eta1
coupled_variables = 'eta2 eta3'
mob_name = L1
f_name = F
[../]
[./ACInterface1]
type = ACMultiInterface
variable = eta1
etas = 'eta1 eta2 eta3'
mob_name = L1
kappa_names = 'kappa11 kappa12 kappa13'
[../]
[./lagrange1]
type = SwitchingFunctionConstraintEta
variable = eta1
h_name = h1
lambda = lambda
[../]
[./deta2dt]
type = TimeDerivative
variable = eta2
[../]
[./ACBulk2]
type = AllenCahn
variable = eta2
coupled_variables = 'eta1 eta3'
mob_name = L2
f_name = F
[../]
[./ACInterface2]
type = ACMultiInterface
variable = eta2
etas = 'eta1 eta2 eta3'
mob_name = L2
kappa_names = 'kappa21 kappa22 kappa23'
[../]
[./lagrange2]
type = SwitchingFunctionConstraintEta
variable = eta2
h_name = h2
lambda = lambda
[../]
[./deta3dt]
type = TimeDerivative
variable = eta3
[../]
[./ACBulk3]
type = AllenCahn
variable = eta3
coupled_variables = 'eta1 eta2'
mob_name = L3
f_name = F
[../]
[./ACInterface3]
type = ACMultiInterface
variable = eta3
etas = 'eta1 eta2 eta3'
mob_name = L3
kappa_names = 'kappa31 kappa32 kappa33'
[../]
[./lagrange3]
type = SwitchingFunctionConstraintEta
variable = eta3
h_name = h3
lambda = lambda
[../]
[./lagrange]
type = SwitchingFunctionConstraintLagrange
variable = lambda
etas = 'eta1 eta2 eta3'
h_names = 'h1 h2 h3'
epsilon = 0
[../]
[]
[BCs]
[./Periodic]
[./All]
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
[./consts]
type = GenericConstantMaterial
prop_names = 'Fx L1 L2 L3 kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
prop_values = '0 1 1 1 1 1 1 1 1 1 1 1 1 '
[../]
[./etasummat]
type = ParsedMaterial
property_name = etasum
coupled_variables = 'eta1 eta2 eta3'
material_property_names = 'h1 h2 h3'
expression = 'h1+h2+h3'
[../]
[./switching1]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
h_order = SIMPLE
[../]
[./switching2]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
h_order = SIMPLE
[../]
[./switching3]
type = SwitchingFunctionMaterial
function_name = h3
eta = eta3
h_order = SIMPLE
[../]
[./barrier]
type = MultiBarrierFunctionMaterial
etas = 'eta1 eta2 eta3'
[../]
[./free_energy]
type = DerivativeMultiPhaseMaterial
property_name = F
# we use a constant free energy (GeneriConstantmaterial property Fx)
fi_names = 'Fx Fx Fx'
hi_names = 'h1 h2 h3'
etas = 'eta1 eta2 eta3'
# the free energy is given by the MultiBarrierFunctionMaterial only
W = 1
derivative_order = 2
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = 'bdf2'
solve_type = 'PJFNK'
#petsc_options = '-snes_ksp -snes_ksp_ew'
#petsc_options = '-ksp_monitor_snes_lg-snes_ksp_ew'
#petsc_options_iname = '-ksp_gmres_restart'
#petsc_options_value = '1000 '
l_max_its = 15
l_tol = 1.0e-6
nl_max_its = 50
nl_rel_tol = 1.0e-8
nl_abs_tol = 1.0e-10
start_time = 0.0
num_steps = 2
dt = 0.2
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(modules/phase_field/test/tests/GrandPotentialPFM/GrandPotentialPFM.i)
# this input file test the implementation of the grand-potential phase-field model based on M.Plapp PRE 84,031601(2011)
# in this simple example, the liquid and solid free energies are parabola with the same curvature and the material properties are constant
# Note that this example also test The SusceptibilityTimeDerivative kernels
[Mesh]
type = GeneratedMesh
dim = 2
nx = 16
ny = 16
xmax = 32
ymax = 32
[]
[GlobalParams]
radius = 20.0
int_width = 4.0
x1 = 0
y1 = 0
[]
[Variables]
[./w]
[../]
[./eta]
[../]
[]
[ICs]
[./w]
type = SmoothCircleIC
variable = w
# note w = A*(c-cleq), A = 1.0, cleq = 0.0 ,i.e., w = c (in the matrix/liquid phase)
outvalue = -0.2
invalue = 0.2
[../]
[./eta]
type = SmoothCircleIC
variable = eta
outvalue = 0.0
invalue = 1.0
[../]
[]
[Kernels]
[./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 = D
args = ''
[../]
[./coupled_etadot]
type = CoupledSusceptibilityTimeDerivative
variable = w
v = eta
f_name = ft
coupled_variables = 'eta'
[../]
[./AC_bulk]
type = AllenCahn
variable = eta
f_name = F
coupled_variables = 'w'
[../]
[./AC_int]
type = ACInterface
variable = eta
[../]
[./e_dot]
type = TimeDerivative
variable = eta
[../]
[]
[Materials]
[./constants]
type = GenericConstantMaterial
prop_names = 'kappa_op D L chi cseq cleq A'
prop_values = '4.0 1.0 1.0 1.0 0.0 1.0 1.0'
[../]
[./liquid_GrandPotential]
type = DerivativeParsedMaterial
expression = '-0.5 * w^2/A - cleq * w'
coupled_variables = 'w'
property_name = f1
material_property_names = 'cleq A'
[../]
[./solid_GrandPotential]
type = DerivativeParsedMaterial
expression = '-0.5 * w^2/A - cseq * w'
coupled_variables = 'w'
property_name = f2
material_property_names = 'cseq A'
[../]
[./switching_function]
type = SwitchingFunctionMaterial
eta = eta
h_order = HIGH
[../]
[./barrier_function]
type = BarrierFunctionMaterial
eta = eta
[../]
[./cs]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = cs
material_property_names = 'A cseq'
expression = 'w/A + cseq' # since w = A*(c-cseq)
derivative_order = 2
[../]
[./cl]
type = DerivativeParsedMaterial
coupled_variables = 'w'
property_name = cl
material_property_names = 'A cleq'
expression = 'w/A + cleq' # since w = A*(c-cleq)
derivative_order = 2
[../]
[./total_GrandPotential]
type = DerivativeTwoPhaseMaterial
coupled_variables = 'w'
eta = eta
fa_name = f1
fb_name = f2
derivative_order = 2
W = 1.0
[../]
[./coupled_eta_function]
type = DerivativeParsedMaterial
expression = '(cs - cl) * dh'
coupled_variables = 'eta w'
property_name = ft
material_property_names = 'cs cl dh:=D[h,eta]'
derivative_order = 1
outputs = exodus
[../]
[./concentration]
type = ParsedMaterial
property_name = c
material_property_names = 'dF:=D[F,w]'
expression = '-dF'
outputs = exodus
[../]
[]
[Postprocessors]
[./C]
type = ElementIntegralMaterialProperty
mat_prop = c
execute_on = 'initial timestep_end'
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = NEWTON
l_max_its = 15
l_tol = 1e-3
nl_max_its = 15
nl_rel_tol = 1e-8
nl_abs_tol = 1e-8
num_steps = 5
dt = 10.0
[]
[Outputs]
exodus = true
csv = true
execute_on = 'TIMESTEP_END'
[]
(modules/phase_field/examples/multiphase/DerivativeMultiPhaseMaterial.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 40
ny = 40
nz = 0
xmin = -12
xmax = 12
ymin = -12
ymax = 12
elem_type = QUAD4
[]
[GlobalParams]
# let's output all material properties for demonstration purposes
outputs = exodus
# prefactor on the penalty function kernels. The higher this value is, the
# more rigorously the constraint is enforced
penalty = 1e3
[]
#
# These AuxVariables hold the directly calculated free energy density in the
# simulation cell. They are provided for visualization purposes.
#
[AuxVariables]
[./local_energy]
order = CONSTANT
family = MONOMIAL
[../]
[./cross_energy]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./local_free_energy]
type = TotalFreeEnergy
variable = local_energy
interfacial_vars = 'c'
kappa_names = 'kappa_c'
additional_free_energy = cross_energy
[../]
#
# Helper kernel to cpompute the gradient contribution from interfaces of order
# parameters evolved using the ACMultiInterface kernel
#
[./cross_terms]
type = CrossTermGradientFreeEnergy
variable = cross_energy
interfacial_vars = 'eta1 eta2 eta3'
#
# The interface coefficient matrix. This should be symmetrical!
#
kappa_names = 'kappa11 kappa12 kappa13
kappa21 kappa22 kappa23
kappa31 kappa32 kappa33'
[../]
[]
[Variables]
[./c]
order = FIRST
family = LAGRANGE
#
# We set up a smooth cradial concentrtaion gradient
# The concentration will quickly change to adapt to the preset order
# parameters eta1, eta2, and eta3
#
[./InitialCondition]
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 5.0
invalue = 1.0
outvalue = 0.01
int_width = 10.0
[../]
[../]
[./eta1]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
#
# Note: this initial conditions sets up a _sharp_ interface. Ideally
# we should start with a smooth interface with a width consistent
# with the kappa parameter supplied for the given interface.
#
function = 'r:=sqrt(x^2+y^2);if(r<=4,1,0)'
[../]
[../]
[./eta2]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
function = 'r:=sqrt(x^2+y^2);if(r>4&r<=7,1,0)'
[../]
[../]
[./eta3]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
function = 'r:=sqrt(x^2+y^2);if(r>7,1,0)'
[../]
[../]
[]
[Kernels]
#
# Cahn-Hilliard kernel for the concentration variable.
# Note that we are not using an interfcae kernel on this variable, but rather
# rely on the interface width enforced on the order parameters. This allows us
# to use a direct solve using the CahnHilliard kernel _despite_ only using first
# order elements.
#
[./c_res]
type = CahnHilliard
variable = c
f_name = F
coupled_variables = 'eta1 eta2 eta3'
[../]
[./time]
type = TimeDerivative
variable = c
[../]
#
# Order parameter eta1
# Each order parameter is acted on by 4 kernels:
# 1. The stock time derivative deta_i/dt kernel
# 2. The Allen-Cahn kernel that takes a Dervative Material for the free energy
# 3. A gradient interface kernel that includes cross terms
# see http://mooseframework.org/wiki/PhysicsModules/PhaseField/DevelopingModels/MultiPhaseModels/ACMultiInterface/
# 4. A penalty contribution that forces the interface contributions h(eta)
# to sum up to unity
#
[./deta1dt]
type = TimeDerivative
variable = eta1
[../]
[./ACBulk1]
type = AllenCahn
variable = eta1
coupled_variables = 'eta2 eta3 c'
mob_name = L1
f_name = F
[../]
[./ACInterface1]
type = ACMultiInterface
variable = eta1
etas = 'eta1 eta2 eta3'
mob_name = L1
kappa_names = 'kappa11 kappa12 kappa13'
[../]
[./penalty1]
type = SwitchingFunctionPenalty
variable = eta1
etas = 'eta1 eta2 eta3'
h_names = 'h1 h2 h3'
[../]
#
# Order parameter eta2
#
[./deta2dt]
type = TimeDerivative
variable = eta2
[../]
[./ACBulk2]
type = AllenCahn
variable = eta2
coupled_variables = 'eta1 eta3 c'
mob_name = L2
f_name = F
[../]
[./ACInterface2]
type = ACMultiInterface
variable = eta2
etas = 'eta1 eta2 eta3'
mob_name = L2
kappa_names = 'kappa21 kappa22 kappa23'
[../]
[./penalty2]
type = SwitchingFunctionPenalty
variable = eta2
etas = 'eta1 eta2 eta3'
h_names = 'h1 h2 h3'
[../]
#
# Order parameter eta3
#
[./deta3dt]
type = TimeDerivative
variable = eta3
[../]
[./ACBulk3]
type = AllenCahn
variable = eta3
coupled_variables = 'eta1 eta2 c'
mob_name = L3
f_name = F
[../]
[./ACInterface3]
type = ACMultiInterface
variable = eta3
etas = 'eta1 eta2 eta3'
mob_name = L3
kappa_names = 'kappa31 kappa32 kappa33'
[../]
[./penalty3]
type = SwitchingFunctionPenalty
variable = eta3
etas = 'eta1 eta2 eta3'
h_names = 'h1 h2 h3'
[../]
[]
[BCs]
[./Periodic]
[./All]
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
# here we declare some of the model parameters: the mobilities and interface
# gradient prefactors. For this example we use arbitrary numbers. In an actual simulation
# physical mobilities would be used, and the interface gradient prefactors would
# be readjusted to the free energy magnitudes.
[./consts]
type = GenericConstantMaterial
prop_names = 'M kappa_c L1 L2 L3 kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
prop_values = '0.2 0.75 1 1 1 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 '
[../]
# This material sums up the individual phase contributions. It is written to the output file
# (see GlobalParams section above) and can be used to check the constraint enforcement.
[./etasummat]
type = ParsedMaterial
property_name = etasum
material_property_names = 'h1 h2 h3'
expression = 'h1+h2+h3'
[../]
# The phase contribution factors for each material point are computed using the
# SwitchingFunctionMaterials. Each phase with an order parameter eta contributes h(eta)
# to the global free energy density. h is a function that switches smoothly from 0 to 1
[./switching1]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
h_order = SIMPLE
[../]
[./switching2]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
h_order = SIMPLE
[../]
[./switching3]
type = SwitchingFunctionMaterial
function_name = h3
eta = eta3
h_order = SIMPLE
[../]
# The barrier function adds a phase transformation energy barrier. It also
# Drives order parameters toward the [0:1] interval to avoid negative or larger than 1
# order parameters (these are set to 0 and 1 contribution by the switching functions
# above)
[./barrier]
type = MultiBarrierFunctionMaterial
etas = 'eta1 eta2 eta3'
[../]
# We use DerivativeParsedMaterials to specify three (very) simple free energy
# expressions for the three phases. All necessary derivatives are built automatically.
# In a real problem these expressions can be arbitrarily complex (or even provided
# by custom kernels).
[./phase_free_energy_1]
type = DerivativeParsedMaterial
property_name = F1
expression = '(c-1)^2'
coupled_variables = 'c'
[../]
[./phase_free_energy_2]
type = DerivativeParsedMaterial
property_name = F2
expression = '(c-0.5)^2'
coupled_variables = 'c'
[../]
[./phase_free_energy_3]
type = DerivativeParsedMaterial
property_name = F3
expression = 'c^2'
coupled_variables = 'c'
[../]
# The DerivativeMultiPhaseMaterial ties the phase free energies together into a global free energy.
# http://mooseframework.org/wiki/PhysicsModules/PhaseField/DevelopingModels/MultiPhaseModels/
[./free_energy]
type = DerivativeMultiPhaseMaterial
property_name = F
# we use a constant free energy (GeneriConstantmaterial property Fx)
fi_names = 'F1 F2 F3'
hi_names = 'h1 h2 h3'
etas = 'eta1 eta2 eta3'
coupled_variables = 'c'
W = 1
[../]
[]
[Postprocessors]
# The total free energy of the simulation cell to observe the energy reduction.
[./total_free_energy]
type = ElementIntegralVariablePostprocessor
variable = local_energy
[../]
# for testing we also monitor the total solute amount, which should be conserved.
[./total_solute]
type = ElementIntegralVariablePostprocessor
variable = c
[../]
[]
[Preconditioning]
# This preconditioner makes sure the Jacobian Matrix is fully populated. Our
# kernels compute all Jacobian matrix entries.
# This allows us to use the Newton solver below.
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = 'bdf2'
# Automatic differentiation provedes a _full_ Jacobian in this example
# so we can safely use NEWTON for a fast solve
solve_type = 'NEWTON'
l_max_its = 15
l_tol = 1.0e-6
nl_max_its = 50
nl_rel_tol = 1.0e-6
nl_abs_tol = 1.0e-6
start_time = 0.0
end_time = 150.0
[./TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 0.1
[../]
[]
[Debug]
# show_var_residual_norms = true
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[./table]
type = CSV
delimiter = ' '
[../]
[]
(modules/phase_field/test/tests/KKS_system/nonlinear.i)
#
# This test checks if the thwo phase and lagrange multiplier solutions can be replicated
# with a two order parameter approach, where the second order parameter eta2 is a
# nonlinear variable that is set as eta2 := 1 - eta1 (using Reaction, CoupledForce, and BodyForce)
# The solution is reproduced.
#
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 5
[]
[AuxVariables]
[Fglobal]
order = CONSTANT
family = MONOMIAL
[]
[]
[Variables]
# concentration
[c]
order = FIRST
family = LAGRANGE
[InitialCondition]
type = FunctionIC
function = x/5
[]
[]
# order parameter 1
[eta1]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
# order parameter 2
[eta2]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
# phase concentration 1
[c1]
order = FIRST
family = LAGRANGE
initial_condition = 0.9
[]
# phase concentration 2
[c2]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[]
[]
[Materials]
# simple toy free energies
[f1] # = fd
type = DerivativeParsedMaterial
property_name = F1
coupled_variables = 'c1'
expression = '(0.9-c1)^2'
[]
[f2] # = fm
type = DerivativeParsedMaterial
property_name = F2
coupled_variables = 'c2'
expression = '(0.1-c2)^2'
[]
# Switching functions for each phase
[h1_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta1
function_name = h1
[]
[h2_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta2
function_name = h2
[]
# Coefficients for diffusion equation
[Dh1]
type = DerivativeParsedMaterial
material_property_names = 'D h1(eta1)'
expression = D*h1
property_name = Dh1
coupled_variables = eta1
[]
[Dh2]
type = DerivativeParsedMaterial
material_property_names = 'D h2(eta2)'
expression = D*h2
property_name = Dh2
coupled_variables = eta2
[]
# Barrier functions for each phase
[g1]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta1
function_name = g1
[]
[g2]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta2
function_name = g2
[]
# constant properties
[constants]
type = GenericConstantMaterial
prop_names = 'D L kappa'
prop_values = '0.7 0.7 0.2'
[]
[]
[Kernels]
#Kernels for diffusion equation
[diff_time]
type = TimeDerivative
variable = c
[]
[diff_c1]
type = MatDiffusion
variable = c
diffusivity = Dh1
v = c1
args = 'eta1'
[]
[diff_c2]
type = MatDiffusion
variable = c
diffusivity = Dh2
v = c2
args = 'eta2'
[]
# Kernels for Allen-Cahn equation for eta1
[deta1dt]
type = TimeDerivative
variable = eta1
[]
[ACBulkF1]
type = KKSMultiACBulkF
variable = eta1
Fj_names = 'F1 F2 '
hj_names = 'h1 h2 '
gi_name = g1
eta_i = eta1
wi = 0.2
coupled_variables = 'c1 c2 eta2'
[]
[ACBulkC1]
type = KKSMultiACBulkC
variable = eta1
Fj_names = 'F1 F2'
hj_names = 'h1 h2'
cj_names = 'c1 c2'
eta_i = eta1
coupled_variables = 'eta2'
[]
[ACInterface1]
type = ACInterface
variable = eta1
kappa_name = kappa
[]
# Phase concentration constraints
[chempot12]
type = KKSPhaseChemicalPotential
variable = c1
cb = c2
fa_name = F1
fb_name = F2
[]
[phaseconcentration]
type = KKSMultiPhaseConcentration
variable = c2
cj = 'c1 c2'
hj_names = 'h1 h2'
etas = 'eta1 eta2'
c = c
[]
# equation for eta2 = 1 - eta1
# 0 = eta2 + eta1 -1
[constraint_eta1] # eta2
type = Reaction
variable = eta2
[]
[constraint_eta2] # + eta1
type = CoupledForce
variable = eta2
coef = -1
v = eta1
[]
[constraint_one] # - 1
type = BodyForce
variable = eta2
[]
[]
[AuxKernels]
[Fglobal_total]
type = KKSMultiFreeEnergy
Fj_names = 'F1 F2 '
hj_names = 'h1 h2 '
gj_names = 'g1 g2 '
variable = Fglobal
w = 0.2
interfacial_vars = 'eta1 eta2 '
kappa_names = 'kappa kappa'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'lu nonzero'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-10
nl_abs_tol = 1.0e-11
end_time = 350
dt = 10
[]
[VectorPostprocessors]
[c]
type = LineValueSampler
variable = c
start_point = '0 0 0'
end_point = '5 0 0'
num_points = 21
sort_by = x
[]
[]
[Outputs]
csv = true
execute_on = FINAL
[]
(modules/phase_field/examples/kim-kim-suzuki/kks_example_noflux.i)
#
# KKS simple example in the split form
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 150
ny = 15
nz = 0
xmin = -25
xmax = 25
ymin = -2.5
ymax = 2.5
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[AuxVariables]
[./Fglobal]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# solute concentration
[./c]
order = FIRST
family = LAGRANGE
[../]
# chemical potential
[./w]
order = FIRST
family = LAGRANGE
[../]
# Liquid phase solute concentration
[./cl]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[../]
# Solid phase solute concentration
[./cs]
order = FIRST
family = LAGRANGE
initial_condition = 0.9
[../]
[]
[Functions]
[./ic_func_eta]
type = ParsedFunction
expression = '0.5*(1.0-tanh((x)/sqrt(2.0)))'
[../]
[./ic_func_c]
type = ParsedFunction
expression = '0.9*(0.5*(1.0-tanh(x/sqrt(2.0))))^3*(6*(0.5*(1.0-tanh(x/sqrt(2.0))))^2-15*(0.5*(1.0-tanh(x/sqrt(2.0))))+10)+0.1*(1-(0.5*(1.0-tanh(x/sqrt(2.0))))^3*(6*(0.5*(1.0-tanh(x/sqrt(2.0))))^2-15*(0.5*(1.0-tanh(x/sqrt(2.0))))+10))'
[../]
[]
[ICs]
[./eta]
variable = eta
type = FunctionIC
function = ic_func_eta
[../]
[./c]
variable = c
type = FunctionIC
function = ic_func_c
[../]
[]
[Materials]
# Free energy of the liquid
[./fl]
type = DerivativeParsedMaterial
property_name = fl
coupled_variables = 'cl'
expression = '(0.1-cl)^2'
[../]
# Free energy of the solid
[./fs]
type = DerivativeParsedMaterial
property_name = fs
coupled_variables = 'cs'
expression = '(0.9-cs)^2'
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L eps_sq'
prop_values = '0.7 0.7 1.0 '
[../]
[]
[Kernels]
active = 'PhaseConc ChemPotSolute CHBulk ACBulkF ACBulkC ACInterface dcdt detadt ckernel'
# enforce c = (1-h(eta))*cl + h(eta)*cs
[./PhaseConc]
type = KKSPhaseConcentration
ca = cl
variable = cs
c = c
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotSolute]
type = KKSPhaseChemicalPotential
variable = cl
cb = cs
fa_name = fl
fb_name = fs
[../]
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = KKSSplitCHCRes
variable = c
ca = cl
fa_name = fl
w = w
[../]
[./dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[./ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = fl
fb_name = fs
w = 1.0
coupled_variables = 'cl cs'
[../]
[./ACBulkC]
type = KKSACBulkC
variable = eta
ca = cl
cb = cs
fa_name = fl
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = eps_sq
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[AuxKernels]
[./GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fl
fb_name = fs
w = 1.0
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'asm ilu nonzero'
l_max_its = 100
nl_max_its = 100
num_steps = 50
dt = 0.1
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[VectorPostprocessors]
[./c]
type = LineValueSampler
start_point = '-25 0 0'
end_point = '25 0 0'
variable = c
num_points = 151
sort_by = id
execute_on = timestep_end
[../]
[./eta]
type = LineValueSampler
start_point = '-25 0 0'
end_point = '25 0 0'
variable = eta
num_points = 151
sort_by = id
execute_on = timestep_end
[../]
[]
[Outputs]
exodus = true
[./csv]
type = CSV
execute_on = final
[../]
[]
(modules/combined/examples/publications/rapid_dev/fig8.i)
#
# Fig. 8 input for 10.1016/j.commatsci.2017.02.017
# D. Schwen et al./Computational Materials Science 132 (2017) 36-45
# Two growing particles with differnet anisotropic Eigenstrains
#
[Mesh]
[./gen]
type = GeneratedMeshGenerator
dim = 2
nx = 80
ny = 40
xmin = -20
xmax = 20
ymin = 0
ymax = 20
elem_type = QUAD4
[../]
[./cnode]
type = ExtraNodesetGenerator
input = gen
coord = '0.0 0.0'
new_boundary = 100
tolerance = 0.1
[../]
[]
[GlobalParams]
# CahnHilliard needs the third derivatives
derivative_order = 3
enable_jit = true
displacements = 'disp_x disp_y'
int_width = 1
[]
# AuxVars to compute the free energy density for outputting
[AuxVariables]
[./local_energy]
order = CONSTANT
family = MONOMIAL
[../]
[./cross_energy]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./local_free_energy]
type = TotalFreeEnergy
variable = local_energy
interfacial_vars = 'c'
kappa_names = 'kappa_c'
additional_free_energy = cross_energy
execute_on = 'INITIAL TIMESTEP_END'
[../]
[./cross_terms]
type = CrossTermGradientFreeEnergy
variable = cross_energy
interfacial_vars = 'eta1 eta2 eta3'
kappa_names = 'kappa11 kappa12 kappa13
kappa21 kappa22 kappa23
kappa31 kappa32 kappa33'
execute_on = 'INITIAL TIMESTEP_END'
[../]
[]
# particle x positions and radius
P1X=8
P2X=-4
PR=2
[Variables]
# Solute concentration variable
[./c]
[./InitialCondition]
type = SpecifiedSmoothCircleIC
x_positions = '${P1X} ${P2X}'
y_positions = '0 0'
z_positions = '0 0'
radii = '${PR} ${PR}'
outvalue = 0.5
invalue = 0.9
[../]
[../]
[./w]
[../]
# Order parameter for the Matrix
[./eta1]
[./InitialCondition]
type = SpecifiedSmoothCircleIC
x_positions = '${P1X} ${P2X}'
y_positions = '0 0'
z_positions = '0 0'
radii = '${PR} ${PR}'
outvalue = 1.0
invalue = 0.0
[../]
[../]
# Order parameters for the 2 different inclusion orientations
[./eta2]
[./InitialCondition]
type = SmoothCircleIC
x1 = ${P2X}
y1 = 0
radius = ${PR}
invalue = 1.0
outvalue = 0.0
[../]
[../]
[./eta3]
[./InitialCondition]
type = SmoothCircleIC
x1 = ${P1X}
y1 = 0
radius = ${PR}
invalue = 1.0
outvalue = 0.0
[../]
[../]
# Lagrange-multiplier
[./lambda]
initial_condition = 1.0
[../]
[]
[Modules]
[./TensorMechanics]
[./Master]
[./all]
add_variables = true
strain = SMALL
eigenstrain_names = eigenstrain
[../]
[../]
[../]
[]
[Kernels]
# Split Cahn-Hilliard kernels
[./c_res]
type = SplitCHParsed
variable = c
f_name = F
args = 'eta1 eta2 eta3'
kappa_name = kappa_c
w = w
[../]
[./wres]
type = SplitCHWRes
variable = w
mob_name = M
[../]
[./time]
type = CoupledTimeDerivative
variable = w
v = c
[../]
# Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 1
[./deta1dt]
type = TimeDerivative
variable = eta1
[../]
[./ACBulk1]
type = AllenCahn
variable = eta1
args = 'eta2 eta3 c'
mob_name = L1
f_name = F
[../]
[./ACInterface1]
type = ACMultiInterface
variable = eta1
etas = 'eta1 eta2 eta3'
mob_name = L1
kappa_names = 'kappa11 kappa12 kappa13'
[../]
[./lagrange1]
type = SwitchingFunctionConstraintEta
variable = eta1
h_name = h1
lambda = lambda
[../]
# Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 2
[./deta2dt]
type = TimeDerivative
variable = eta2
[../]
[./ACBulk2]
type = AllenCahn
variable = eta2
args = 'eta1 eta3 c'
mob_name = L2
f_name = F
[../]
[./ACInterface2]
type = ACMultiInterface
variable = eta2
etas = 'eta1 eta2 eta3'
mob_name = L2
kappa_names = 'kappa21 kappa22 kappa23'
[../]
[./lagrange2]
type = SwitchingFunctionConstraintEta
variable = eta2
h_name = h2
lambda = lambda
[../]
# Allen-Cahn and Lagrange-multiplier constraint kernels for order parameter 3
[./deta3dt]
type = TimeDerivative
variable = eta3
[../]
[./ACBulk3]
type = AllenCahn
variable = eta3
args = 'eta1 eta2 c'
mob_name = L3
f_name = F
[../]
[./ACInterface3]
type = ACMultiInterface
variable = eta3
etas = 'eta1 eta2 eta3'
mob_name = L3
kappa_names = 'kappa31 kappa32 kappa33'
[../]
[./lagrange3]
type = SwitchingFunctionConstraintEta
variable = eta3
h_name = h3
lambda = lambda
[../]
# Lagrange-multiplier constraint kernel for lambda
[./lagrange]
type = SwitchingFunctionConstraintLagrange
variable = lambda
etas = 'eta1 eta2 eta3'
h_names = 'h1 h2 h3'
epsilon = 1e-6
[../]
[]
[Materials]
# declare a few constants, such as mobilities (L,M) and interface gradient prefactors (kappa*)
[./consts]
type = GenericConstantMaterial
block = 0
prop_names = 'M kappa_c L1 L2 L3 kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
prop_values = '0.2 0.5 1 1 1 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00 '
[../]
# We use this to output the level of constraint enforcement
# ideally it should be 0 everywhere, if the constraint is fully enforced
[./etasummat]
type = ParsedMaterial
property_name = etasum
coupled_variables = 'eta1 eta2 eta3'
material_property_names = 'h1 h2 h3'
expression = 'h1+h2+h3-1'
outputs = exodus
[../]
# This parsed material creates a single property for visualization purposes.
# It will be 0 for phase 1, -1 for phase 2, and 1 for phase 3
[./phasemap]
type = ParsedMaterial
property_name = phase
coupled_variables = 'eta2 eta3'
expression = 'if(eta3>0.5,1,0)-if(eta2>0.5,1,0)'
outputs = exodus
[../]
# global mechanical properties
[./elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '400 400'
fill_method = symmetric_isotropic
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
# eigenstrain
[./eigenstrain_2]
type = GenericConstantRankTwoTensor
tensor_name = s2
tensor_values = '0 -0.05 0 0 0 0'
[../]
[./eigenstrain_3]
type = GenericConstantRankTwoTensor
tensor_name = s3
tensor_values = '-0.05 0 0 0 0 0'
[../]
[./eigenstrain]
type = CompositeEigenstrain
weights = 'h2 h3'
tensors = 's2 s3'
args = 'eta2 eta3'
eigenstrain_name = eigenstrain
[../]
# switching functions
[./switching1]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
h_order = SIMPLE
[../]
[./switching2]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
h_order = SIMPLE
[../]
[./switching3]
type = SwitchingFunctionMaterial
function_name = h3
eta = eta3
h_order = SIMPLE
[../]
[./barrier]
type = MultiBarrierFunctionMaterial
etas = 'eta1 eta2 eta3'
[../]
# chemical free energies
[./chemical_free_energy_1]
type = DerivativeParsedMaterial
property_name = Fc1
expression = '4*c^2'
coupled_variables = 'c'
derivative_order = 2
[../]
[./chemical_free_energy_2]
type = DerivativeParsedMaterial
property_name = Fc2
expression = '(c-0.9)^2-0.4'
coupled_variables = 'c'
derivative_order = 2
[../]
[./chemical_free_energy_3]
type = DerivativeParsedMaterial
property_name = Fc3
expression = '(c-0.9)^2-0.5'
coupled_variables = 'c'
derivative_order = 2
[../]
# global chemical free energy
[./chemical_free_energy]
type = DerivativeMultiPhaseMaterial
f_name = Fc
fi_names = 'Fc1 Fc2 Fc3'
hi_names = 'h1 h2 h3'
etas = 'eta1 eta2 eta3'
coupled_variables = 'c'
W = 3
[../]
# global elastic free energy
[./elastic_free_energy]
type = ElasticEnergyMaterial
f_name = Fe
args = 'eta2 eta3'
outputs = exodus
output_properties = Fe
derivative_order = 2
[../]
# Penalize phase 2 and 3 coexistence
[./multi_phase_penalty]
type = DerivativeParsedMaterial
property_name = Fp
expression = '50*(eta2*eta3)^2'
coupled_variables = 'eta2 eta3'
derivative_order = 2
outputs = exodus
output_properties = Fp
[../]
# free energy
[./free_energy]
type = DerivativeSumMaterial
property_name = F
sum_materials = 'Fc Fe Fp'
coupled_variables = 'c eta1 eta2 eta3'
derivative_order = 2
[../]
[]
[BCs]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = 100
variable = disp_x
value = 0
[../]
# fix side point x coordinate to inhibit rotation
[./angularfix]
type = DirichletBC
boundary = bottom
variable = disp_y
value = 0
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
# We monitor the total free energy and the total solute concentration (should be constant)
[Postprocessors]
[./total_free_energy]
type = ElementIntegralVariablePostprocessor
variable = local_energy
execute_on = 'INITIAL TIMESTEP_END'
[../]
[./total_solute]
type = ElementIntegralVariablePostprocessor
variable = c
execute_on = 'INITIAL TIMESTEP_END'
[../]
[]
[Executioner]
type = Transient
scheme = bdf2
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type'
petsc_options_value = 'asm lu'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-8
nl_abs_tol = 1.0e-10
start_time = 0.0
end_time = 12.0
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 8
iteration_window = 1
dt = 0.01
[../]
[]
[Outputs]
print_linear_residuals = false
execute_on = 'INITIAL TIMESTEP_END'
exodus = true
[./table]
type = CSV
delimiter = ' '
[../]
[]
[Debug]
# show_var_residual_norms = true
[]
(modules/phase_field/test/tests/KKS_system/kks_example_split.i)
#
# KKS toy problem in the split form
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 15
ny = 15
nz = 0
xmin = -2.5
xmax = 2.5
ymin = -2.5
ymax = 2.5
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[AuxVariables]
[./Fglobal]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# hydrogen concentration
[./c]
order = FIRST
family = LAGRANGE
[../]
# chemical potential
[./w]
order = FIRST
family = LAGRANGE
[../]
# hydrogen phase concentration (matrix)
[./cm]
order = FIRST
family = LAGRANGE
initial_condition = 0.0
[../]
# hydrogen phase concentration (delta phase)
[./cd]
order = FIRST
family = LAGRANGE
initial_condition = 0.0
[../]
[]
[ICs]
[./eta]
variable = eta
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 1.5
invalue = 0.2
outvalue = 0.1
int_width = 0.75
[../]
[./c]
variable = c
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 1.5
invalue = 0.6
outvalue = 0.4
int_width = 0.75
[../]
[]
[BCs]
[./Periodic]
[./all]
variable = 'eta w c cm cd'
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
# Free energy of the matrix
[./fm]
type = DerivativeParsedMaterial
property_name = fm
coupled_variables = 'cm'
expression = '(0.1-cm)^2'
[../]
# Free energy of the delta phase
[./fd]
type = DerivativeParsedMaterial
property_name = fd
coupled_variables = 'cd'
expression = '(0.9-cd)^2'
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L kappa'
prop_values = '0.7 0.7 0.4 '
[../]
[]
[Kernels]
# full transient
active = 'PhaseConc ChemPotVacancies CHBulk ACBulkF ACBulkC ACInterface dcdt detadt ckernel'
# enforce c = (1-h(eta))*cm + h(eta)*cd
[./PhaseConc]
type = KKSPhaseConcentration
ca = cm
variable = cd
c = c
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotVacancies]
type = KKSPhaseChemicalPotential
variable = cm
cb = cd
fa_name = fm
fb_name = fd
[../]
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = KKSSplitCHCRes
variable = c
ca = cm
fa_name = fm
w = w
[../]
[./dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[./ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = fm
fb_name = fd
coupled_variables = 'cm cd'
w = 0.4
[../]
[./ACBulkC]
type = KKSACBulkC
variable = eta
ca = cm
cb = cd
fa_name = fm
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[AuxKernels]
[./GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fm
fb_name = fd
w = 0.4
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pctype -sub_pc_type -sub_pc_factor_shift_type -pc_factor_shift_type'
petsc_options_value = ' asm lu nonzero nonzero'
l_max_its = 100
nl_max_its = 100
num_steps = 3
dt = 0.1
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Outputs]
file_base = kks_example_split
exodus = true
[]
(modules/phase_field/examples/kim-kim-suzuki/kks_example_ternary.i)
#
# KKS ternary (3 chemical component) system example in the split form
# We track c1 and c2 only, since c1 + c2 + c3 = 1
#
[Mesh]
type = GeneratedMesh
dim = 2
nx = 150
ny = 15
nz = 0
xmin = -25
xmax = 25
ymin = -2.5
ymax = 2.5
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[AuxVariables]
[./Fglobal]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# solute 1 concentration
[./c1]
order = FIRST
family = LAGRANGE
[../]
# solute 2 concentration
[./c2]
order = FIRST
family = LAGRANGE
[../]
# chemical potential solute 1
[./w1]
order = FIRST
family = LAGRANGE
[../]
# chemical potential solute 2
[./w2]
order = FIRST
family = LAGRANGE
[../]
# Liquid phase solute 1 concentration
[./c1l]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[../]
# Liquid phase solute 2 concentration
[./c2l]
order = FIRST
family = LAGRANGE
initial_condition = 0.05
[../]
# Solid phase solute 1 concentration
[./c1s]
order = FIRST
family = LAGRANGE
initial_condition = 0.8
[../]
# Solid phase solute 2 concentration
[./c2s]
order = FIRST
family = LAGRANGE
initial_condition = 0.1
[../]
[]
[Functions]
[./ic_func_eta]
type = ParsedFunction
expression = '0.5*(1.0-tanh((x)/sqrt(2.0)))'
[../]
[./ic_func_c1]
type = ParsedFunction
expression = '0.8*(0.5*(1.0-tanh(x/sqrt(2.0))))^3*(6*(0.5*(1.0-tanh(x/sqrt(2.0))))^2-15*(0.5*(1.0-tanh(x/sqrt(2.0))))+10)+0.1*(1-(0.5*(1.0-tanh(x/sqrt(2.0))))^3*(6*(0.5*(1.0-tanh(x/sqrt(2.0))))^2-15*(0.5*(1.0-tanh(x/sqrt(2.0))))+10))'
[../]
[./ic_func_c2]
type = ParsedFunction
expression = '0.1*(0.5*(1.0-tanh(x/sqrt(2.0))))^3*(6*(0.5*(1.0-tanh(x/sqrt(2.0))))^2-15*(0.5*(1.0-tanh(x/sqrt(2.0))))+10)+0.05*(1-(0.5*(1.0-tanh(x/sqrt(2.0))))^3*(6*(0.5*(1.0-tanh(x/sqrt(2.0))))^2-15*(0.5*(1.0-tanh(x/sqrt(2.0))))+10))'
[../]
[]
[ICs]
[./eta]
variable = eta
type = FunctionIC
function = ic_func_eta
[../]
[./c1]
variable = c1
type = FunctionIC
function = ic_func_c1
[../]
[./c2]
variable = c2
type = FunctionIC
function = ic_func_c2
[../]
[]
[Materials]
# Free energy of the liquid
[./fl]
type = DerivativeParsedMaterial
property_name = fl
coupled_variables = 'c1l c2l'
expression = '(0.1-c1l)^2+(0.05-c2l)^2'
[../]
# Free energy of the solid
[./fs]
type = DerivativeParsedMaterial
property_name = fs
coupled_variables = 'c1s c2s'
expression = '(0.8-c1s)^2+(0.1-c2s)^2'
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L eps_sq'
prop_values = '0.7 0.7 1.0 '
[../]
[]
[Kernels]
# enforce c1 = (1-h(eta))*c1l + h(eta)*c1s
[./PhaseConc1]
type = KKSPhaseConcentration
ca = c1l
variable = c1s
c = c1
eta = eta
[../]
# enforce c2 = (1-h(eta))*c2l + h(eta)*c2s
[./PhaseConc2]
type = KKSPhaseConcentration
ca = c2l
variable = c2s
c = c2
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotSolute1]
type = KKSPhaseChemicalPotential
variable = c1l
cb = c1s
fa_name = fl
fb_name = fs
args_a = 'c2l'
args_b = 'c2s'
[../]
[./ChemPotSolute2]
type = KKSPhaseChemicalPotential
variable = c2l
cb = c2s
fa_name = fl
fb_name = fs
args_a = 'c1l'
args_b = 'c1s'
[../]
#
# Cahn-Hilliard Equations
#
[./CHBulk1]
type = KKSSplitCHCRes
variable = c1
ca = c1l
fa_name = fl
w = w1
args_a = 'c2l'
[../]
[./CHBulk2]
type = KKSSplitCHCRes
variable = c2
ca = c2l
fa_name = fl
w = w2
args_a = 'c1l'
[../]
[./dc1dt]
type = CoupledTimeDerivative
variable = w1
v = c1
[../]
[./dc2dt]
type = CoupledTimeDerivative
variable = w2
v = c2
[../]
[./w1kernel]
type = SplitCHWRes
mob_name = M
variable = w1
[../]
[./w2kernel]
type = SplitCHWRes
mob_name = M
variable = w2
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = fl
fb_name = fs
w = 1.0
coupled_variables = 'c1l c1s c2l c2s'
[../]
[./ACBulkC1]
type = KKSACBulkC
variable = eta
ca = c1l
cb = c1s
fa_name = fl
coupled_variables = 'c2l'
[../]
[./ACBulkC2]
type = KKSACBulkC
variable = eta
ca = c2l
cb = c2s
fa_name = fl
coupled_variables = 'c1l'
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = eps_sq
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[AuxKernels]
[./GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fl
fb_name = fs
w = 1.0
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'asm ilu nonzero'
l_max_its = 100
nl_max_its = 100
num_steps = 50
dt = 0.1
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Outputs]
exodus = true
[]
(modules/combined/examples/phase_field-mechanics/kks_mechanics_KHS.i)
# KKS phase-field model coupled with elasticity using Khachaturyan's scheme as
# described in L.K. Aagesen et al., Computational Materials Science, 140, 10-21 (2017)
# Original run #170403a
[Mesh]
type = GeneratedMesh
dim = 3
nx = 640
ny = 1
nz = 1
xmin = -10
xmax = 10
ymin = 0
ymax = 0.03125
zmin = 0
zmax = 0.03125
elem_type = HEX8
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# solute concentration
[./c]
order = FIRST
family = LAGRANGE
[../]
# chemical potential
[./w]
order = FIRST
family = LAGRANGE
[../]
# solute phase concentration (matrix)
[./cm]
order = FIRST
family = LAGRANGE
[../]
# solute phase concentration (precipitate)
[./cp]
order = FIRST
family = LAGRANGE
[../]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[./disp_z]
order = FIRST
family = LAGRANGE
[../]
[]
[ICs]
[./eta_ic]
variable = eta
type = FunctionIC
function = ic_func_eta
block = 0
[../]
[./c_ic]
variable = c
type = FunctionIC
function = ic_func_c
block = 0
[../]
[./w_ic]
variable = w
type = ConstantIC
value = 0.00991
block = 0
[../]
[./cm_ic]
variable = cm
type = ConstantIC
value = 0.131
block = 0
[../]
[./cp_ic]
variable = cp
type = ConstantIC
value = 0.236
block = 0
[../]
[]
[Functions]
[./ic_func_eta]
type = ParsedFunction
expression = '0.5*(1.0+tanh((x)/delta_eta/sqrt(2.0)))'
symbol_names = 'delta_eta'
symbol_values = '0.8034'
[../]
[./ic_func_c]
type = ParsedFunction
expression = '0.2389*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^3*(6*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^2-15*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))+10)+0.1339*(1-(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^3*(6*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))^2-15*(0.5*(1.0+tanh(x/delta/sqrt(2.0))))+10))'
symbol_names = 'delta'
symbol_values = '0.8034'
[../]
[./psi_eq_int]
type = ParsedFunction
expression = 'volume*psi_alpha'
symbol_names = 'volume psi_alpha'
symbol_values = 'volume psi_alpha'
[../]
[./gamma]
type = ParsedFunction
expression = '(psi_int - psi_eq_int) / dy / dz'
symbol_names = 'psi_int psi_eq_int dy dz'
symbol_values = 'psi_int psi_eq_int 0.03125 0.03125'
[../]
[]
[AuxVariables]
[./sigma11]
order = CONSTANT
family = MONOMIAL
[../]
[./sigma22]
order = CONSTANT
family = MONOMIAL
[../]
[./sigma33]
order = CONSTANT
family = MONOMIAL
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
[../]
[./e33]
order = CONSTANT
family = MONOMIAL
[../]
[./e_el11]
order = CONSTANT
family = MONOMIAL
[../]
[./e_el12]
order = CONSTANT
family = MONOMIAL
[../]
[./e_el22]
order = CONSTANT
family = MONOMIAL
[../]
[./f_el]
order = CONSTANT
family = MONOMIAL
[../]
[./eigen_strain00]
order = CONSTANT
family = MONOMIAL
[../]
[./Fglobal]
order = CONSTANT
family = MONOMIAL
[../]
[./psi]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./matl_sigma11]
type = RankTwoAux
rank_two_tensor = stress
index_i = 0
index_j = 0
variable = sigma11
[../]
[./matl_sigma22]
type = RankTwoAux
rank_two_tensor = stress
index_i = 1
index_j = 1
variable = sigma22
[../]
[./matl_sigma33]
type = RankTwoAux
rank_two_tensor = stress
index_i = 2
index_j = 2
variable = sigma33
[../]
[./matl_e11]
type = RankTwoAux
rank_two_tensor = total_strain
index_i = 0
index_j = 0
variable = e11
[../]
[./f_el]
type = MaterialRealAux
variable = f_el
property = f_el_mat
execute_on = timestep_end
[../]
[./GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fm
fb_name = fp
w = 0.0264
kappa_names = kappa
interfacial_vars = eta
[../]
[./psi_potential]
variable = psi
type = ParsedAux
coupled_variables = 'Fglobal w c f_el sigma11 e11'
expression = 'Fglobal - w*c + f_el - sigma11*e11'
[../]
[]
[BCs]
[./left_x]
type = DirichletBC
variable = disp_x
boundary = left
value = 0
[../]
[./right_x]
type = DirichletBC
variable = disp_x
boundary = right
value = 0
[../]
[./front_y]
type = DirichletBC
variable = disp_y
boundary = front
value = 0
[../]
[./back_y]
type = DirichletBC
variable = disp_y
boundary = back
value = 0
[../]
[./top_z]
type = DirichletBC
variable = disp_z
boundary = top
value = 0
[../]
[./bottom_z]
type = DirichletBC
variable = disp_z
boundary = bottom
value = 0
[../]
[]
[Materials]
# Chemical free energy of the matrix
[./fm]
type = DerivativeParsedMaterial
property_name = fm
coupled_variables = 'cm'
expression = '6.55*(cm-0.13)^2'
[../]
# Chemical Free energy of the precipitate phase
[./fp]
type = DerivativeParsedMaterial
property_name = fp
coupled_variables = 'cp'
expression = '6.55*(cp-0.235)^2'
[../]
# Elastic energy of the precipitate
[./elastic_free_energy_p]
type = ElasticEnergyMaterial
f_name = f_el_mat
args = 'eta'
outputs = exodus
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# 1- h(eta), putting in function explicitly
[./one_minus_h_eta_explicit]
type = DerivativeParsedMaterial
property_name = one_minus_h_explicit
coupled_variables = eta
expression = 1-eta^3*(6*eta^2-15*eta+10)
outputs = exodus
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L kappa misfit'
prop_values = '0.7 0.7 0.01704 0.00377'
[../]
#Mechanical properties
[./Stiffness_matrix]
type = ComputeElasticityTensor
base_name = C_matrix
C_ijkl = '103.3 74.25 74.25 103.3 74.25 103.3 46.75 46.75 46.75'
fill_method = symmetric9
[../]
[./Stiffness_ppt]
type = ComputeElasticityTensor
C_ijkl = '100.7 71.45 71.45 100.7 71.45 100.7 50.10 50.10 50.10'
base_name = C_ppt
fill_method = symmetric9
[../]
[./C]
type = CompositeElasticityTensor
args = eta
tensors = 'C_matrix C_ppt'
weights = 'one_minus_h_explicit h'
[../]
[./stress]
type = ComputeLinearElasticStress
[../]
[./strain]
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = 'eigenstrain_ppt'
[../]
[./eigen_strain]
type = ComputeVariableEigenstrain
eigen_base = '0.00377 0.00377 0.00377 0 0 0'
prefactor = h
args = eta
eigenstrain_name = 'eigenstrain_ppt'
[../]
[]
[Kernels]
[./TensorMechanics]
displacements = 'disp_x disp_y disp_z'
[../]
# enforce c = (1-h(eta))*cm + h(eta)*cp
[./PhaseConc]
type = KKSPhaseConcentration
ca = cm
variable = cp
c = c
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotVacancies]
type = KKSPhaseChemicalPotential
variable = cm
cb = cp
fa_name = fm
fb_name = fp
[../]
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = KKSSplitCHCRes
variable = c
ca = cm
fa_name = fm
w = w
[../]
[./dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[./ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = fm
fb_name = fp
w = 0.0264
args = 'cp cm'
[../]
[./ACBulkC]
type = KKSACBulkC
variable = eta
ca = cm
cb = cp
fa_name = fm
[../]
[./ACBulk_el] #This adds df_el/deta for strain interpolation
type = AllenCahn
variable = eta
f_name = f_el_mat
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'asm ilu nonzero'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-8
nl_abs_tol = 1.0e-11
num_steps = 200
[./TimeStepper]
type = SolutionTimeAdaptiveDT
dt = 0.5
[../]
[]
[Postprocessors]
[./f_el_int]
type = ElementIntegralMaterialProperty
mat_prop = f_el_mat
[../]
[./c_alpha]
type = SideAverageValue
boundary = left
variable = c
[../]
[./c_beta]
type = SideAverageValue
boundary = right
variable = c
[../]
[./e11_alpha]
type = SideAverageValue
boundary = left
variable = e11
[../]
[./e11_beta]
type = SideAverageValue
boundary = right
variable = e11
[../]
[./s11_alpha]
type = SideAverageValue
boundary = left
variable = sigma11
[../]
[./s22_alpha]
type = SideAverageValue
boundary = left
variable = sigma22
[../]
[./s33_alpha]
type = SideAverageValue
boundary = left
variable = sigma33
[../]
[./s11_beta]
type = SideAverageValue
boundary = right
variable = sigma11
[../]
[./s22_beta]
type = SideAverageValue
boundary = right
variable = sigma22
[../]
[./s33_beta]
type = SideAverageValue
boundary = right
variable = sigma33
[../]
[./f_el_alpha]
type = SideAverageValue
boundary = left
variable = f_el
[../]
[./f_el_beta]
type = SideAverageValue
boundary = right
variable = f_el
[../]
[./f_c_alpha]
type = SideAverageValue
boundary = left
variable = Fglobal
[../]
[./f_c_beta]
type = SideAverageValue
boundary = right
variable = Fglobal
[../]
[./chem_pot_alpha]
type = SideAverageValue
boundary = left
variable = w
[../]
[./chem_pot_beta]
type = SideAverageValue
boundary = right
variable = w
[../]
[./psi_alpha]
type = SideAverageValue
boundary = left
variable = psi
[../]
[./psi_beta]
type = SideAverageValue
boundary = right
variable = psi
[../]
[./total_energy]
type = ElementIntegralVariablePostprocessor
variable = Fglobal
[../]
# Get simulation cell size from postprocessor
[./volume]
type = ElementIntegralMaterialProperty
mat_prop = 1
[../]
[./psi_eq_int]
type = FunctionValuePostprocessor
function = psi_eq_int
[../]
[./psi_int]
type = ElementIntegralVariablePostprocessor
variable = psi
[../]
[./gamma]
type = FunctionValuePostprocessor
function = gamma
[../]
[./int_position]
type = FindValueOnLine
start_point = '-10 0 0'
end_point = '10 0 0'
v = eta
target = 0.5
[../]
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Outputs]
[./exodus]
type = Exodus
time_step_interval = 20
[../]
checkpoint = true
[./csv]
type = CSV
execute_on = 'final'
[../]
[]
(modules/phase_field/test/tests/KKS_system/kks_example_offset.i)
#
# KKS toy problem in the split form
# This has an offset in the minima of the free energies so there will be a shift
# in equilibrium composition
[Mesh]
type = GeneratedMesh
dim = 2
nx = 15
ny = 15
nz = 0
xmin = -2.5
xmax = 2.5
ymin = -2.5
ymax = 2.5
zmin = 0
zmax = 0
elem_type = QUAD4
[]
[AuxVariables]
[./Fglobal]
order = CONSTANT
family = MONOMIAL
[../]
[]
[Variables]
# order parameter
[./eta]
order = FIRST
family = LAGRANGE
[../]
# hydrogen concentration
[./c]
order = FIRST
family = LAGRANGE
[../]
# chemical potential
[./w]
order = FIRST
family = LAGRANGE
[../]
# hydrogen phase concentration (matrix)
[./cm]
order = FIRST
family = LAGRANGE
initial_condition = 0.0
[../]
# hydrogen phase concentration (delta phase)
[./cd]
order = FIRST
family = LAGRANGE
initial_condition = 0.0
[../]
[]
[ICs]
[./eta]
variable = eta
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 1.5
invalue = 0.2
outvalue = 0.1
int_width = 0.75
[../]
[./c]
variable = c
type = SmoothCircleIC
x1 = 0.0
y1 = 0.0
radius = 1.5
invalue = 0.6
outvalue = 0.4
int_width = 0.75
[../]
[]
[BCs]
[./Periodic]
[./all]
variable = 'eta w c cm cd'
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
# Free energy of the matrix
[./fm]
type = DerivativeParsedMaterial
property_name = fm
coupled_variables = 'cm'
expression = '(0.1-cm)^2'
[../]
# Free energy of the delta phase
[./fd]
type = DerivativeParsedMaterial
property_name = fd
coupled_variables = 'cd'
expression = '(0.9-cd)^2+0.5'
[../]
# h(eta)
[./h_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta
[../]
# g(eta)
[./g_eta]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta
[../]
# constant properties
[./constants]
type = GenericConstantMaterial
prop_names = 'M L kappa'
prop_values = '0.7 0.7 0.4 '
[../]
[]
[Kernels]
# full transient
active = 'PhaseConc ChemPotVacancies CHBulk ACBulkF ACBulkC ACInterface dcdt detadt ckernel'
# enforce c = (1-h(eta))*cm + h(eta)*cd
[./PhaseConc]
type = KKSPhaseConcentration
ca = cm
variable = cd
c = c
eta = eta
[../]
# enforce pointwise equality of chemical potentials
[./ChemPotVacancies]
type = KKSPhaseChemicalPotential
variable = cm
cb = cd
fa_name = fm
fb_name = fd
[../]
#
# Cahn-Hilliard Equation
#
[./CHBulk]
type = KKSSplitCHCRes
variable = c
ca = cm
fa_name = fm
w = w
[../]
[./dcdt]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[./ckernel]
type = SplitCHWRes
mob_name = M
variable = w
[../]
#
# Allen-Cahn Equation
#
[./ACBulkF]
type = KKSACBulkF
variable = eta
fa_name = fm
fb_name = fd
coupled_variables = 'cm cd'
w = 0.4
[../]
[./ACBulkC]
type = KKSACBulkC
variable = eta
ca = cm
cb = cd
fa_name = fm
[../]
[./ACInterface]
type = ACInterface
variable = eta
kappa_name = kappa
[../]
[./detadt]
type = TimeDerivative
variable = eta
[../]
[]
[AuxKernels]
[./GlobalFreeEnergy]
variable = Fglobal
type = KKSGlobalFreeEnergy
fa_name = fm
fb_name = fd
w = 0.4
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pctype -sub_pc_type -sub_pc_factor_shift_type -pc_factor_shift_type'
petsc_options_value = ' asm lu nonzero nonzero'
l_max_its = 100
nl_max_its = 100
num_steps = 3
dt = 0.1
[]
#
# Precondition using handcoded off-diagonal terms
#
[Preconditioning]
[./full]
type = SMP
full = true
[../]
[]
[Outputs]
file_base = kks_example_offset
exodus = true
[]
(modules/phase_field/test/tests/MultiPhase/penalty.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 14
ny = 10
nz = 0
xmin = 10
xmax = 40
ymin = 15
ymax = 35
elem_type = QUAD4
[]
[GlobalParams]
penalty = 5
[]
[Variables]
[./c]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 25.0
y1 = 25.0
radius = 6.0
invalue = 0.9
outvalue = 0.1
int_width = 3.0
[../]
[../]
[./w]
order = FIRST
family = LAGRANGE
[../]
[./eta1]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = SmoothCircleIC
x1 = 30.0
y1 = 25.0
radius = 4.0
invalue = 0.9
outvalue = 0.1
int_width = 2.0
[../]
[../]
[./eta2]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[../]
[]
[Kernels]
[./deta1dt]
type = TimeDerivative
variable = eta1
[../]
[./ACBulk1]
type = AllenCahn
variable = eta1
coupled_variables = 'c eta2'
f_name = F
[../]
[./ACInterface1]
type = ACInterface
variable = eta1
kappa_name = kappa_eta
[../]
[./penalty1]
type = SwitchingFunctionPenalty
variable = eta1
etas = 'eta1 eta2'
h_names = 'h1 h2'
[../]
[./deta2dt]
type = TimeDerivative
variable = eta2
[../]
[./ACBulk2]
type = AllenCahn
variable = eta2
coupled_variables = 'c eta1'
f_name = F
[../]
[./ACInterface2]
type = ACInterface
variable = eta2
kappa_name = kappa_eta
[../]
[./penalty2]
type = SwitchingFunctionPenalty
variable = eta2
etas = 'eta1 eta2'
h_names = 'h1 h2'
[../]
[./c_res]
type = SplitCHParsed
variable = c
f_name = F
kappa_name = kappa_c
w = w
coupled_variables = 'eta1 eta2'
[../]
[./w_res]
type = SplitCHWRes
variable = w
mob_name = M
[../]
[./time1]
type = CoupledTimeDerivative
variable = w
v = c
[../]
[]
[BCs]
[./Periodic]
[./All]
auto_direction = 'x y'
[../]
[../]
[]
[Materials]
[./consts]
type = GenericConstantMaterial
prop_names = 'L kappa_eta'
prop_values = '1 1 '
[../]
[./consts2]
type = GenericConstantMaterial
prop_names = 'M kappa_c'
prop_values = '1 1'
[../]
[./hsum]
type = ParsedMaterial
expression = h1+h2
property_name = hsum
material_property_names = 'h1 h2'
coupled_variables = 'c'
outputs = exodus
[../]
[./switching1]
type = SwitchingFunctionMaterial
function_name = h1
eta = eta1
h_order = SIMPLE
[../]
[./switching2]
type = SwitchingFunctionMaterial
function_name = h2
eta = eta2
h_order = SIMPLE
[../]
[./barrier]
type = MultiBarrierFunctionMaterial
etas = 'eta1 eta2'
[../]
[./free_energy_A]
type = DerivativeParsedMaterial
property_name = Fa
coupled_variables = 'c'
expression = '(c-0.1)^2'
derivative_order = 2
[../]
[./free_energy_B]
type = DerivativeParsedMaterial
property_name = Fb
coupled_variables = 'c'
expression = '(c-0.9)^2'
derivative_order = 2
[../]
[./free_energy]
type = DerivativeMultiPhaseMaterial
property_name = F
fi_names = 'Fa Fb'
hi_names = 'h1 h2'
etas = 'eta1 eta2'
coupled_variables = 'c'
derivative_order = 2
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
scheme = 'bdf2'
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -sub_pc_type'
petsc_options_value = 'asm lu'
l_max_its = 15
l_tol = 1.0e-6
nl_max_its = 50
nl_rel_tol = 1.0e-7
nl_abs_tol = 1.0e-9
start_time = 0.0
num_steps = 2
dt = 0.05
dtmin = 0.01
[]
[Debug]
# show_var_residual_norms = true
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(modules/combined/test/tests/linear_elasticity/extra_stress.i)
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
type = GeneratedMesh
dim = 2
nx = 128
ny = 1
xmax = 3.2
ymax = 0.025
elem_type = QUAD4
[]
[Modules/TensorMechanics/Master/All]
add_variables = true
generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_xy strain_yy'
[]
[AuxVariables]
[./c]
[../]
[]
[ICs]
[./c_IC]
type = BoundingBoxIC
variable = c
x1 = -1
y1 = -1
x2 = 1.6
y2 = 1
inside = 0
outside = 1
block = 0
[../]
[]
[Materials]
[./elasticity_tensor]
type = ComputeElasticityTensor
block = 0
C_ijkl = '104 74 74 104 74 104 47.65 47.65 47.65'
fill_method = symmetric9
base_name = matrix
[../]
[./stress]
type = ComputeLinearElasticStress
block = 0
base_name = matrix
[../]
[./strain]
type = ComputeSmallStrain
block = 0
base_name = matrix
[../]
[./elasticity_tensor_ppt]
type = ComputeElasticityTensor
block = 0
C_ijkl = '0.104 0.074 0.074 0.104 0.074 0.104 0.04765 0.04765 0.04765'
fill_method = symmetric9
base_name = ppt
[../]
[./stress_ppt]
type = ComputeLinearElasticStress
block = 0
base_name = ppt
[../]
[./strain_ppt]
type = ComputeSmallStrain
block = 0
base_name = ppt
[../]
[./const_stress]
type = ComputeExtraStressConstant
block = 0
base_name = ppt
extra_stress_tensor = '-0.288 -0.373 -0.2747 0 0 0'
[../]
[./global_stress]
type = TwoPhaseStressMaterial
base_A = matrix
base_B = ppt
[../]
[./switching]
type = SwitchingFunctionMaterial
eta = c
[../]
[]
[BCs]
active = 'left_x right_x bottom_y top_y'
[./bottom_y]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0
[../]
[./left_x]
type = DirichletBC
variable = disp_x
boundary = left
value = 0
[../]
[./right_x]
type = DirichletBC
variable = disp_x
boundary = right
value = 0
[../]
[./top_y]
type = DirichletBC
variable = disp_y
boundary = top
value = 0
[../]
[]
[Preconditioning]
[./SMP]
type = SMP
full = true
[../]
[]
[Executioner]
type = Steady
solve_type = 'NEWTON'
[]
[Outputs]
exodus = true
[]
(modules/phase_field/test/tests/slkks/full_solve.i)
#
# SLKKS two phase example for the BCC and SIGMA phases. The sigma phase contains
# multiple sublattices. Free energy from
# Jacob, Aurelie, Erwin Povoden-Karadeniz, and Ernst Kozeschnik. "Revised thermodynamic
# description of the Fe-Cr system based on an improved sublattice model of the sigma phase."
# Calphad 60 (2018): 16-28.
#
# In this simulation we consider diffusion (Cahn-Hilliard) and phase transformation.
#
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 1
nx = 30
ny = 1
xmin = -25
xmax = 25
[]
[]
[AuxVariables]
[Fglobal]
order = CONSTANT
family = MONOMIAL
[]
[]
[Variables]
# order parameters
[eta1]
initial_condition = 0.5
[]
[eta2]
initial_condition = 0.5
[]
# solute concentration
[cCr]
order = FIRST
family = LAGRANGE
[InitialCondition]
type = FunctionIC
function = '(x+25)/50*0.5+0.1'
[]
[]
# sublattice concentrations (good guesses are needed here! - they can be obtained
# form a static solve like in sublattice_concentrations.i)
[BCC_CR]
[InitialCondition]
type = FunctionIC
function = '(x+25)/50*0.5+0.1'
[]
[]
[SIGMA_0CR]
[InitialCondition]
type = FunctionIC
function = '(x+25)/50*0.17+0.01'
[]
[]
[SIGMA_1CR]
[InitialCondition]
type = FunctionIC
function = '(x+25)/50*0.36+0.02'
[]
[]
[SIGMA_2CR]
[InitialCondition]
type = FunctionIC
function = '(x+25)/50*0.33+0.20'
[]
[]
# Lagrange multiplier
[lambda]
[]
[]
[Materials]
# CALPHAD free energies
[F_BCC_A2]
type = DerivativeParsedMaterial
property_name = F_BCC_A2
outputs = exodus
output_properties = F_BCC_A2
expression = 'BCC_FE:=1-BCC_CR; G := 8.3145*T*(1.0*if(BCC_CR > 1.0e-15,BCC_CR*log(BCC_CR),0) + '
'1.0*if(BCC_FE > 1.0e-15,BCC_FE*plog(BCC_FE,eps),0) + 3.0*if(BCC_VA > '
'1.0e-15,BCC_VA*log(BCC_VA),0))/(BCC_CR + BCC_FE) + 8.3145*T*if(T < '
'548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - 932.5*BCC_CR*BCC_FE*BCC_VA + '
'311.5*BCC_CR*BCC_VA - '
'1043.0*BCC_FE*BCC_VA,-8.13674105561218e-49*T^15/(0.525599232981783*BCC_CR*BCC_FE*BCC_'
'VA*(BCC_CR - BCC_FE) - 0.894055608820709*BCC_CR*BCC_FE*BCC_VA + '
'0.298657718120805*BCC_CR*BCC_VA - BCC_FE*BCC_VA + 9.58772770853308e-13)^15 - '
'4.65558036243985e-30*T^9/(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^9 - '
'1.3485349181899e-10*T^3/(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^3 + 1 - '
'0.905299382744392*(548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'932.5*BCC_CR*BCC_FE*BCC_VA + 311.5*BCC_CR*BCC_VA - 1043.0*BCC_FE*BCC_VA + '
'1.0e-9)/T,if(T < -548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'932.5*BCC_CR*BCC_FE*BCC_VA - 311.5*BCC_CR*BCC_VA + '
'1043.0*BCC_FE*BCC_VA,-8.13674105561218e-49*T^15/(-0.525599232981783*BCC_CR*BCC_FE*BCC'
'_VA*(BCC_CR - BCC_FE) + 0.894055608820709*BCC_CR*BCC_FE*BCC_VA - '
'0.298657718120805*BCC_CR*BCC_VA + BCC_FE*BCC_VA + 9.58772770853308e-13)^15 - '
'4.65558036243985e-30*T^9/(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) '
'+ 0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^9 - '
'1.3485349181899e-10*T^3/(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^3 + 1 - '
'0.905299382744392*(-548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'932.5*BCC_CR*BCC_FE*BCC_VA - 311.5*BCC_CR*BCC_VA + 1043.0*BCC_FE*BCC_VA + '
'1.0e-9)/T,if(T > -548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'932.5*BCC_CR*BCC_FE*BCC_VA - 311.5*BCC_CR*BCC_VA + 1043.0*BCC_FE*BCC_VA & '
'548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - 932.5*BCC_CR*BCC_FE*BCC_VA + '
'311.5*BCC_CR*BCC_VA - 1043.0*BCC_FE*BCC_VA < '
'0,-79209031311018.7*(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^5/T^5 - '
'3.83095660520737e+42*(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^15/T^15 - '
'1.22565886734485e+72*(-0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) + '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA - 0.298657718120805*BCC_CR*BCC_VA + '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^25/T^25,if(T > '
'548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - 932.5*BCC_CR*BCC_FE*BCC_VA + '
'311.5*BCC_CR*BCC_VA - 1043.0*BCC_FE*BCC_VA & 548.2*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - '
'BCC_FE) - 932.5*BCC_CR*BCC_FE*BCC_VA + 311.5*BCC_CR*BCC_VA - 1043.0*BCC_FE*BCC_VA > '
'0,-79209031311018.7*(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^5/T^5 - '
'3.83095660520737e+42*(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^15/T^15 - '
'1.22565886734485e+72*(0.525599232981783*BCC_CR*BCC_FE*BCC_VA*(BCC_CR - BCC_FE) - '
'0.894055608820709*BCC_CR*BCC_FE*BCC_VA + 0.298657718120805*BCC_CR*BCC_VA - '
'BCC_FE*BCC_VA + 9.58772770853308e-13)^25/T^25,0))))*log((2.15*BCC_CR*BCC_FE*BCC_VA - '
'0.008*BCC_CR*BCC_VA + 2.22*BCC_FE*BCC_VA)*if(2.15*BCC_CR*BCC_FE*BCC_VA - '
'0.008*BCC_CR*BCC_VA + 2.22*BCC_FE*BCC_VA <= 0,-1.0,1.0) + 1)/(BCC_CR + BCC_FE) + '
'1.0*(BCC_CR*BCC_VA*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + '
'BCC_FE*BCC_VA*if(T >= 298.15 & T < 1811.0,77358.5*1/T - 23.5143*T*log(T) + 124.134*T '
'- 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= 1811.0 & T < '
'6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - 25383.581,0)))/(BCC_CR '
'+ BCC_FE) + 1.0*(BCC_CR*BCC_FE*BCC_VA*(500.0 - 1.5*T)*(BCC_CR - BCC_FE) + '
'BCC_CR*BCC_FE*BCC_VA*(24600.0 - 14.98*T) + BCC_CR*BCC_FE*BCC_VA*(9.15*T - '
'14000.0)*(BCC_CR - BCC_FE)^2)/(BCC_CR + BCC_FE); G/100000'
coupled_variables = 'BCC_CR'
constant_names = 'BCC_VA T eps'
constant_expressions = '1 1000 0.01'
[]
[F_SIGMA]
type = DerivativeParsedMaterial
property_name = F_SIGMA
outputs = exodus
output_properties = F_SIGMA
expression = 'SIGMA_0FE := 1-SIGMA_0CR; SIGMA_1FE := 1-SIGMA_1CR; SIGMA_2FE := 1-SIGMA_2CR; G := '
'8.3145*T*(10.0*if(SIGMA_0CR > 1.0e-15,SIGMA_0CR*plog(SIGMA_0CR,eps),0) + '
'10.0*if(SIGMA_0FE > 1.0e-15,SIGMA_0FE*plog(SIGMA_0FE,eps),0) + 4.0*if(SIGMA_1CR > '
'1.0e-15,SIGMA_1CR*plog(SIGMA_1CR,eps),0) + 4.0*if(SIGMA_1FE > '
'1.0e-15,SIGMA_1FE*plog(SIGMA_1FE,eps),0) + 16.0*if(SIGMA_2CR > '
'1.0e-15,SIGMA_2CR*plog(SIGMA_2CR,eps),0) + 16.0*if(SIGMA_2FE > '
'1.0e-15,SIGMA_2FE*plog(SIGMA_2FE,eps),0))/(10.0*SIGMA_0CR + 10.0*SIGMA_0FE + '
'4.0*SIGMA_1CR + 4.0*SIGMA_1FE + 16.0*SIGMA_2CR + 16.0*SIGMA_2FE) + '
'(SIGMA_0FE*SIGMA_1CR*SIGMA_2CR*SIGMA_2FE*(-70.0*T - 170400.0) + '
'SIGMA_0FE*SIGMA_1FE*SIGMA_2CR*SIGMA_2FE*(-10.0*T - 330839.0))/(10.0*SIGMA_0CR + '
'10.0*SIGMA_0FE + 4.0*SIGMA_1CR + 4.0*SIGMA_1FE + 16.0*SIGMA_2CR + 16.0*SIGMA_2FE) + '
'(SIGMA_0CR*SIGMA_1CR*SIGMA_2CR*(30.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - '
'26.908*T*log(T) + 157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= '
'2180.0 & T < 6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) '
'+ 132000.0) + SIGMA_0CR*SIGMA_1CR*SIGMA_2FE*(-110.0*T + 16.0*if(T >= 298.15 & T < '
'1811.0,77358.5*1/T - 23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - '
'5.89269e-8*T^3.0 + 1225.7,if(T >= 1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - '
'46.0*T*log(T) + 299.31255*T - 25383.581,0)) + 14.0*if(T >= 298.15 & T < '
'2180.0,139250.0*1/T - 26.908*T*log(T) + 157.48*T + 0.00189435*T^2.0 - '
'1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < 6000.0,-2.88526e+32*T^(-9.0) - '
'50.0*T*log(T) + 344.18*T - 34869.344,0)) + 123500.0) + '
'SIGMA_0CR*SIGMA_1FE*SIGMA_2CR*(4.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 26.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 140486.0) '
'+ SIGMA_0CR*SIGMA_1FE*SIGMA_2FE*(20.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 10.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 148800.0) '
'+ SIGMA_0FE*SIGMA_1CR*SIGMA_2CR*(10.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 20.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 56200.0) + '
'SIGMA_0FE*SIGMA_1CR*SIGMA_2FE*(26.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 4.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 152700.0) '
'+ SIGMA_0FE*SIGMA_1FE*SIGMA_2CR*(14.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 16.0*if(T >= 298.15 & T < 2180.0,139250.0*1/T - 26.908*T*log(T) + '
'157.48*T + 0.00189435*T^2.0 - 1.47721e-6*T^3.0 - 8856.94,if(T >= 2180.0 & T < '
'6000.0,-2.88526e+32*T^(-9.0) - 50.0*T*log(T) + 344.18*T - 34869.344,0)) + 46200.0) + '
'SIGMA_0FE*SIGMA_1FE*SIGMA_2FE*(30.0*if(T >= 298.15 & T < 1811.0,77358.5*1/T - '
'23.5143*T*log(T) + 124.134*T - 0.00439752*T^2.0 - 5.89269e-8*T^3.0 + 1225.7,if(T >= '
'1811.0 & T < 6000.0,2.2960305e+31*T^(-9.0) - 46.0*T*log(T) + 299.31255*T - '
'25383.581,0)) + 173333.0))/(10.0*SIGMA_0CR + 10.0*SIGMA_0FE + 4.0*SIGMA_1CR + '
'4.0*SIGMA_1FE + 16.0*SIGMA_2CR + 16.0*SIGMA_2FE); G/100000'
coupled_variables = 'SIGMA_0CR SIGMA_1CR SIGMA_2CR'
constant_names = 'T eps'
constant_expressions = '1000 0.01'
[]
# h(eta)
[h1]
type = SwitchingFunctionMaterial
function_name = h1
h_order = HIGH
eta = eta1
[]
[h2]
type = SwitchingFunctionMaterial
function_name = h2
h_order = HIGH
eta = eta2
[]
# g(eta)
[g1]
type = BarrierFunctionMaterial
function_name = g1
g_order = SIMPLE
eta = eta1
[]
[g2]
type = BarrierFunctionMaterial
function_name = g2
g_order = SIMPLE
eta = eta2
[]
# constant properties
[constants]
type = GenericConstantMaterial
prop_names = 'D L kappa'
prop_values = '10 1 0.1 '
[]
# Coefficients for diffusion equation
[Dh1]
type = DerivativeParsedMaterial
material_property_names = 'D h1(eta1)'
expression = D*h1
property_name = Dh1
coupled_variables = eta1
derivative_order = 1
[]
[Dh2a]
type = DerivativeParsedMaterial
material_property_names = 'D h2(eta2)'
expression = D*h2*10/30
property_name = Dh2a
coupled_variables = eta2
derivative_order = 1
[]
[Dh2b]
type = DerivativeParsedMaterial
material_property_names = 'D h2(eta2)'
expression = D*h2*4/30
property_name = Dh2b
coupled_variables = eta2
derivative_order = 1
[]
[Dh2c]
type = DerivativeParsedMaterial
material_property_names = 'D h2(eta2)'
expression = D*h2*16/30
property_name = Dh2c
coupled_variables = eta2
derivative_order = 1
[]
[]
[Kernels]
#Kernels for diffusion equation
[diff_time]
type = TimeDerivative
variable = cCr
[]
[diff_c1]
type = MatDiffusion
variable = cCr
diffusivity = Dh1
v = BCC_CR
args = eta1
[]
[diff_c2a]
type = MatDiffusion
variable = cCr
diffusivity = Dh2a
v = SIGMA_0CR
args = eta2
[]
[diff_c2b]
type = MatDiffusion
variable = cCr
diffusivity = Dh2b
v = SIGMA_1CR
args = eta2
[]
[diff_c2c]
type = MatDiffusion
variable = cCr
diffusivity = Dh2c
v = SIGMA_2CR
args = eta2
[]
# enforce pointwise equality of chemical potentials
[chempot1a2a]
# The BCC phase has only one sublattice
# we tie it to the first sublattice with site fraction 10/(10+4+16) in the sigma phase
type = KKSPhaseChemicalPotential
variable = BCC_CR
cb = SIGMA_0CR
kb = '${fparse 10/30}'
fa_name = F_BCC_A2
fb_name = F_SIGMA
args_b = 'SIGMA_1CR SIGMA_2CR'
[]
[chempot2a2b]
# This kernel ties the first two sublattices in the sigma phase together
type = SLKKSChemicalPotential
variable = SIGMA_0CR
a = 10
cs = SIGMA_1CR
as = 4
F = F_SIGMA
coupled_variables = 'SIGMA_2CR'
[]
[chempot2b2c]
# This kernel ties the remaining two sublattices in the sigma phase together
type = SLKKSChemicalPotential
variable = SIGMA_1CR
a = 4
cs = SIGMA_2CR
as = 16
F = F_SIGMA
coupled_variables = 'SIGMA_0CR'
[]
[phaseconcentration]
# This kernel ties the sum of the sublattice concentrations to the global concentration cCr
type = SLKKSMultiPhaseConcentration
variable = SIGMA_2CR
c = cCr
ns = '1 3'
as = '1 10 4 16'
cs = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR'
h_names = 'h1 h2'
eta = 'eta1 eta2'
[]
# Kernels for Allen-Cahn equation for eta1
[deta1dt]
type = TimeDerivative
variable = eta1
[]
[ACBulkF1]
type = KKSMultiACBulkF
variable = eta1
Fj_names = 'F_BCC_A2 F_SIGMA'
hj_names = 'h1 h2'
gi_name = g1
eta_i = eta1
wi = 0.1
coupled_variables = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR eta2'
[]
[ACBulkC1]
type = SLKKSMultiACBulkC
variable = eta1
F = F_BCC_A2
c = BCC_CR
ns = '1 3'
as = '1 10 4 16'
cs = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR'
h_names = 'h1 h2'
eta = 'eta1 eta2'
[]
[ACInterface1]
type = ACInterface
variable = eta1
kappa_name = kappa
[]
[lagrange1]
type = SwitchingFunctionConstraintEta
variable = eta1
h_name = h1
lambda = lambda
coupled_variables = 'eta2'
[]
# Kernels for Allen-Cahn equation for eta1
[deta2dt]
type = TimeDerivative
variable = eta2
[]
[ACBulkF2]
type = KKSMultiACBulkF
variable = eta2
Fj_names = 'F_BCC_A2 F_SIGMA'
hj_names = 'h1 h2'
gi_name = g2
eta_i = eta2
wi = 0.1
coupled_variables = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR eta1'
[]
[ACBulkC2]
type = SLKKSMultiACBulkC
variable = eta2
F = F_BCC_A2
c = BCC_CR
ns = '1 3'
as = '1 10 4 16'
cs = 'BCC_CR SIGMA_0CR SIGMA_1CR SIGMA_2CR'
h_names = 'h1 h2'
eta = 'eta1 eta2'
[]
[ACInterface2]
type = ACInterface
variable = eta2
kappa_name = kappa
[]
[lagrange2]
type = SwitchingFunctionConstraintEta
variable = eta2
h_name = h2
lambda = lambda
coupled_variables = 'eta1'
[]
# Lagrange-multiplier constraint kernel for lambda
[lagrange]
type = SwitchingFunctionConstraintLagrange
variable = lambda
h_names = 'h1 h2'
etas = 'eta1 eta2'
epsilon = 1e-6
[]
[]
[AuxKernels]
[GlobalFreeEnergy]
type = KKSMultiFreeEnergy
variable = Fglobal
Fj_names = 'F_BCC_A2 F_SIGMA'
hj_names = 'h1 h2'
gj_names = 'g1 g2'
interfacial_vars = 'eta1 eta2'
kappa_names = 'kappa kappa'
w = 0.1
[]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
line_search = none
petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -ksp_gmres_restart'
petsc_options_value = 'asm lu nonzero 30'
l_max_its = 100
nl_max_its = 20
nl_abs_tol = 1e-10
end_time = 1000
[TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 12
iteration_window = 2
growth_factor = 2
cutback_factor = 0.5
dt = 0.1
[]
[]
[Postprocessors]
[F]
type = ElementIntegralVariablePostprocessor
variable = Fglobal
[]
[cmin]
type = NodalExtremeValue
value_type = min
variable = cCr
[]
[cmax]
type = NodalExtremeValue
value_type = max
variable = cCr
[]
[]
[Outputs]
exodus = true
print_linear_residuals = false
# exclude lagrange multiplier from output, it can diff more easily
hide = lambda
[]
(modules/phase_field/test/tests/KKS_system/lagrange_multiplier.i)
#
# This test ensures that the equilibrium solution using two order parameters with a
# Lagrange multiplier constraint is identical to the dedicated two phase formulation
# in two_phase.i
#
[Mesh]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 5
[]
[AuxVariables]
[Fglobal]
order = CONSTANT
family = MONOMIAL
[]
[]
[Variables]
# concentration
[c]
order = FIRST
family = LAGRANGE
[InitialCondition]
type = FunctionIC
function = x/5
[]
[]
# order parameter 1
[eta1]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
# order parameter 2
[eta2]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
# phase concentration 1
[c1]
order = FIRST
family = LAGRANGE
initial_condition = 0.2
[]
# phase concentration 2
[c2]
order = FIRST
family = LAGRANGE
initial_condition = 0.5
[]
# Lagrange multiplier
[lambda]
order = FIRST
family = LAGRANGE
initial_condition = 0.0
[]
[]
[Materials]
# simple toy free energies
[f1] # = fd
type = DerivativeParsedMaterial
property_name = F1
coupled_variables = 'c1'
expression = '(0.9-c1)^2'
[]
[f2] # = fm
type = DerivativeParsedMaterial
property_name = F2
coupled_variables = 'c2'
expression = '(0.1-c2)^2'
[]
# Switching functions for each phase
[h1_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta1
function_name = h1
[]
[h2_eta]
type = SwitchingFunctionMaterial
h_order = HIGH
eta = eta2
function_name = h2
[]
# Coefficients for diffusion equation
[Dh1]
type = DerivativeParsedMaterial
material_property_names = 'D h1'
expression = D*h1
property_name = Dh1
[]
[Dh2]
type = DerivativeParsedMaterial
material_property_names = 'D h2'
expression = D*h2
property_name = Dh2
[]
# Barrier functions for each phase
[g1]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta1
function_name = g1
[]
[g2]
type = BarrierFunctionMaterial
g_order = SIMPLE
eta = eta2
function_name = g2
[]
# constant properties
[constants]
type = GenericConstantMaterial
prop_names = 'D L kappa'
prop_values = '0.7 0.7 0.2'
[]
[]
[Kernels]
#Kernels for diffusion equation
[diff_time]
type = TimeDerivative
variable = c
[]
[diff_c1]
type = MatDiffusion
variable = c
diffusivity = Dh1
v = c1
[]
[diff_c2]
type = MatDiffusion
variable = c
diffusivity = Dh2
v = c2
[]
# Kernels for Allen-Cahn equation for eta1
[deta1dt]
type = TimeDerivative
variable = eta1
[]
[ACBulkF1]
type = KKSMultiACBulkF
variable = eta1
Fj_names = 'F1 F2 '
hj_names = 'h1 h2 '
gi_name = g1
eta_i = eta1
wi = 0.2
coupled_variables = 'c1 c2 eta2'
[]
[ACBulkC1]
type = KKSMultiACBulkC
variable = eta1
Fj_names = 'F1 F2'
hj_names = 'h1 h2'
cj_names = 'c1 c2'
eta_i = eta1
coupled_variables = 'eta2'
[]
[ACInterface1]
type = ACInterface
variable = eta1
kappa_name = kappa
[]
[multipler1]
type = MatReaction
variable = eta1
v = lambda
mob_name = L
[]
# Kernels for the Lagrange multiplier equation
[mult_lambda]
type = MatReaction
variable = lambda
mob_name = 2
[]
[mult_ACBulkF_1]
type = KKSMultiACBulkF
variable = lambda
Fj_names = 'F1 F2 '
hj_names = 'h1 h2 '
gi_name = g1
eta_i = eta1
wi = 0.2
mob_name = 1
coupled_variables = 'c1 c2 eta2 '
[]
[mult_ACBulkC_1]
type = KKSMultiACBulkC
variable = lambda
Fj_names = 'F1 F2'
hj_names = 'h1 h2'
cj_names = 'c1 c2'
eta_i = eta1
coupled_variables = 'eta2 '
mob_name = 1
[]
[mult_CoupledACint_1]
type = SimpleCoupledACInterface
variable = lambda
v = eta1
kappa_name = kappa
mob_name = 1
[]
[mult_ACBulkF_2]
type = KKSMultiACBulkF
variable = lambda
Fj_names = 'F1 F2 '
hj_names = 'h1 h2 '
gi_name = g2
eta_i = eta2
wi = 0.2
mob_name = 1
coupled_variables = 'c1 c2 eta1 '
[]
[mult_ACBulkC_2]
type = KKSMultiACBulkC
variable = lambda
Fj_names = 'F1 F2'
hj_names = 'h1 h2'
cj_names = 'c1 c2'
eta_i = eta2
coupled_variables = 'eta1 '
mob_name = 1
[]
[mult_CoupledACint_2]
type = SimpleCoupledACInterface
variable = lambda
v = eta2
kappa_name = kappa
mob_name = 1
[]
# Kernels for constraint equation eta1 + eta2 = 1
# eta2 is the nonlinear variable for the constraint equation
[eta2reaction]
type = MatReaction
variable = eta2
mob_name = 1
[]
[eta1reaction]
type = MatReaction
variable = eta2
v = eta1
mob_name = 1
[]
[one]
type = BodyForce
variable = eta2
value = -1.0
[]
# Phase concentration constraints
[chempot12]
type = KKSPhaseChemicalPotential
variable = c1
cb = c2
fa_name = F1
fb_name = F2
[]
[phaseconcentration]
type = KKSMultiPhaseConcentration
variable = c2
cj = 'c1 c2'
hj_names = 'h1 h2'
etas = 'eta1 eta2'
c = c
[]
[]
[AuxKernels]
[Fglobal_total]
type = KKSMultiFreeEnergy
Fj_names = 'F1 F2 '
hj_names = 'h1 h2 '
gj_names = 'g1 g2 '
variable = Fglobal
w = 0.2
interfacial_vars = 'eta1 eta2 '
kappa_names = 'kappa kappa'
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type -sub_pc_factor_shift_type'
petsc_options_value = 'lu nonzero'
l_max_its = 30
nl_max_its = 10
l_tol = 1.0e-4
nl_rel_tol = 1.0e-10
nl_abs_tol = 1.0e-11
num_steps = 35
dt = 10
[]
[VectorPostprocessors]
[c]
type = LineValueSampler
variable = c
start_point = '0 0 0'
end_point = '5 0 0'
num_points = 21
sort_by = x
[]
[]
[Outputs]
csv = true
execute_on = FINAL
[]