Compute Surface Tension KKS

Surface tension of an interface defined by the gradient of an order parameter

Description

This material generates a surface tension tensor in the interface between two phases in the Kim-Kim-Suzuki (KKS) phase-field model (Kim et al. (1999)). The surface tension consists of two orthogonal principal components lying in plane of the interface. The approach is based on the formulation described in Levitas and Samani (2011). (For solid surfaces where the surface stress is to be specified rather than surface tension, use the ComputeInterfaceStress kernel instead.)

The surface tension tensor $var element = document.getElementById("moose-equation-5e72b0ec-6301-40f9-b945-747fe294f231");katex.render("\\mathbf{\\sigma}_{st}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is given by $var element = document.getElementById("moose-equation-ac83072f-94ba-44e8-8ecd-b7425aec2619");katex.render("\\mathbf{\\sigma}_{st} = \\left[W g(\\eta) + \\frac \\kappa 2 |\\nabla \\eta|^2 \\right] \\bf{I} - \\kappa \\nabla \\eta \\otimes \\nabla \\eta", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ where $var element = document.getElementById("moose-equation-3b16a3ef-d744-4ce7-8470-b4fdb5cc3e83");katex.render("W", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the free energy barrier between phases, $var element = document.getElementById("moose-equation-d504324c-1f32-4a1a-a4c5-25cd0827080f");katex.render("g(\\eta)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the double-well function, $var element = document.getElementById("moose-equation-7733489d-57fb-4a26-8f3a-492fbb8fe44f");katex.render("\\eta", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the order parameter, $var element = document.getElementById("moose-equation-ec3210c2-1ecc-4520-b68b-1eca584c0df5");katex.render("\\kappa", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the gradient energy coefficient, $var element = document.getElementById("moose-equation-1bc54821-302f-4c8e-bdb1-58bbe27aa8dd");katex.render("\\bf{I}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the identity tensor, and $var element = document.getElementById("moose-equation-ed71f8a6-2925-4315-b95b-7cac939951ee");katex.render("\\otimes", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ denotes the dyadic product.

Example Input File Syntax

[./interface_stress]
  type = ComputeSurfaceTensionKKS
  v = eta
  kappa_name = kappa
  w = 0.0033
[../]

Input Parameters

wDouble well height parameter
C++ Type:double
Controllable:No
Description:Double well height parameter

Required Parameters

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
Options:NONE, ELEMENT, SUBDOMAIN
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.
ggBarrier Function Material that provides g(eta)
Default:g
C++ Type:MaterialPropertyName
Controllable:No
Description:Barrier Function Material that provides g(eta)
kappa_namekappa_opGradient energy coefficient
Default:kappa_op
C++ Type:MaterialPropertyName
Controllable:No
Description:Gradient energy coefficient
planar_stress_nameextra_stressMaterial property name for the interfacial planar stress
Default:extra_stress
C++ Type:MaterialPropertyName
Controllable:No
Description:Material property name for the interfacial planar stress
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.
vOrder parameter that defines the interface, assumed to vary from 0 to 1.
C++ Type:std::vector<VariableName>
Controllable:No
Description:Order parameter that defines the interface, assumed to vary from 0 to 1.

Optional 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/test/tests/surface_tension_KKS/surface_tension_KKS.i)
(modules/combined/test/tests/surface_tension_KKS/surface_tension_VDWgas.i)

References

Seong Gyoon Kim, Won Tae Kim, and Toshio Suzuki. Phase-field model for binary alloys. Physical Review E, 60(6):7186–7197, December 1999. URL: http://link.aps.org/doi/10.1103/PhysRevE.60.7186 (visited on 2014-03-31), doi:10.1103/PhysRevE.60.7186.

@article{kim_phase-field_1999,
    author = "Kim, Seong Gyoon and Kim, Won Tae and Suzuki, Toshio",
    title = "Phase-field model for binary alloys",
    volume = "60",
    url = "http://link.aps.org/doi/10.1103/PhysRevE.60.7186",
    doi = "10.1103/PhysRevE.60.7186",
    abstract = "We present a phase-field model ({PFM}) for solidification in binary alloys, which is found from the phase-field model for a pure material by direct comparison of the variables for a pure material solidification and alloy solidification. The model appears to be equivalent with the Wheeler-Boettinger-{McFadden} ({WBM}) model [A.A. Wheeler, W. J. Boettinger, and G. B. {McFadden}, Phys. Rev. A 45, 7424 (1992)], but has a different definition of the free energy density for interfacial region. An extra potential originated from the free energy density definition in the {WBM} model disappears in this model. At a dilute solution limit, the model is reduced to the Tiaden et al. model [Physica D 115, 73 (1998)] for a binary alloy. A relationship between the phase-field mobility and the interface kinetics coefficient is derived at a thin-interface limit condition under an assumption of negligible diffusivity in the solid phase. For a dilute alloy, a steady-state solution of the concentration profile across the diffuse interface is obtained as a function of the interface velocity and the resultant partition coefficient is compared with the previous solute trapping model. For one dimensional steady-state solidification, where the classical sharp-interface model is exactly soluble, we perform numerical simulations of the phase-field model: At low interface velocity, the simulated results from the thin-interface {PFM} are in excellent agreement with the exact solutions. As the partition coefficient becomes close to unit at high interface velocities, whereas, the sharp-interface {PFM} yields the correct answer.",
    number = "6",
    urldate = "2014-03-31",
    journal = "Physical Review E",
    month = "December",
    year = "1999",
    pages = "7186--7197",
    annote = "This is the {KKS} model paper",
    file = "Kim et al. - 1999 - Phase-field model for binary alloys.html:/Users/schwd/Library/Application Support/Firefox/Profiles/i2j5gt6k.default/zotero/storage/IQMMGVCR/Kim et al. - 1999 - Phase-field model for binary alloys.html:text/html;Kim et al. - 1999 - Phase-field model for binary alloys.pdf:/Users/schwd/Library/Application Support/Firefox/Profiles/i2j5gt6k.default/zotero/storage/R6VMMT2S/Kim et al. - 1999 - Phase-field model for binary alloys.pdf:application/pdf"
}

V. I. Levitas and K. Samani. Size and mechanics effects in surface-induced melting of nanoparticles. Nature Communications, 2:284, 2011. URL: https://www.nature.com/articles/ncomms1275, doi:10.1038/ncomms1275.

@article{Levitas2011,
    author = "Levitas, V. I. and Samani, K.",
    journal = "Nature Communications",
    url = "https://www.nature.com/articles/ncomms1275",
    doi = "10.1038/ncomms1275",
    title = "Size and mechanics effects in surface-induced melting of nanoparticles",
    volume = "2",
    pages = "284",
    year = "2011"
}

(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/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/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
[]

Description
Example Input File Syntax
Input Parameters
Input Files
References