Specific Heat Ratio Does Not Accept Input from UDF

Clearly defined bug reports and their fixes
Post Reply
Posts: 17
Joined: 21 Aug 2016, 05:40
Antispam: Yes

Specific Heat Ratio Does Not Accept Input from UDF

Post by gschrank »

I am using Elmer 8.4 on Ubuntu 18.04 LTS. Elmer is complied using gfortran.

I have written a User Defined Function to calculate the specific heat ratio for a multi-component system for use with the Navier-Stokes "Flow" Solver. The solver will run just fine when I put a number for Specific Heat Ratio. However, when I attempt to call the UDF, the Navier-Stokes Solver returns the error:

Code: Select all

ERROR:: ComputeChange: Numerical Error: Norm of solution appears to be NaN
This happens even when the value that the UDF returns and the value put manually put in are the same. I know that the the UDF returns a value because I can get it to be called by "SaveData" "SaveMaterials" solver.

The sif I am using is:

Code: Select all

  Mesh DB "." "."
  Include Path ""
  Results Directory ""

  Max Output Level = 5
  Coordinate System = Cartesian 2D
  Simulation Type = Steady State
  Steady State Max Iterations = 1
  Output Intervals = 1
  Timestepping Method = BDF
  BDF Order = 1
  Solver Input File = case.sif
  Post File = case.vtu

  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19

Body 1
  Target Bodies(1) = 1
  Name = "Body 1"
  Equation = 1
  Material = 1

Solver 1
  Equation = Navier-Stokes
  Variable = Flow Solution[Velocity:2 Pressure:1]
  Procedure = "FlowSolve" "FlowSolver"
  Exec Solver = Always
  Stabilize = False
  Bubbles = True
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-5
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Direct
  Linear System Direct Method = Umfpack

Equation 1
  Name = "Equation 1"
  NS Convect = False
  Active Solvers(1) = 1

Material 1
  Name = "Air (room temperature)"
  Reference Temperature = 300
  Viscosity = 1.983e-5
  Heat expansion Coefficient = 3.43e-3
  Compressibility Model = Perfect Gas
  Reference Pressure = 0

  xe fraction = Real 0.25
  n2 fraction = Real 0.25
  he fraction = Real 0.5

  Specific Heat Ratio = Variable Concentration, Pressure, Temperature; Real Procedure "OPUtil" "calculateheatcapratio"
  Heat Conductivity = 0.0257
  Relative Permittivity = 1.00059
  Sound speed = 343.0
  Heat Capacity = 1005.0
  Density = 1.205

Boundary Condition 1
  Target Boundaries(2) = 3 4 
  Name = "Walls"
  Noslip wall BC = True

Boundary Condition 2
  Target Boundaries(1) = 1 
  Name = "Inlet"
  Velocity 1 = 1
  Velocity 2 = 0

Boundary Condition 3
  Target Boundaries(1) = 2 
  Name = "Outlet"
  External Pressure = -1e5
  Velocity 2 = 0
For these inputs, the UDF outputs a value of 1.6. The solver runs when that value is directly put in for Specific Heat Ratio.

The code for the UDF is:

Code: Select all

FUNCTION calculateheatcapratio(Model,n,arguments)RESULT(heatcapratio)
    !Defines the Heat Capacity ratio as a function of gas fraction, Note* Xe and He (gamma = 1.666) are assumed to be perfect
    !monotonic gasses and N2 is assumed to be a perfect diatomic gas (gamma = 1.4).
    USE DefUtils
    TYPE(Model_t) :: Model
    INTEGER :: n
    REAL(KIND=dp) :: arguments(3)
    REAL(KIND=dp) :: heatcapratio
    TYPE(ValueList_t), POINTER :: Materials
    REAL(KIND=dp) :: n2_fraction, xe_fraction, he_fraction
    REAL(KIND=dp) :: avgmolarmass, xemassfrac, n2massfrac, hemassfrac
    LOGICAL :: found


    xe_fraction=GetConstReal(Materials, 'xe fraction', found)
    CALL FoundCheck(found, 'xe fraction', 'fatal')

    n2_fraction = GetConstReal(Materials, 'n2 fraction', found)
    CALL FoundCheck(found, 'n2 fraction' , 'fatal')

    he_fraction = GetConstReal(Materials, 'he fraction', found)
    CALL FoundCheck(found, 'he fraction', 'fatal')

    !Call fatal if the gas fractions don't add to 1

    IF (ABS(1-he_fraction-n2_fraction-xe_fraction)>1e-5) THEN
        CALL Fatal('GetSpinDestructionRate', &
            'Gas fractions do not add to 1')
    END IF

    avgmolarmass = xe_fraction*131.293D0+n2_fraction*28.0134D0+he_fraction*4.002602D0

    xemassfrac = xe_fraction*131.293D0/avgmolarmass
    n2massfrac = n2_fraction*28.0134D0/avgmolarmass
    hemassfrac = he_fraction*4.0026202D0/avgmolarmass


END FUNCTION calculateheatcapratio
I've attached the source file for the UDF. There are several functions in there that are used to calculate various quantities, and the "calculateheatcapratio" is the only one that appears to not be working with the Elmer infrastructure.
(31.25 KiB) Downloaded 98 times
Site Admin
Posts: 4173
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland

Re: Specific Heat Ratio Does Not Accept Input from UDF

Post by raback »


Whereas almost all material parameters in Elmer may depend on other variables the specific heat ratio cannot. Technically it depends whether the code uses ListGetConstReal to fetch just one value or ListGetReal to get values for all elemental nodes.

So unfortunately some coding would be needed. I'll add some error detection for this to the ListGet routines. Altering the solver module would be a somewhat bigger task.

Post Reply