Chemical potential calcluted with Split formulation Cahn-Hilliard Kernel

[ZhouXD1989]

Hello MOOSE members,

I am using the Split formulation Cahn-Hilliard Kernel to simulate coupled chemo-mechanical diffusing problems (with kappa=0) with additional eigenstrains, e.g. dislocations.

When I consider the additional eigenstrain, the calculated chemical potential 'w' is quite different from the material properties directly calculated from mu = dF/dc (saved in axuvariable 'mu' through MaterialRealAux), as shown in the following figure. If there is no additional eigenstrain, the solutions of 'w' and 'mu' are the same.
So what could be the reason for this difference?

Following is my input file:

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = -25
  xmax = 25
  ymin = -25
  ymax = 25
  nx = 50
  ny = 50
  elem_type = QUAD4
[]

[Variables]
  [./c]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.05
  [../]
  [./w]
    order = FIRST
    family = LAGRANGE
  [../]
  [./disp_x]
    order = FIRST
    family = LAGRANGE
  [../]
  [./disp_y]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[Kernels]
  [./TensorMechanics]
    displacements = 'disp_x disp_y'
  [../]
  [./c_res]
    type = SplitCHParsed
    variable = c
    f_name = F
    kappa_name = kappa_c
    w = w
  [../]
  [./w_res]
    type = SplitCHWRes
    variable = w
    mob_name = M
    coupled_variables = c
  [../]
  [./time]
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]
[]

[BCs]
  [./flux]
    type = FunctionNeumannBC
    variable = w
    function = 0.1
    boundary = 'bottom'
  [../]
[]

[AuxVariables]
  [./mu]
    order = FIRST
    family = MONOMIAL
  [../]
  [./d]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./sigma11_aux]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./sigma22_aux]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./sigma12_aux]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]
[AuxKernels]
  [./mu]
    type = MaterialRealAux
    property = dF_dc 
    variable = mu
  [../]
  [./d]
    type = FunctionAux
    variable = d
    function = 'if(x>-30&x<0.1&y>-0.6&y<0.6,1,0)'
  [../]
  [./matl_sigma11]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 0
    index_j = 0
    variable = sigma11_aux
  [../]
  [./matl_sigma22]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 1
    index_j = 1
    variable = sigma22_aux
  [../]
  [./matl_sigma12]
    type = RankTwoAux
    rank_two_tensor = stress
    index_i = 0
    index_j = 1
    variable = sigma11_aux
  [../]
[]

[Materials]
  [./interface]
    type = GenericConstantMaterial
    prop_names  = 'kappa_c'
    prop_values = '0'
  [../]
  [./mobility]
    type = DerivativeParsedMaterial
    property_name = 'M'
    coupled_variables = 'c'
    constant_names       = 'D '
    constant_expressions = '1000'
    expression = 'D*c*(1-c)'
    enable_jit = true
    derivative_order = 1
  [../]
  [./chemical_free_energy]
    type = DerivativeParsedMaterial
    property_name = Fc
    coupled_variables = 'c'
    constant_names       = 'X '
    constant_expressions = '1'
    expression = 'c*log(c)+(1-c)*log(1-c)+X*c*(1-c)'
    enable_jit = true
    derivative_order = 2
  [../]
  [./elastic_free_energy]
    type = ElasticEnergyMaterial
    f_name = Fe
    args = 'c'
    derivative_order = 2
  [../]
  [./free_energy]
    type = DerivativeSumMaterial
    property_name = F
    sum_materials = 'Fc Fe'
    coupled_variables = 'c'
    derivative_order = 2
  [../]
  [./dFdc]
    type = ParsedMaterial
    property_name = dF_dc
    material_property_names = 'dF:=D[F,c]'
    expression = 'dF'
  [../]
  [./Cijkl_alpha]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1523.1
    poissons_ratio = 0.3
  [../]
  [./stress]
    type = ComputeLinearElasticStress
  [../]
  [./strain]
    type = ComputeSmallStrain
    displacements = 'disp_x disp_y'
    eigenstrain_names = 'eigenstrain_c eigenstrain_D'
  [../]
  [./var_dependence]
    type = DerivativeParsedMaterial
    expression = 0.08/3*c
    coupled_variables = c
    f_name = var_dep
    enable_jit = true
    derivative_order = 2
  [../]
  [./eigenstrain]
    type = ComputeVariableEigenstrain
    eigen_base = '1 1 1 0 0 0'
    prefactor = var_dep
    args = 'c'
    eigenstrain_name = eigenstrain_c
  [../]
  [./var_dis]
    type = DerivativeParsedMaterial
    expression = d
    coupled_variables = d
    f_name = var_dis
    enable_jit = true
    derivative_order = 2
  [../]
  [./Dislocation_Eigenstrain]
    type = ComputeVariableEigenstrain
    eigen_base = '0 0 0 0 0 0.2'
    prefactor = var_dis
    args = 'd'
    eigenstrain_name = eigenstrain_D
  [../]
[]

[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 = 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 = 10000

  [./TimeStepper]
    type = SolutionTimeAdaptiveDT
    dt = 1
  [../]
[]

[Outputs]
  exodus = true
[]

[GiudGiud]

@laagesen @amjokisaari