Discrete Nucleation Module inserting nuclei with negative value order parameter in grand potential model #29081

aclin4 · 2024-11-14T20:01:45Z

aclin4
Nov 14, 2024

Hi all,

I'm working on a grand potential phase field model to simulate the formation of gas bubbles in a fuel matrix. In my model there are four variables: the matrix phase order parameter, the bubble phase order parameter, chemical potential of gas, and chemical potential of vacancies. The full input file is included at the end of this post.

I'm inserting the nuclei using the DiscreteNucleationForce kernel acting on the bubble order parameter:

[nucleation_b]
    type = DiscreteNucleationForce
    variable = etab0
    map = B_map
    nucleus_value = 100
    no_nucleus_value = 0
[]

where the nucleation probability is based on the supersaturation of the gas and vacancies in the matrix:

[G_nuc]
   type = DerivativeParsedMaterial
   property_name = G_nuc
   material_property_names = 'km cveq_m cgeq_m cgeq_b cv(wv) cg(wg)'
   expression = '(cv - cveq_m) + (cg-cgeq_m)'
   outputs = exodus
   output_properties = G_nuc
 []
 [nucl_matl]
   type = ParsedMaterial
   property_name = nucl_matl
   coupled_variables = 'etab0'
   material_property_names = 'cg(wg) T G_nuc'
   constant_names = 'k1star k2'
   constant_expressions = '1e-6 	10'
   expression = 'if(cg>1.5e-6,if(etab0<0.0001, k1star * exp(-k2/G_nuc/T), 0),0)'
   outputs = exodus
 []
[../]

The UserObjects for the Discrete Nucleation Module are set up like so:

[UserObjects]
  [B_inserter]
    type = DiscreteNucleationInserter
    hold_time = 1000
    probability = nucl_matl
    radius = 2
  []
  [B_map]
    type = DiscreteNucleationMap
    periodic = etab0
    inserter = B_inserter
    int_width = 0.5
  []
[]

However, when I run my simulation and a nuclei is inserted, the order parameter for the bubble phase takes a negative value. (case A)

In certain instances, if I modify the nucleation rate to be higher at a lower supersaturation, the nuclei will be inserted initially with a positive value, but in the subsequent timestep, the value becomes negative. (case B)

Does anyone have any insight as to why this is happening?

Update
I managed to get a positive value inserted for the nuclei by using the DiscreteNucleationTimeStep to limit the timestep near the nuclei insertion, however, now I get a strange propagating numerical error surrounding the nuclei (see Case C figures below). I'm unsure if this is related to the above issue or is a completely different numerical problem. Any insight into this would be greatly appreciated!

Thanks,
Albert

Case A:

Case B:

Case C:

Full Input File:

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    elem_type = QUAD4
    nx = 80
    ny = 80
    nz = 0
    xmin = -20
    xmax = 20
    ymin = -20
    ymax = 20
    zmin = 0
    zmax = 0
    #uniform_refine = 2
  []
[]
[GlobalParams]
  file_base = 'just_nuc'
[]

[Variables]
  [./wv]
  [../]
  [./wg]
  [../]
  [./etam0]
  [../]
  [./etab0]
  [../]
[]

[AuxVariables]
  [./n]
    order = CONSTANT
    family = MONOMIAL
  []
[]
[AuxKernels]
  [./nucc]
    type = DiscreteNucleationAux
    map = B_map
    variable = n
    no_nucleus_value = 0
    nucleus_value = 10
    execute_on = TIMESTEP_END
  [../]
[]

[BCs]
  [./Periodic]
    [./all]
      variable = 'etam0 etab0 wv wg'
      auto_direction = 'x y'
    [../]
  [../]

[]

[ICs]
  [./IC_etam0]
    type = SpecifiedSmoothCircleIC
    variable = etam0
    invalue = 1
    outvalue = 1
    x_positions = '  -8.0   7.0 -10.0   7.0'
    y_positions = '  -8.0   8.0   6.0  -8.0'
    z_positions = '   0.0   0.0   0.0   0.0'
    radii = '   1.0   1.0   1.0   1.0'
  []
  [./IC_etab0]
    type = SpecifiedSmoothCircleIC
    variable = etab0
    invalue = 0
    outvalue = 0
    x_positions = '  -8.0   7.0 -10.0   7.0'
    y_positions = '  -8.0   8.0   6.0  -8.0'
    z_positions = '   0.0   0.0   0.0   0.0'
    radii = '   1.0   1.0   1.0   1.0'
  []
  [./IC_wv]
    type = RandomIC
    variable = wv
    min  = 5
    max  = 6
    seed = 349
  [../]
  [./IC_wg]
    type = RandomIC
    variable = wg
    min = 5
    max = 6
    seed = 234
  [../]
[]

[Kernels]
# Order parameter eta_m0
  [./em0_dot]
    type = TimeDerivative
    variable = etam0
  [../]
  [./ACm0_bulk]
    type = ACGrGrMulti
    variable = etam0
    v =           'etab0'
    gamma_names = 'gmb'
    mob_name = L
  [../]
  [./ACm0_int]
    type = ACInterface
    variable = etam0
    kappa_name = kappa
  [../]
  [./ACm0_sw]
    type = ACSwitching
    variable = etam0
    Fj_names  = 'omegam omegab'
    hj_names  = 'hm     hb'
    coupled_variables = 'etab0 wv wg'
  [../]
  # [nucleation_m]
  #   type = DiscreteNucleationForce
  #   variable = etam0
  #   map = B_map
  #   nucleus_value = -5
  #   no_nucleus_value = 0
  # []
 # Order parameter eta_b0
  [./eb0_dot]
    type = TimeDerivative
    variable = etab0
  [../]
  [./ACb0_bulk]
    type = ACGrGrMulti
    variable = etab0
    v =           'etam0'
    gamma_names = 'gmb'
    mob_name = L
  [../]
  [./ACb0_int]
    type = ACInterface
    variable = etab0
    kappa_name = kappa
  [../]
  [./ACb0_sw]
    type = ACSwitching
    variable = etab0
    Fj_names  = 'omegam omegab'
    hj_names  = 'hm     hb'
    coupled_variables = 'etam0 wv wg'
  [../]
  [nucleation_b]
    type = DiscreteNucleationForce
    variable = etab0
    map = B_map
    nucleus_value = 100
    no_nucleus_value = 0
  []
#Chemical potential vacancies
  [./wv_dot]
    type = SusceptibilityTimeDerivative
    variable = wv
    f_name = chiv
    coupled_variables = 'etam0 etab0'
  [../]
  [./wvDiffusion]
    type = MatDiffusion
    variable = wv
    diffusivity = Dchiv
    args = ''
  [../]
  [./wvcoupled_etam0dot]
    type = CoupledSwitchingTimeDerivative
    variable = wv
    v = etam0
    Fj_names = 'rhovm rhovb'
    hj_names = 'hm   hb'
    coupled_variables = 'etam0 etab0'
  [../]
  [./wvcoupled_etab0dot]
    type = CoupledSwitchingTimeDerivative
    variable = wv
    v = etab0
    Fj_names = 'rhovm rhovb'
    hj_names = 'hm   hb'
    coupled_variables = 'etam0 etab0'
  [../]
  [vac_prod]
    type = MaskedBodyForce
    variable = wv
    mask = prod
    value = 0.2156
  []
#Chemical potential gas
  [./wg_dot]
    type = SusceptibilityTimeDerivative
    variable = wg
    f_name = chig
    coupled_variables = 'etam0 etab0'
  [../]
  [./wgDiffusion]
    type = MatDiffusion
    variable = wg
    diffusivity = Dchig
    args = ''
  [../]
  [./wgcoupled_etam0dot]
    type = CoupledSwitchingTimeDerivative
    variable = wg
    v = etam0
    Fj_names = 'rhogm rhogb'
    hj_names = 'hm   hb'
    coupled_variables = 'etam0 etab0'
  [../]
  [./wgcoupled_etab0dot]
    type = CoupledSwitchingTimeDerivative
    variable = wg
    v = etab0
    Fj_names = 'rhogm rhogb'
    hj_names = 'hm   hb'
    coupled_variables = 'etam0 etab0'
  [../]
  [gas_prod]
    type = MaskedBodyForce
    variable = wg
    mask = prod
    value = 0.2716
  []
[]

[Materials]
[prod]
  type = ParsedMaterial
  property_name = prod
  coupled_variables = 'etab0'
  material_property_names = 'h_psi F_dot'
  expression = '(etab0<0.5)*F_dot'
  outputs = exodus
  output_properties = prod
[]
[./const]
  type = GenericConstantMaterial
  prop_names =  'L  kappa_p    mu_p     gmb Vm        Va         Egf Evf   cveq_b cgeq_b T    kbz       F_dot  cv0   l_star t_star E_star'
  prop_values = '10 1.5e-10    3e10     1.5 2.462e22  4.090e-2    5.6 3.27  0.4425 0.5575 1200 8.617e-5  1e-8   0.005 1e-9   1      1e9'
[../]
[hvk_m]
  type = DerivativeParsedMaterial
  property_name = hvk_m
  material_property_names = 'hm(etam0,etab0) Va km'
  coupled_variables = 'etam0 etab0'
  expression = 'hm/Va/km'
[../]
[hvk_b]
  type = DerivativeParsedMaterial
  property_name = hvk_b
  material_property_names = 'hb(etam0,etab0) Va kb'
  coupled_variables = 'etam0 etab0'
  expression = 'hb/Va/kb'
[../]
# Switching Functions
  [./ha_test]
    type = SwitchingFunctionMultiPhaseMaterial
    h_name = hm
    all_etas = 'etam0 etab0'
    phase_etas = 'etam0'
  [../]
  [./hb_test]
    type = SwitchingFunctionMultiPhaseMaterial
    h_name = hb
    all_etas = 'etam0 etab0'
    phase_etas = 'etab0'
  [../]
 # Chemical Densities
  [./rhovm]
    type = DerivativeParsedMaterial
    coupled_variables = 'wv'
    property_name = rhovm
    material_property_names = 'Va km cveq_m'
    expression = 'wv/Va^2/km + cveq_m/Va'
    derivative_order = 2
    outputs = exodus
    output_properties = rhovm
  [../]
  [./rhovb]
    type = DerivativeParsedMaterial
    coupled_variables = 'wv'
    property_name = rhovb
    material_property_names = 'Va kb cveq_b'
    expression = 'wv/Va^2/kb + cveq_b/Va'
    derivative_order = 2
    outputs = exodus
    output_properties = rhovb
  [../]
  [./rhogm]
    type = DerivativeParsedMaterial
    coupled_variables = 'wg'
    property_name = rhogm
    material_property_names = 'Va km cgeq_m'
    expression = 'wg/Va^2/km + cgeq_m/Va'
    derivative_order = 2
    outputs = exodus
    output_properties = rhogm
  [../]
  [./rhogb]
    type = DerivativeParsedMaterial
    coupled_variables = 'wg'
    property_name = rhogb
    material_property_names = 'Va kb cgeq_b'
    expression = 'wg/Va^2/kb + cgeq_b/Va'
    derivative_order = 2
    outputs = exodus
    output_properties = rhogb
  [../]
  [./rhog]
    type = DerivativeParsedMaterial
    property_name = rhog
    material_property_names = 'rhogm(wg) rhogb(wg) hm(etam0,etab0) hb(etam0,etab0)'
    expression = 'rhogm*hm + rhogb*hb'
    derivative_order = 2
    outputs = exodus
    output_properties = rhog
  [../]
  [./rhov]
    type = DerivativeParsedMaterial
    property_name = rhov
    material_property_names = 'rhovm(wv) rhovb(wv) hm(etam0,etab0) hb(etam0,etab0)'
    expression = 'rhovm*hm + rhovb*hb'
    derivative_order = 2
    outputs = exodus
    output_properties = rhov
  [../]
  [./cv]
    type = DerivativeParsedMaterial
    material_property_names = 'Va rhovm(wv) rhovb(wv) hm(etam0,etab0) hb(etam0,etab0) l_star'
    expression = 'Va * (hm * rhovm + hb * rhovb)'
    property_name = cv
    outputs = exodus
    output_properties = cv
  [../]
  [./cg]
    type = DerivativeParsedMaterial
    material_property_names = 'Va rhogm(wg) rhogb(wg) hm(etam0,etab0) hb(etam0,etab0) l_star'
    expression = 'Va * (hm * rhogm + hb * rhogb)'
    property_name = cg
    outputs = exodus
    output_properties = cg
  [../]
  [./cveq_m]
    property_name = cveq_m
    type = DerivativeParsedMaterial
    material_property_names = 'Evf T kbz'
    expression = '1e-6'#'exp(-Evf/kbz/T)'
    outputs = exodus
    output_properties = cveq_m
  [../]
  [./cgeq_m]
    property_name = cgeq_m
    type = DerivativeParsedMaterial
    material_property_names = 'Egf T kbz'
    expression = '1e-6'#'exp(-Egf/kbz/T)'
    outputs = exodus
    output_properties = cgeq_m
  [../]

  [./kappa_nondim]
    type = DerivativeParsedMaterial
    property_name = kappa
    material_property_names = 'kappa_p l_star E_star'
    expression = 'kappa_p/E_star/(l_star^2)'
    derivative_order = 2
    outputs = exodus
    output_properties = kappa
  [../]
  [./mu_nondim]
    type = DerivativeParsedMaterial
    property_name = mu
    material_property_names = 'mu_p E_star'
    expression = 'mu_p/E_star'
    derivative_order = 2
    outputs = exodus
    output_properties = mu
  [../]
 # Diffusivity
  [./vMobility]
    type = DerivativeParsedMaterial
    property_name = Dchiv
    material_property_names = 'Dv chiv'
    expression = 'Dv*chiv'
    derivative_order = 2
    outputs = exodus
    output_properties = Dchiv
  [../]
  [./vchi] #Non-dim
    type = DerivativeParsedMaterial
    property_name = chiv
    material_property_names = 'Va hm(etam0,etab0) km hb(etam0,etab0) kb'
    expression = '(hm/km + hb/kb) / Va^2'
    derivative_order = 2
    outputs = exodus
    output_properties = chiv
  [../]
  [./Dv] #Non-dim
    type = DerivativeParsedMaterial
    property_name = Dv
    material_property_names = 'T kbz F_dot cv0'
    expression = 10#'8.24e-10*exp(-4.96/kbz/T) + 1.27e-27*exp(-1.73/kbz/T)*sqrt(F_dot) + 4.5e-44*F_dot'
    derivative_order = 2
    outputs = exodus
    output_properties = Dv
  [../]
  [./gMobility]
    type = DerivativeParsedMaterial
    property_name = Dchig
    material_property_names = 'Dg chig'
    expression = 'Dg*chig'
    derivative_order = 2
    outputs = exodus
    output_properties = Dchig
  [../]
  [./gchi] #Non-dim
    type = DerivativeParsedMaterial
    property_name = chig
    material_property_names = 'Va hm km hb kb'
    expression = '(hm/km + hb/kb) / Va^2'
    derivative_order = 2
    outputs = exodus
    output_properties = chig
  [../]
  [./Dg]
    type = DerivativeParsedMaterial
    property_name = Dg
    material_property_names = 'T kbz F_dot t_star l_star'
    expression = 10 #'(1.288e-3*exp(-5.842/kbz/T) + 2.2837e-79*exp(-(1.32271/kbz/T)^2 - (28.425/kbz/T))*sqrt(F_dot) + 2.623e-35*exp(-1.588/kbz/T)*F_dot)'
    derivative_order = 2
    outputs = exodus
    output_properties = Dg
  [../]
# Grand Potential Densities
  [./omegam]
    type = DerivativeParsedMaterial
    coupled_variables = 'wv wg'
    property_name = omegam
    material_property_names = 'Va km cveq_m cgeq_m'
    expression = '-0.5*wv^2/Va^2/km-wv/Va*cveq_m - 0.5*wg^2/Va^2/km-wg/Va*cgeq_m'
    derivative_order = 2
    outputs = exodus
    output_properties = omegam
  [../]
  [./km] #nondim
    type = ParsedMaterial
    property_name = km
    material_property_names = 'cv0 kbz Evf T E_star'
    constant_names =       'k1        k2'
    constant_expressions = '-1.787e6  1.287e10'
    expression = '100*((k1*T + k2)/(cv0 - exp(-Evf/kbz/T))) / E_star'
    outputs = exodus
    output_properties = km
  [../]
  [./omegab]
    type = DerivativeParsedMaterial
    coupled_variables = 'wv wg'
    property_name = omegab
    material_property_names = 'Va kb f_0'
    constant_names = 'cveq_b cgeq_b'
    constant_expressions = '0.5575 0.4425'
    expression = '-0.5*wv^2/Va^2/kb - wv/Va*cveq_b - 0.5*wg^2/Va^2/kb - wg/Va*cgeq_b + f_0'
    derivative_order = 2
    outputs = exodus
    output_properties = omegab
  [../]
  [./kb]
    type = ParsedMaterial
    property_name = kb
    material_property_names = 'T E_star'
    expression = '(2.320e8*T + 8.6e7)/E_star'
    outputs = exodus
    output_properties = kb
  [../]
  [./f_0]
    type = ParsedMaterial
    property_name = f_0
    material_property_names = 'T E_star'
    expression = '(-6.28976e3*T + 1.28704e8)/E_star'
  [../]
  [isbubble]
    type = ParsedMaterial
    property_name = isbubble
    coupled_variables = 'etab0'
    expression = '(etab0>0.5)'
  []
  [f_chem]
    type = ParsedMaterial
    property_name = f_chem
    material_property_names = 'hb hm omegab omegam'
    expression = 'hb*omegab + hm*omegam'
    outputs = exodus
    output_properties = f_chem
  []

  ##############################################################################
  # Nucleation Parameters                        ###############################
  ##############################################################################
  [G_nuc]
    type = DerivativeParsedMaterial
    property_name = G_nuc
    material_property_names = 'km cveq_m cgeq_m cgeq_b cv(wv) cg(wg)'
    expression = '(cv - cveq_m) + (cg-cgeq_m)'
    outputs = exodus
    output_properties = G_nuc
  []
  [nucl_matl]
    type = ParsedMaterial
    property_name = nucl_matl
    coupled_variables = 'etab0'
    material_property_names = 'cg(wg) T G_nuc'
    constant_names = 'k1star k2'
    constant_expressions = '1e-6 	0.005'
    expression = 'if(cg>1.5e-6,if(etab0<0.0001, k1star * exp(-k2/G_nuc/T), 0),0)'
    outputs = exodus
  []
[../]
[UserObjects]
  [B_inserter]
    type = DiscreteNucleationInserter
    hold_time = 1000
    probability = nucl_matl
    radius = 2
  []
  [B_map]
    type = DiscreteNucleationMap
    periodic = etab0
    inserter = B_inserter
    int_width = 0.5
  []
[]
[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Postprocessors]
  [./dt]
    type = TimestepSize
  [../]
  [./DOF]
    type = NumDOFs
  [../]
  [./rhov_total]
    type = ElementIntegralMaterialProperty
    mat_prop = rhov
    execute_on = 'initial timestep_end'
  [../]
  [./rhog_total]
    type = ElementIntegralMaterialProperty
    mat_prop = rhog
    execute_on = 'initial timestep_end'
  [../]
  [./cgtotal]
    type = ElementIntegralMaterialProperty
    mat_prop = cg
  [../]
  [./cvtotal]
    type = ElementIntegralMaterialProperty
    mat_prop = cv
  [../]
  [bubblearea]
    type = ElementIntegralMaterialProperty
    mat_prop = isbubble
  []
  [bubblepresence]
    type = ElementIntegralVariablePostprocessor
    variable = etab0
  []
  [num_bubbles]
    type = FeatureFloodCount
    variable = etab0
    threshold = 0.9
  []
  [nuc_rate]
    type = DiscreteNucleationData
    inserter = B_inserter
    value = RATE
  []
  [nuc_insertions]
    type = DiscreteNucleationData
    inserter = B_inserter
    value = INSERTIONS
  []
  [nuc_count]
    type = DiscreteNucleationData
    inserter = B_inserter
    value = COUNT
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'#'-pc_type -ksp_gmres_restart -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'lu       superlu_dist'#'asm      31                 lu           1'
  line_search = 'none'
  automatic_scaling = true
  nl_max_its = 15
  l_max_its = 50
  l_tol = 1e-3
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-12
  end_time = 1e7
  dtmin = 2e-20
  dtmax = 50
  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.1
    cutback_factor = 0.5
    growth_factor = 1.5
    optimal_iterations = 15
  [../]
[]

[Adaptivity]
  marker = combo
  max_h_level = 3
  [Indicators]
    [etaind]
      type = GradientJumpIndicator
      variable = etab0
    []
  []
  [Markers]
    [combo]
      type = ComboMarker
      markers = 'etaerror nuc'
    []
    [etaerror]
      type = ErrorFractionMarker
      coarsen = 0.2
      indicator = etaind
      refine = 0.7
    []
    [nuc]
      type = DiscreteNucleationMarker
      map = B_map
    []
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[Outputs]
  [./exodus]
    type = Exodus
    execute_on = 'initial timestep_end final'
    min_simulation_time_interval = 0
  [../]
  [./csv]
    type = CSV
  [../]
  print_linear_residuals = true
[]

GiudGiud · 2024-11-15T00:27:54Z

GiudGiud
Nov 15, 2024
Collaborator

Hello

The results seem very oscillatory. Is the mesh sufficiently discretized? Are the time steps small enough?
I think you are supposed to resolve the interface of the grain with the mesh size, which is not happening here. I think adaptivity is trying to address it but it's not clear to me that it is doing enough

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Discrete Nucleation Module inserting nuclei with negative value order parameter in grand potential model #29081

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment

{{title}}

Select a reply

Discrete Nucleation Module inserting nuclei with negative value order parameter in grand potential model #29081

aclin4 Nov 14, 2024

Replies: 1 comment

GiudGiud Nov 15, 2024 Collaborator

aclin4
Nov 14, 2024

GiudGiud
Nov 15, 2024
Collaborator