Specific Heat Ratio Does Not Accept Input from UDF

Clearly defined bug reports and their fixes
Post Reply
gschrank
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

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  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
End

Constants
  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
End

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

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
End

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

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
End

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

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

Boundary Condition 3
  Target Boundaries(1) = 2 
  Name = "Outlet"
  External Pressure = -1e5
  Velocity 2 = 0
End
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
    IMPLICIT NONE
    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
    !------------------------------------------------------------------------------

    Materials=>GetMaterial()

    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

    heatcapratio=(5.0D0/2.0D0*(xemassfrac+hemassfrac)+7.0D0/2.0D0*n2massfrac)/&
        (3.0D0/2.0D0*(xemassfrac+hemassfrac)+5.0D0/2.0D0*n2massfrac)

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.
Attachments
OPUtil.F90
(31.25 KiB) Downloaded 245 times
raback
Site Admin
Posts: 4801
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Specific Heat Ratio Does Not Accept Input from UDF

Post by raback »

Hi

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.

-Peter
Post Reply