Input File Changes:

We are going to make two changes to our simulation. First, we will make the mobility a function instead of a constant. Then, we will check if the surface is decomposing to 39.5% chromium phase as the material balance says it should. Both of these changes are primarily in the materials block. We will use a DerivativeParsedMaterial to define the equation for mobility, similar to how we defined the free energy density. We will have to remove it from the constants block. Then we will add another material to check if each element is in the chromium phase. The materials block now looks like this:

  # 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 first three
  # materials and must have the same value in each one.
  [./kappa]                  # Gradient energy coefficient (eV nm^2/mol)
    type = GenericFunctionMaterial
    prop_names = 'kappa_c'
    prop_values = '8.125e-16*6.24150934e+18*1e+09^2*1e-27'
                  # kappa_c*eV_J*nm_m^2*d
  [./mobility]               # Mobility (nm^2 mol/eV/s)
    type = DerivativeParsedMaterial
    f_name = M
    args = c
    constant_names =       'Acr    Bcr    Ccr    Dcr
                            Ecr    Fcr    Gcr
                            Afe    Bfe    Cfe    Dfe
                            Efe    Ffe    Gfe
                            nm_m   eV_J   d'
    constant_expressions = '-32.770969 -25.8186669 -3.29612744 17.669757
                            37.6197853 20.6941796  10.8095813
                            -31.687117 -26.0291774 0.2286581   24.3633544
                            44.3334237 8.72990497  20.956768
                            1e+09      6.24150934e+18          1e-27'
    function = 'nm_m^2/eV_J/d*((1-c)^2*c*10^
    derivative_order = 1
    outputs = exodus
  [./local_energy]           # Local free energy function (eV/mol)
    type = DerivativeParsedMaterial
    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)+
    derivative_order = 2
  [./precipitate_indicator]  # Returns 1/625 if precipitate
      type = ParsedMaterial
      f_name = prec_indic
      args = c
      function = if(c>0.6,0.0016,0)

Note that the area of our surface is 625 $$$nm^2$$$, and $$$\frac {1} {625} = 0.0016$$$. Within the different materials, the argument "block = 0" is not included. It can save time if we can define common arguments once rather than in each material. To do this, we include a global parameters block.

  block = 0           # The generated mesh is used for all materials and kernels

The mobility equation is now ready to run. However, the chromium phase fraction will not be output. In order to output it we add a post-processor.

  [./precipitate_area]      # Fraction of surface devoted to precipitates
    type = ElementIntegralMaterialProperty
    mat_prop = prec_indic

This brings the number of post-processors up to 6.


simulation result

chromium phase fraction

The first figure to the right shows the results of this simulation. As you can see, changing the mobility to a function reduced the number of features. The second figure shows the fraction of the surface covered with the chromium phase as a function of time. It is very close to the 39.5% we calculated, which is promising. It is not exact because of the concentration gradients at the phase interfaces.

If you still have the residual magnitudes outputting to the screen, you may have noticed that the magnitudes changed this time. The residual for the concentration is now several orders of magnitude smaller than the residual for chemical potential. In our next simulation we will scale the concentration variable so that this is not the case.


We are almost done. Next, we will calculate the total free energy of the surface and look for its s-curve.

step 5: Check the Surface Energy Curve