--GAS TURBINE TRANSIENT SIMULATION in BeamNG
--v0.2 14/7/16 
--Author: Davded http://www.beamng.com/members/davded.188/
--With thanks to Max98, Fabian Enkler, Dummiesman

--TODO:
--auto/calibration of turbine/massflow/area properties at init 			DONE
--fix or add minimum throttle/engine dying/engine start conditions				Needs to be done in seperate iterative calibration run
--relate throttle to fuel flow, rather then C/C temperature ratio		DONE
--ISA Atmosphere model															Ingame atmosphere model?
--mass flow rate should be independent of RPM									Probably not nessesary
--3D compressor/turbine map interpolation										"
--above two points will enable compressor surge/stall & choke
--multiple engine/independent throttle support									Fairly simple, run loop for n. engines, index variable like temperature[2][4], or [24]

--add other parts of engine
	--afterburner														EASY	Just add q between stations 4 & 5
	--water/methanol injection
	--sonic exhuast flow												DONE
	--Shock diamonds from over/under expansion									Only really for visuals
	--variable geometry nozzles											EASY	Already halfway with sonic exhuast flow
	--bypass flow
	--fan simulation
	--multiple spools
	--supersonic intakes												EASY 	Normal Shock calc's
	--variable geometry/oblique shock diffusion									Hard if multiple shock reflections											
	--bleed air, muliple exhuast points, eg Harrier Jump Jet

local M = {}

--[[
--engine properties for Whittle Turbojet/RR Welland
--engine 'design' points describe engine conditions at maximum throttle, 100% RPM and usually sea level ISA conditions
local designCPR = 4.4           	--Max Compressor Pressure ratio
local designRPM = 12000     		--Max engine RPM
local idleRPM = 2000				--engine start/idle RPM
local designTemperature = 650+273   --K Turbine entry temperature
local designThrust = 7200 			--N Engine Thrust

local inertia = 4  					--kg/m^2 total engine inertia
local angularFriction = 20    		--Nm   angular drag torque

local calorificFuel = 46000000		--J/kg -calorific value of the fuel used (kerosene) 46000kJ/kg
--actual mass flow rate: 13kg/s, simulation gives 14kg/s
--]]

--Input properties for RR RB44 Tay
local designCPR = 4.9        		--Max Compressor Pressure ratio
local designRPM = 11000     		--Max engine RPM
local idleRPM = 2000				--engine start/idle RPM
local designTemperature = 1033   	--K Turbine entry temperature
local designThrust = 34000 			--N Engine Thrust

local inertia = 14  				--kg/m^2 total engine inertia
local angularFriction = 60    		--Nm   angular drag torque

local calorificFuel = 46000000		--J/kg -calorific value of the fuel used (kerosene) 46000kJ/kg
--actual mass flow rate: 60kg/s, this simulation gives 60.4kg/s

local idleThrottle = 0.02			--2% thottle, well be calculated on loading later
local designTTR = 0
local designmdot = 0 
local designmdotFuel = 0

--turbomachinery (adiabatic) efficiencies (if 0 or undefined assume ideal of 1, check that 0 < n < 1 in code later)
local aCompressor = 0.8			
local aTurbine = 0.8 

--isentropic efficiencies	ideal = 1
local nCompressor = 0.9			
local nTurbine = 0.8 
--assume no losses in C/C
local nIntake = 0.85					
local nNozzle = 0.9

--variables:
local input_throttle = 0
local input_airspeed = 0

--ambient conditions
local Ta = 293			--K ambient temperature (overwritten by obj:getEnvTemperature() later)
local Pa = 100000		--Pa ambient pressure

--velocities
local Va = 0			--m/s Airspeed
local Vj = 0			--m/s jet exit velocity

--Mach numbers
local Ma = 0			
local Mj = 0

--pressures
local pressure = {}

--temperatures
local temperature = {}

--enthalpys
local enthalpy = {}

--gas data
local RAir = 287    		--specific gas constant for air kJ kg−1 K−1
local rhoAir = 1.225 		--air density 

--transient properties initalisation
local throttle = idleThrottle	--starts engine running at idle conditions
local RPM = idleRPM				--engine RPM
local RPMratio = 0				--RPM ratio to design point rpm
local mdot = 0					--kg/s mass flow rate
local mdotFuel = 0				--kg/s fuel flow
local cWork = 0					--Compressor work
local tWork = 0					--Turbine work
local excessWork = 0			--Power excess or defecit
local ts = 0					--specific thrust N/kg/s

local i_dt = 0


local function Cp_get(temperature)
	--specific heat capacity for constant pressure
	return ((0.00019*temperature+0.94)*1000)
end

local function a_get(g,R,T)
	--speed of sound
	return math.sqrt(g*R*T)
end

local function onReset()
	log('I', "wellandTurbojet.lua", "RESET")
	RPM = idleRPM
	throttle = idleThrottle
	--Ta = obj:getEnvTemperature() + 273
end

local function onExtensionLoaded()

	--Calibration at design point, 100% throttle, design RPM, design Thurst
	temperature[0] = Ta
	pressure[0] = Pa
	local Cp = Cp_get(temperature[0])
	local gamma = Cp/(Cp-RAir)
	
	Ma = 0		--static thurst design point -- intake has no effect, p1 = pa, t1 = ta
	mdot = 10	--[kg/s] initial assumption for mass flow rate
	
	local mdotError = 1
	
	--iteratation to find correct Turbine Temperature Ratio (TTR) & massFlow (mdot)
	while mdotError > 0.001 do
	
		--frictional drag at designRPM
		--rotation friction is expressed as a torque
		--all other works are specific (mass flow independent) [m^2/s^2] but we need the actual work
		local designAV = designRPM*0.1047		--convert to rad/s
		local fPower = angularFriction*designAV --[kg*m^2/s^3] or [W] (SI base units)
		local fWork = fPower/mdot				--[m^2/s^2] changed to specific work
		
		
		
		--Intake/Diffuser (0) -> (1)
		--static conditions T.F.
		temperature[1] = temperature[0]
		pressure[1] = pressure[0]
		
		
		
		--Compressor (1) -> (2)
		local CPR = designCPR			--CPR = designCP
		pressure[2] = pressure[1]*CPR	--P2
		
		cWork = Cp*temperature[1]*(CPR^((gamma-1)/gamma)-1)/aCompressor --doesn't 'know' about T2 yet, T.F. have to apply efficiency term again
		
		local CTR = (CPR)^((gamma-1)/gamma)
		local flowCTR = 1 + nCompressor*(CTR - 1)
		temperature[2] = temperature[1]*flowCTR
		
		--Combustion Chamber (2) -> (3)
		--No pressure change
		--100% thottle, max TIT used
		pressure[3] = pressure[2]
		temperature[3] = designTemperature --Heat added
		enthalpy[2] = temperature[2]*Cp_get(temperature[2])
		enthalpy[3] = temperature[3]*Cp_get(temperature[3])
		
		--mdotFuel is added to mdot
		local qCombustion = enthalpy[3] - enthalpy[2]		--heat added is difference in enthalpy
		mdotFuel = mdot*qCombustion/(calorificFuel)	--calculate fuel flow required
		
		local mdotBurnt = mdot + mdotFuel
		
		--Turbine (3) -> (4)
		--turbine work = compressor work + frictional drag
		--get temperature drop therefore designTTR
		local CpBurn = Cp_get(temperature[3])
		local gammaBurn = CpBurn/(CpBurn-RAir)
		
		--power extracted by turbine
		tWork = (fWork + cWork)*(mdot/(mdotBurnt))	--add mass flow of burnt fuel to engine flow
		local TTR = (1-tWork/(aTurbine*CpBurn*temperature[3]))
		temperature[4] = temperature[3]*TTR
		
		local flowTTR = 1+nTurbine*(TTR-1)
		local TPR = (flowTTR)^(gammaBurn/(gammaBurn-1))
		pressure[4] = pressure[3]*TPR
		
		
		
		--Nozzle (4) -> (5)
		--assume expansion to ambient pressure (Pa)
		--Get gas properties
		local CpExhaust = Cp_get(temperature[4])
		local gammaExhaust = CpExhaust/(CpExhaust-RAir)
		enthalpy[4] = temperature[4]*CpExhaust
		
		--Check if flow is sonic in nozzle
		local rcrit = (1 / (1 - (1/nNozzle)*((gammaExhaust-1)/(gammaExhaust+1))) )^(gammaExhaust/(gammaExhaust-1))
		
		if rcrit > pressure[4]/Pa then
		
			pressure[5] = Pa --able to expand to ambient pressure (idealised thrust)
			local NPR = pressure[5]/pressure[4]
			local NTR = (NPR)^((gammaExhaust-1)/gammaExhaust)
			temperature[5] = temperature[4]*(1+nNozzle*(NTR-1))
			enthalpy[5] = temperature[5]*Cp_get(temperature[5])
			
			if enthalpy[5] < enthalpy[4] then	--condition to deal with low/minimum power settings
				Vj = math.sqrt(2*(enthalpy[4]-enthalpy[5]))
			else
				Vj = 0
			end
			
			ts = Vj*(mdotBurnt/mdot)-Va --specific thrust N/kg^s (no pressure contribution, P5 = Pa)
		else
			--sonic flow, can't expand through acceleration to ambient pressure without flow going supersonic, T.F. need pressure contribution
			--efficiency term is included in rcrit
			pressure[5] = pressure[4]/rcrit
			pressure[50] = pressure[5]*(1+((gammaExhaust-1)/2))^(gammaExhaust/(gammaExhaust-1))

			temperature[50] = temperature[4] 								--T[50] is hacky stagnation temp, no heat added
			temperature[5] = temperature[50]/(1+((gammaExhaust-1)/2))

			Vj = a_get(gammaExhaust,RAir,temperature[5])   					--speed of sound at exhaust exit, since M = 1, Vj = a

			local MFPcrit = (math.sqrt(gammaExhaust/RAir)) / ((1 + ((gammaExhaust-1)/2)) ^ ((gammaExhaust+1) / (2*(gammaExhaust-1))))

			local As = (math.sqrt(temperature[50])/(pressure[50]))*(1/MFPcrit) 	--specific area (can be used to determine nozzle throat area)
			ts = (Vj*(mdotBurnt/mdot)-Va)+As*(pressure[5]-Pa)*(mdotBurnt/mdot) 	--specific thrust with contribution from pressure
		end
		
		local mdotNew = designThrust/ts
		mdotError = math.abs(mdotNew-mdot)
		mdot = mdotNew
		
		--output design values
		designTTR = TTR
		designmdot = mdot
		designmdotFuel = mdotFuel
		
	end
	log('I', "wellandTurbojet.lua", "CALIBRATION COMPLETE")
	log('I', "wellandTurbojet.lua", "Max Flow Rate = " .. tostring(designmdot) .. " kg/s")
	log('I', "wellandTurbojet.lua", "Max Fuel Flow = " .. tostring(designmdotFuel) .. " kg/s")
	log('I', "wellandTurbojet.lua", "TTR Maximium  = " .. tostring(designTTR))				
	
	log('I', "wellandTurbojet.lua", "LOADED")
end

local function updateGFX(dt)

	Va = obj:getVelocity():length()

	if dt ~= nil then
		i_dt = dt
		if input.keys["Y"] then
				throttle = throttle + (0.5 * dt)
		end

		if input.keys["H"] then
				throttle = throttle - (0.5 * dt)
		end
	end
	
	RPMratio = RPM/designRPM
	mdot = designmdot*RPMratio			--transient mass flow rate kg/s
	
	--NOMENCLATURE:
	--Split thermodynamics into stages:
	--Diffuser - Compressor - C/C - Turbine - Nozzle
	--0--------1------------2-----3---------4-------5
	--a----I----------C--------B-------T--------N---j
	
	--Intake/Diffuser (0) -> (1)
	--Free Stream -> Diffuser Duct -> Compressor Front
	
	--ITR - Intake Temperature Ratio
	--CCR - Compressor Compression Ratio... etc.
	
	temperature[0] = Ta
	pressure[0] = Pa
	
	local Cp = Cp_get(temperature[0])
	local gamma = Cp/(Cp-RAir)
	
	Ma = Va/a_get(gamma,RAir,temperature[0])	--Mach No.
	
		--temperatures
	local T1ideal = temperature[0]*(1+(gamma-1)/2*Ma*Ma) 	--idealised T1 (insentropic process of a perfect gas)
	local ITR = 1 + nIntake*(T1ideal/temperature[0] - 1)			--isentropic efficiency applied to idealised T1
	temperature[1] = temperature[0]*ITR
	
		--pressures
	local IPR = (ITR)^(gamma/(gamma-1))
	pressure[1] = pressure[0]*IPR
	
	
	
	--Compressor (1) -> (2)
	--Compression ratio is proportional to RPM
	
	local CPR = 1+RPMratio*(designCPR-1)	--linear compressor map
	pressure[2] = pressure[1]*CPR			--P2
	
	cWork = Cp*temperature[1]*(CPR^((gamma-1)/gamma)-1)/aCompressor --doesn't 'know' about T2 yet, T.F. have to apply efficiency term again
	
	local CTR = (CPR)^((gamma-1)/gamma)
	local flowCTR = 1 + nCompressor*(CTR - 1)
	temperature[2] = temperature[1]*flowCTR
	
	
	
	--Combustion Chamber (2) -> (3)
	--heat addition proportional to throttle setting
	--if 100% throttle, 100% exit temp is used
	--if 0%, exit temp (temperature[3]) ~= entry temp (temperature[2])
	--idle will be at some minimum throttle %, so temperature[3] > temperature[2]
	--pressure kept roughly the same in combustion chamber so
	pressure[3] = pressure[2]
	
	throttle = math.max(idleThrottle, math.min(1, throttle))	--0 <= throttle <= 1
	
	local addedFuel	= throttle*(designmdotFuel)		--throttle proportional to fuel flow
	local addedEnthalpy = addedFuel*calorificFuel/mdot
	local exitEnthalpy = addedEnthalpy + enthalpy[2]
	local maxExitEnthalpy = designTemperature*Cp_get(designTemperature)
	
	if exitEnthalpy > maxExitEnthalpy then
		--prevents exceeding the max turbine inlet tempera
		enthalpy[3] = maxExitEnthalpy
		temperature[3] = designTemperature
	else
		enthalpy[3] = exitEnthalpy
		temperature[3] = enthalpy[3]/Cp_get(designTemperature)
	end
		
	local qCombustion = enthalpy[3] - enthalpy[2]				
	mdotFuel = mdot*qCombustion/(calorificFuel)
	local mdotBurnt = mdot + mdotFuel

	--Turbine   (3) -> (4)
	--calculating how much energy the turbine needs to extract from flow to power compressor
	--tWork = enthalpy[3]-enthalpy[4]
	local CpBurn = Cp_get(temperature[3])
	local gammaBurn = CpBurn/(CpBurn-RAir)
	
	local TTR = (1-(RPMratio^0.2)*(1-designTTR)) 	--turbine map
	if RPM < idleRPM*0.75 then
		TTR = 1
	end
	
	tWork = aTurbine*CpBurn*temperature[3]*(1-TTR)	--work extracted by turbine
	temperature[4] = temperature[3]*TTR
	enthalpy[4] = temperature[4]*Cp_get(temperature[4])

	local flowTTR = 1+nTurbine*(TTR-1)
	local TPR = (flowTTR)^(gammaBurn/(gammaBurn-1))
	pressure[4] = pressure[3]*TPR
	
	
	
	--Nozzle (4) -> (5)
	--assume expansion to ambient pressure (Pa)
	--Get gas properties
	local CpExhaust = Cp_get(temperature[4])
	local gammaExhaust = CpExhaust/(CpExhaust-RAir)
	enthalpy[4] = temperature[4]*CpExhaust
	
	--Check if flow is sonic in nozzle
	local rcrit = (1 / (1 - (1/nNozzle)*((gammaExhaust-1)/(gammaExhaust+1))) )^(gammaExhaust/(gammaExhaust-1))
	
	if rcrit > pressure[4]/Pa then
		--exhuast flow is < Mach 1
		pressure[5] = Pa --able to expand to ambient pressure (idealised thrust)
		local NPR = pressure[5]/pressure[4]
		local NTR = (NPR)^((gammaExhaust-1)/gammaExhaust)
		temperature[5] = temperature[4]*(1+nNozzle*(NTR-1))
		enthalpy[5] = temperature[5]*Cp_get(temperature[5])
		
		if enthalpy[5] < enthalpy[4] then	--condition to deal with low/minimum power settings
			Vj = math.sqrt(2*(enthalpy[4]-enthalpy[5]))
		else
			Vj = 0
		end
		
		ts = Vj*(mdotBurnt/mdot)-Va --specific thrust N/kg^s (no pressure contribution, P5 = Pa)
	else
		--exhuast gas is choked at Mach 1
		pressure[5] = pressure[4]/rcrit
		pressure[50] = pressure[5]*(1+((gammaExhaust-1)/2))^(gammaExhaust/(gammaExhaust-1))

		temperature[50] = temperature[4] 								--Because No heat added, T[50] is hacky stagnation temp
		temperature[5] = temperature[50]/(1+((gammaExhaust-1)/2))

		Vj = a_get(gammaExhaust,RAir,temperature[5])   					--speed of sound at exhaust exit, since M = 1, Vj = a

		local MFPcrit = (math.sqrt(gammaExhaust/RAir)) / ((1 + ((gammaExhaust-1)/2)) ^ ((gammaExhaust+1) / (2*(gammaExhaust-1))))
		local As = (math.sqrt(temperature[50])/(pressure[50]))*(1/MFPcrit) 	--specific area
		ts = (Vj*(mdotBurnt/mdot)-Va)+As*(pressure[5]-Pa)*(mdotBurnt/mdot)	--specific thrust with contribution from pressure
	end
	
	
	
	--Transient variables
	local angularVelocity = RPM*0.1047				--conversion ro rad/s
	
	local fPower = angularFriction*angularVelocity 	--[kg*m^2/s^3] or [W] (SI base units)
	local fWork = fPower/mdot						--specific work
	excessWork = tWork*(mdotBurnt/mdot)-cWork-fWork					--all work is specific work, W/kg/s. Actual work = specificWork*mdot
	
	local angularAccel = excessWork*mdot/(angularVelocity*inertia)	--excess work is used to accelerate the engine
	angularVelocity = angularVelocity + i_dt*angularAccel		--new angular velocity after timestep
	
	RPM = (angularVelocity/0.1047)
	RPM = math.max(RPM, 0)	--prevent -ve.
   
	local finalThrust = math.max(mdot*ts, 0)
	electrics.values['turbineThrust'] = finalThrust
	
	gui.send('gasTurbineData', {
					--essential
					throttle = (throttle * 100),
					fuelflow = mdotFuel*1000,
					massflow = mdot,
					rpm = RPM,
					thrust = finalThrust,
					
					--optional
					p0 = pressure[0]/100000,	--Conversion from Pa to bar
					p1 = pressure[1]/100000,
					p2 = pressure[2]/100000,
					p3 = pressure[3]/100000,
					p4 = pressure[4]/100000,
					p5 = pressure[5]/100000,
					t0 = temperature[0],
					t1 = temperature[1],
					t2 = temperature[2],
					t3 = temperature[3],
					t4 = temperature[4],
					t5 = temperature[5],
					Wturbine = tWork*mdotBurnt/1000,
					Wcompressor = cWork*mdot/1000,
					Wexcess = excessWork*mdot/1000,
					Tspecific = ts,
					Ua = Va,
					Uj = Vj
			})
end

M.update = update
M.updateGFX = updateGFX
M.onExtensionLoaded = onExtensionLoaded
M.onReset = onReset
return M