Specific Heat Ratio Does Not Accept Input from UDF
Posted: 12 Sep 2020, 01:56
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:
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:
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:
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.
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