Input File Blocks:

We are going to start with the simplest input file we can. We will add on more features and make the model more accurate as we go. We will be using the split Cahn-Hilliard method because it converges faster than the direct method. The basic blocks we need are: Mesh, Variables, ICs, BCs, Kernels, Materials, Precoditioning, Executioner, and Outputs.


The mesh block is going to create a two dimensional mesh that is 25nm × 25nm and has 100 × 100 elements.

  # generate a 2D, 25nm x 25nm mesh
  type = GeneratedMesh
  dim = 2
  elem_type = QUAD4
  nx = 100
  ny = 100
  nz = 0
  xmin = 0
  xmax = 25
  ymin = 0
  ymax = 25
  zmin = 0
  zmax = 0


Because we are using the split form, two variables are defined, the mole fraction of chromium and the chemical potential.

  [./c]   # Mole fraction of Cr (unitless)
    order = FIRST
    family = LAGRANGE
  [./w]   # Chemical potential (eV/mol)
    order = FIRST
    family = LAGRANGE


To start off we want to check our second expectation. That is, that the surface will minimize its interface region by forming circular or straight line regions. To do this, we will use a Bounding Box Initial Condition with the equilibrium concentrations we calculated from the free energy density curve.

  # Use a bounding box IC at equilibrium concentrations to make sure the
  # model behaves as expected.
    type = BoundingBoxIC
    variable = c
    x1 = 5
    x2 = 20
    y1 = 5
    y2 = 20
    inside = 0.823
    outside = 0.236


It is common in phase field simulations to use periodic boundary conditions.

  # periodic BC as is usually done on phase-field models
      auto_direction = 'x y'


We will use the kernels specified for the split Cahn-Hilliard equation. Note that here we name the free energy function, the mobility, and the gradient energy coefficient "f_loc", "M", and "kappa_c" respectively. These have the same names in the materials block.

    variable = w
    v = c
    type = CoupledTimeDerivative
    variable = w
    type = SplitCHWRes
    mob_name = M
    variable = c
    type = SplitCHParsed
    f_name = f_loc
    kappa_name = kappa_c
    w = w


The materials block defines the equations and constants in the model. This is where we will define $$$f_{loc}(c)$$$, $$$M(c)$$$, and $$$\kappa$$$ in MOOSE. To simplify things right now, we are going to define $$$M(c)$$$ as a constant value rather than a function. We will also be doing unit conversions in the materials block. Note that each value has the same name as in the kernels block.

In addition to the unit conversions, it was accidentally discovered that the convergence properties can be improved by using a scaling factor, $$$d$$$. Both $$$f_{loc}(c)$$$ and $$$\kappa$$$ are multiplied by $$$d$$$, while $$$M(c)$$$ is multiplied by the inverse of $$$d$$$. This does not change the solution, but it does seem to help the convergence.

  # d is a scaling factor that makes it easier for the solution to converge
  # without changing the results. It is defined in each of the materials and
  # must have the same value in each one.
    # Define constant values kappa_c and M. Eventually M will be replaced with
    # an equation rather than a constant.
    type = GenericFunctionMaterial
    block = 0
    prop_names = 'kappa_c M'
    prop_values = '8.125e-16*6.24150934e+18*1e+09^2*1e-27
                   # kappa_c*eV_J*nm_m^2*d
                   # M*nm_m^2/eV_J/d
    # Defines the function for the local free energy density as given in the
    # problem, then converts units and adds scaling factor.
    type = DerivativeParsedMaterial
    block = 0
    f_name = f_loc
    args = c
    constant_names = 'A   B   C   D   E   F   G  eV_J  d'
    constant_expressions = '-2.446831e+04 -2.827533e+04 4.167994e+03 7.052907e+03
                            1.208993e+04 2.568625e+03 -2.354293e+03
                            6.24150934e+18 1e-27'
    function = 'eV_J*d*(A*c+B*(1-c)+C*c*log(c)+D*(1-c)*log(1-c)+


The split Cahn-Hilliard has the best convergence properties when we use the Newton solver. The newton solver always requires preconditioning.

    type = SMP
    full = true


The executioner block tells MOOSE how to solve the problem. Many of the values in this block are arbitrary and you can experiment with changing them to see how they change the solution. More information on the executioner and preconditioning can be found on the wiki.

Right now our main interest is in seeing if the model is working. So we shortened the time that the simulation runs to just be long enough to see if the grain is becoming circular.

  type = Transient
  solve_type = NEWTON
  l_max_its = 30
  l_tol = 1e-6
  nl_max_its = 50
  nl_abs_tol = 1e-9
  end_time = 86400   # 1 day. We only need to run this long enough to verify
                     # the model is working properly.
  petsc_options_iname = '-pc_type -ksp_grmres_restart -sub_ksp_type
                         -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  preonly
                         ilu          1'
  dt = 100


The outputs block lets us decide what MOOSE tells us about the simulation. We will output the initial condition so we can compare it to the final condition.

  exodus = true
  console = true
  print_perf_log = true
  output_initial = true

Simulation Results:

initial condition of Simple Test Model

final result of Simple Test Model

The first image to the right shows the initial condition of this simulation. The second image shows the end result of the simulation. The result shows that the chromium phase's corners are rounding out and the shape is becoming more circular. This is exactly what we would expect to happen as the surface tries to reduce its free energy. This result is a promising sign that our simulation is working correctly.


step 2: Make a Faster Test Model