I am trying to solve non linear equation systems in gas dynamics using nlsolve, but the solver seems to breakdown when I call it at h = 0.1R_A2Theta_D (see code below)
using NLsolve
using DifferentialEquations
## About this code module
# Julia Program to solve Prandtl-Meyer flow considering high temperature effects
## GAS PROPERTIES
#comment out the appropriate section/lines to select
##Oxygen
R_A2 = 8.314*1000/32; #the molar gas constant for Oxygen
Θ_d = 59500; #dissociation temperature in Kelvin
ρ_d = 150000; #dissociation density
#flag = 1; #to display on the plot
##Nitrogen
#R_A2 = 8.314*1000/28; #the molar gas constant for Oxygen
#Θ_d = 113000; #dissociation temperature in Kelvin
#ρ_d = 130000; #dissociation density
#flag =2; #to display on the plot
##SOLUTION
## STATE EQUATION for high temperature gas dynamics in-terms of pressure and temperature input
#Intended to calculate inlet properties
function state(P,T)
#x[1]-->α_*, x[2]-->log10(ρ_d/ρ)
function g!(G,x)
G[1] = P - (1+x[1])*R_A2*T*ρ_d/10^x[2]
k = (10^(x[2]))*exp(-Θ_d/T)
G[2] = x[1] - 0.5*(sqrt(k^2 + 4k) -k)
end
res = nlsolve(g!,[0.5,6])
x = res.zero[1];#dissociation
y = res.zero[2];#density ratio
s = R_A2*((3*log(T/Θ_d) + (x*(1-2*log(x)))-((1-x)*log(1-x)) - (1+x)*log(1/10^y)))
a = sqrt(R_A2*(T/Θ_d)*Θ_d*(x*(1-x^2)*(1+2*(T/Θ_d)) + (8+3*x-x^2)*(T/Θ_d)^2)/(x*(1-x) + 3*(2-x)*(T/Θ_d)^2))
h = R_A2*((4+x)*T + x*Θ_d)
#comment out the following two lines to get absolute enthalpies
#h = h/(R_A2*Θ_d)
#s = s/R_A2
vec = [s,a,h,T/Θ_d,log10(R_A2*ρ_d*Θ_d/P),y,x]
return vec
end
##Mach function calculator for High-Temperature Gas Dynamics
#NON-LINEAR EQUATION SYSTEM FOR STATE
#x[1]-->α_*, x[2]-->T/Θ_d, x[3]-->log10(ρ_d/ρ)
function mach(h ,s)
s_c = s/R_A2#reduced entropy
h_c = h/(R_A2*Θ_d)#reduced enthalpy
function f!(F, x)
F[1] = h_c-((4+x[1])*x[2]+x[1])
F[2] = exp(x[1]-s_c) - (x[1]^(2*x[1]))*((1-x[1])^(1-x[1]))*((1/10^x[3])^(1+x[1]))/(x[2]^3)
k = ((10^x[3]))*exp(-1/x[2])
F[3] = x[1] - 0.5*(sqrt(k^2+4k)-k)
end
results = nlsolve(f!, [ 0.5, 0.01,4.5])
x = results.zero[1]#dissociation
y = results.zero[2]#temperature ratio
z = ρ_d/(10^results.zero[3])#density
a = sqrt(R_A2*y*Θ_d*(x*(1-x^2)*(1+2*y) + (8+3*x-x^3)*y^2)/(x*(1-x) + 3*(2-x)*y^2))
vec = [x,y,z,a]
return vec
end
## Main program for Prandtl-Meyer Expansion
##Upstream properties
p_a = 121590;#upstream pressure in pascal
T_a = 6140;#upstream temperature in kelvin
M_a = 1;#upstream Mach number
state_info = state(p_a,T_a);
s_a = state_info[1];#inlet entropy
u_a = M_a*state_info[2];#inlet velocity
h_a = state_info[3];#inlet enthalpy
h_t_a = h_a + 0.5*u_a^2;#inlet stagnaion enthalpy
##Prandtl-Meyer Solver
res = mach(0.1*R_A2*Θ_d,s_a)
println(" ");
println(res);
I get the following error:
LoadError: DomainError with -6.055454134369954e-6:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
in expression starting at C:\Users\SURYA\Desktop\IIST Internship\Programs\Julia programs\MachNumber.jl:99
throw_exp_domainerror(::Float64) at math.jl:37
^ at math.jl:872 [inlined]
(::var"#f!#66"{Float64,Float64})(::Array{Float64,1}, ::Array{Float64,1}) at MachNumber.jl:62
finite_difference_jacobian!(::Array{Float64,2}, ::var"#f!#66"{Float64,Float64}, ::Array{Float64,1}, ::FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:central},Float64}, ::Nothing; relstep::Float64, absstep::Float64, colorvec::UnitRange{Int64}, sparsity::Nothing, dir::Bool) at jacobians.jl:370
finite_difference_jacobian!(::Array{Float64,2}, ::Function, ::Array{Float64,1}, ::FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:central},Float64}, ::Nothing) at jacobians.jl:295
finite_difference_jacobian! at jacobians.jl:295 [inlined]
fj_finitediff! at oncedifferentiable.jl:139 [inlined]
(::NLSolversBase.var"#j_finitediff!#22"{Array{Float64,1},NLSolversBase.var"#fj_finitediff!#21"{var"#f!#66"{Float64,Float64},FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:central},Float64}}})(::Array{Float64,2}, ::Array{Float64,1}) at oncedifferentiable.jl:144
jacobian!!(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,2}, ::Array{Float64,1}) at interface.jl:142
jacobian!! at interface.jl:140 [inlined]
jacobian! at interface.jl:135 [inlined]
trust_region_(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::NLsolve.NewtonTrustRegionCache{Array{Float64,1}}) at trust_region.jl:181
trust_region at trust_region.jl:229 [inlined]
trust_region at trust_region.jl:229 [inlined]
nlsolve(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64) at nlsolve.jl:26
nlsolve at nlsolve.jl:18 [inlined]
nlsolve(::Function, ::Array{Float64,1}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at nlsolve.jl:52
nlsolve(::Function, ::Array{Float64,1}) at nlsolve.jl:46
mach(::Float64, ::Float64) at MachNumber.jl:66
top-level scope at MachNumber.jl:99