Discrepancy in analytic solution

[makeclean]

Hi All

Got a weird one here, sure its a simple one but I can't seem to figure it. Im trying to reproduce a validation from Code_Aster EDF's thermo-mechanical code, which can easily be viewed here - https://www.simscale.com/docs/validation-cases/hollow-sphere-convection-radiation/

I reproduced the geometry in Trelis, applied the same boundary conditions, and setup the problem, correctly as far as I can tell.

However, running the simulation with 1st and 2nd order hexes, I see a consistent difference between the analytic solution from - https://www.code-aster.org/doc/default/fr/man_v/v4/v4.43.001.pdf - of about 3 degrees on the hot surface and 2 degrees on the cooler surface. Ive done this with a tabualated and explicit function to describe the heat flux on the radiative surface but give very similar results.

Any tips/suggestions appreciated.

reactor.zip

Thanks

Andy

[hugary1995]

Your radiation BC is wrong. I modified your input and it works:

[Mesh]
  type = FileMesh
  file = '../mesh/reactor.e'
[]

[Variables]
  [temp]
  []
[]

[AuxVariables]
  [temp_in_C]
  []
  [radiation_flux]
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
  []
[]

[AuxKernels]
  [K_to_C]
    type = ParsedAux
    variable = 'temp_in_C'
    function = 'temp-273.15'
    args = 'temp'
  []
  [radiation]
    type = ParsedAux
    variable = radiation_flux
    args = 'temp'
    function = 'sigma*epsilon*(start_t^4 - temp^4)'
    constant_names = 'sigma epsilon start_t'
    constant_expressions = '5.73e-8 0.6 773.15'
  []
[]

[BCs]
  [abcd]
    type = CoupledConvectiveHeatFluxBC
    boundary = 'abcd'
    variable = temp
    htc = '133.5'
    T_infinity = 293.15
  []
  [efgh]
    type = CoupledNeumannBC
    variable = temp
    boundary = 'efgh'
    coupled_var = 'radiation_flux'
  []
[]

[Materials]
  [density]
    type = GenericConstantMaterial
    prop_names = 'density  thermal_conductivity cp'
    prop_values = '1 40 1'
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  automatic_scaling = true
  type = Steady
  petsc_options = '-snes_converged_reason'
  petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap -ksp_gmres_restart'
  petsc_options_value = 'asm lu 1 101'
  line_search = none
[]

[Postprocessors]
  [external_surf_ave]
    type = SideAverageValue
    variable = temp_in_C
    boundary = 'abcd'
  []
  [internal_surf_ave]
    type = SideAverageValue
    variable = temp_in_C
    boundary = 'efgh'
  []
[]

[Outputs]
  exodus = true
[]

[makeclean]

Thanks @hugary1995 thats what I needed! I need to make one change to your input, for whatever reason I dont see that BC type being available, neither does the documentation. I had to use

    type = CoupledVarNeumannBC
    v = 'radiation_flux'

But other than that, can confirm that this gives what I was after! Superb!

Only one question really, what was incorrect about my original functionneumanbc, shouldn't it behave the same way? What am I missing?

[hugary1995]

Oh my bad, that CoupledNeumannBC is something from my app, sorry for the confusion.

I'm on my phone so cannot confirm this, but I remember you originally used sigma epsilon (T0^4 - x^4), but the equation should really be sigma epsilon (T0^4 - T^4). This is a robin BC strictly speaking.