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`

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
```

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
```