diff --git a/docs/src/sizing/sizing.md b/docs/src/sizing/sizing.md index 9e787a1a..66047621 100644 --- a/docs/src/sizing/sizing.md +++ b/docs/src/sizing/sizing.md @@ -20,7 +20,7 @@ TASOPT.woper(pari, parg, parm, para, pare, A sized aircraft's mission performance can be obtained (`mission!`), along with operation constraints via a pitch trim calculation (`balance`) and balanced field length calculation (`takeoff!`). ```@docs -TASOPT.mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1) +TASOPT.mission!(pari, parg, parm, para, pare, Ldebug) TASOPT.takeoff!(pari, parg, parm, para, pare, initeng, diff --git a/src/mission/mission.jl b/src/mission/mission.jl index 9aefd8a2..d32d2adc 100644 --- a/src/mission/mission.jl +++ b/src/mission/mission.jl @@ -1,5 +1,5 @@ """ - mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1) + mission!(pari, parg, parm, para, pare, Ldebug) Runs aircraft through mission, calculating fuel burn and other mission variables. @@ -30,20 +30,12 @@ NOTE: In an upcoming revision, an `aircraft` struct and auxiliary indices will be passed in lieu of pre-sliced `par` arrays. """ -function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, iairf, initeng, ipc1) +function mission!(pari, parg, parm, para, pare, Ldebug)#, iairf, initeng, ipc1) t_prop = 0.0 calc_ipc1 = true ifirst = true - if pari[iiengmodel] == 0 - # Drela engine model - use_NPSS = false - else - # NPSS - use_NPSS = true - end - # HACK TODO add the para back # iairf initeng = 0 @@ -59,19 +51,6 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i ifclose = pari[iifclose] ifuel = pari[iifuel] - if use_NPSS - mofWpay = parg[igmofWpay] - mofWMTO = parg[igmofWMTO] - PofWpay = parg[igPofWpay] - PofWMTO = parg[igPofWMTO] - - - # BLI - fBLIf = parg[igfBLIf] - DAfsurf = para[iaDAfsurf, ipcruise1] - KAfTE = para[iaKAfTE, ipcruise1] - end - # Mission range Rangetot = parm[imRange] para[iaRange, ipdescentn] = Rangetot @@ -335,16 +314,6 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i FoW = zeros(Float64, iptotal) FFC = zeros(Float64, iptotal) Vgi = zeros(Float64, iptotal) - if use_NPSS - dPdt = zeros(Float64, iptotal) - dygdt = zeros(Float64, iptotal) - - - pare[iePLH2, ipclimb1] = 1.5 - pare[ieyg, ipclimb1] = 0.1 - end - - # integrate trajectory over climb @inbounds for ip = ipclimb1:ipclimbn @@ -363,29 +332,18 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i Wpay = parg[igWpay] neng = parg[igneng] - if use_NPSS - - if pari[iiengtype] == 1 - initfanPCT = 80.0 - finalfanPCT = parg[igfanPCT] - frac = float(ip - ipclimb1) / float(ipclimbn - ipclimb1) - fanPCT = initfanPCT * (1.0 - frac) + finalfanPCT * frac + dgamV = 1.0 + Ftotal = 0.0 + DoL = 0.0 + mdotf = 0.0 + V = 0.0 - mofft = (mofWpay * Wpay + mofWMTO * W) / neng - Pofft = (PofWpay * Wpay + PofWMTO * W) / neng - end - - dgamV = 1.0 - Ftotal = 0.0 - DoL = 0.0 - mdotf = 0.0 - V = 0.0 - if (Ldebug) - printstyled(@sprintf("\t%5s %10s %10s %10s %10s %10s %10s %10s \n", - "iterg", "dgamV", "gamV", "cosg", "BW", "Ftotal", "DoL", "V"); color=:light_green) - end + if (Ldebug) + printstyled(@sprintf("\t%5s %10s %10s %10s %10s %10s %10s %10s \n", + "iterg", "dgamV", "gamV", "cosg", "BW", "Ftotal", "DoL", "V"); color=:light_green) end + @inbounds for iterg = 1:itergmax V = sqrt(2.0 * BW * cosg / (ρ * S * CL)) Mach = V / Vsound @@ -409,52 +367,15 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i end cdsum!(pari, parg, view(para, :, ip), view(pare, :, ip), icdfun) + icall = 1 + icool = 1 - if use_NPSS - ρ0 = pare[ierho0, ip] - u0 = pare[ieu0, ip] - Φinl = 0.5 * ρ0 * u0^3 * (DAfsurf * fBLIf) / 2.0 - Kinl = 0.5 * ρ0 * u0^3 * (KAfTE * fBLIf) / 2.0 # Assume 2 engines - - if pari[iiengtype] == 0 - NPSS_success, Ftotal, η, P, Hrej, heatexcess, - mdotf, deNOx, EINOx1, EINOx2, FAR, Tt3, OPR, Wc3, Tt41, EGT = NPSS_TEsysOD(NPSS, para[iaalt, ip], Mach, 0.0, pare[ieTt4, ip], Kinl, Φinl, 0.0, 0.0, ifirst, parg, parpt, pare, ip) - - pare[iedeNOx, ip] = deNOx - pare[ieEINOx1, ip] = EINOx1 - pare[ieEINOx2, ip] = EINOx2 - pare[ieemot:ieethermal, ip] .= η - pare[ieHrejmot:ieHrejtot, ip] .= Hrej - pare[ieHexcess, ip] = heatexcess - else - NPSS_success, Ftotal, heatexcess, - mdotf, EINOx1, FAR, Tt3, OPR, Wc3, Tt41, EGT = NPSS_TFsysOD(NPSS, para[iaalt, ip], Mach, 0.0, pare[ieTt4, ip], ifirst, parg, parpt, pare, ip, fanPCT, mofft, Pofft) - - pare[ieEINOx1, ip] = EINOx1 - - pare[ieemot:ieethermal, ip] .= 0.0 - pare[ieHrejmot:ieHrejtot, ip] .= 0.0 - pare[ieHexcess, ip] = heatexcess - end - - ifirst = false - pare[ieOPR, ip] = OPR - pare[ieTt3, ip] = Tt3 - pare[ieWc3, ip] = Wc3 - pare[iemdotf, ip] = mdotf - - DoL = para[iaCD, ip] / para[iaCL, ip] - else - icall = 1 - icool = 1 - - ichoke5, ichoke7 = tfcalc!(pari, parg, view(para, :, ip), view(pare, :, ip), ip, icall, icool, initeng) + ichoke5, ichoke7 = tfcalc!(pari, parg, view(para, :, ip), view(pare, :, ip), ip, icall, icool, initeng) - Ftotal = pare[ieFe, ip] * parg[igneng] - TSFC = pare[ieTSFC, ip] - DoL = para[iaCD, ip] / para[iaCL, ip] + Ftotal = pare[ieFe, ip] * parg[igneng] + TSFC = pare[ieTSFC, ip] + DoL = para[iaCD, ip] / para[iaCL, ip] - end # Calculate improved flight angle ϕ = Ftotal / BW @@ -484,20 +405,9 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i pare[ieFe, ip] = Ftotal # Store integrands for range and weight integration using a predictor-corrector scheme FoW[ip] = Ftotal / (BW * cosg) - DoL - if (!use_NPSS) - FFC[ip] = Ftotal * TSFC / (W * V * cosg) - else - FFC[ip] = mdotf * gee / (W * V * cosg) - end - Vgi[ip] = 1.0 / (V * cosg) + FFC[ip] = Ftotal * TSFC / (W * V * cosg) - if (use_NPSS) - ΔT = pare[ieT0, ip] - 20.0 - P = pare[iePLH2, ip] - yg = pare[ieyg, ip] - dygdt[ip] = mdotf / ρmix(yg, P) - dPdt[ip] = dPdt_LH2(ΔT, mdotf, P, yg) - end + Vgi[ip] = 1.0 / (V * cosg) Mach = para[iaMach, ip] CL = para[iaCL, ip] @@ -513,29 +423,14 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i FFCavg = 0.5 * (FFC[ip] + FFC[ip-1]) Vgiavg = 0.5 * (Vgi[ip] + Vgi[ip-1]) - if (use_NPSS) - dPdtavg = 0.5 * (dPdt[ip] + dPdt[ip-1]) - dygdtavg = 0.5 * (dygdt[ip] + dygdt[ip-1]) - end - dR = (dh + 0.5 * dVsq / gee) / FoWavg dt = dR * Vgiavg rW = exp(-dR * FFCavg) #Ratio of weights Wᵢ₊₁/Wᵢ - if (use_NPSS) - dP = dPdtavg * Vgiavg * dR - dyg = dygdtavg * Vgiavg * dR - end - para[iaRange, ip] = para[iaRange, ip-1] + dR para[iatime, ip] = para[iatime, ip-1] + dt para[iafracW, ip] = para[iafracW, ip-1] * rW - if (use_NPSS) - pare[iePLH2, ip] = pare[iePLH2, ip-1] + dP - pare[ieyg, ip] = pare[ieyg, ip-1] + dyg - end - ifirst = false end @@ -562,19 +457,11 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i dt = dR * Vgi[ip] rW = exp(-dR * FFC[ip]) - if (use_NPSS) - dP = dPdt[ip] * Vgi[ip] * dR - dyg = dygdt[ip] * Vgi[ip] * dR - end para[iaRange, ip+1] = para[iaRange, ip] + dR para[iatime, ip+1] = para[iatime, ip] + dt para[iafracW, ip+1] = para[iafracW, ip] * rW - if (use_NPSS) - pare[iePLH2, ip+1] = pare[iePLH2, ip] + dP - pare[ieyg, ip+1] = pare[ieyg, ip] + dyg - end end end # done integrating climb @@ -609,68 +496,14 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i BW = W + para[iaWbuoy, ip] F = BW * (DoL + para[iagamV, ip]) Wpay = parg[igWpay] + icall = 2 + icool = 1 + ichoke5, ichoke7 = tfcalc!(pari, parg, view(para, :, ip), view(pare, :, ip), ip, icall, icool, initeng) - if use_NPSS - if pari[iiengtype] == 1 - mofWpay = parg[igmofWpay] - mofWMTO = parg[igmofWMTO] - PofWpay = parg[igPofWpay] - PofWMTO = parg[igPofWMTO] - neng = parg[igneng] - - mofft = (mofWpay * Wpay + mofWMTO * W) / neng - Pofft = (PofWpay * Wpay + PofWMTO * W) / neng - end - - Mach = pare[ieM0, ip] - # Run powertrain [TODO] actually should run this to a required thrust not Tt4 - ρ0 = pare[ierho0, ip] - u0 = pare[ieu0, ip] - Φinl = 0.5 * ρ0 * u0^3 * (DAfsurf * fBLIf) / 2.0 - Kinl = 0.5 * ρ0 * u0^3 * (KAfTE * fBLIf) / 2.0 # Assume 2 engines - - if pari[iiengtype] == 0 - NPSS_success, Ftotal, η, P, Hrej, heatexcess, - mdotf, deNOx, EINOx1, EINOx2, FAR, Tt3, OPR, Wc3, Tt41, EGT = NPSS_TEsysOD(NPSS, para[iaalt, ip], Mach, F, 0 * pare[ieTt4, ip], Kinl, Φinl, 0.0, 0.0, ifirst, parg, parpt, pare, ip) - - pare[ieOPR, ip] = OPR - pare[ieTt3, ip] = Tt3 - pare[ieWc3, ip] = Wc3 - pare[iedeNOx, ip] = deNOx - pare[ieEINOx1, ip] = EINOx1 - pare[ieEINOx2, ip] = EINOx2 - pare[iemdotf, ip] = mdotf - pare[ieemot:ieethermal, ip] .= η - pare[ieHrejmot:ieHrejtot, ip] .= Hrej - pare[ieHexcess, ip] = heatexcess - pare[ieFe, ip] = Ftotal - - else - NPSS_success, Ftotal, heatexcess, - mdotf, EINOx1, FAR, Tt3, OPR, Wc3, Tt41, EGT = NPSS_TFsysOD(NPSS, para[iaalt, ip], Mach, F, pare[ieTt4, ip], ifirst, parg, parpt, pare, ip, 0.0, mofft, Pofft) - - ifirst = false - - pare[ieOPR, ip] = OPR - pare[ieTt3, ip] = Tt3 - pare[ieWc3, ip] = Wc3 - pare[ieEINOx1, ip] = EINOx1 - pare[iemdotf, ip] = mdotf - pare[ieemot:ieethermal, ip] .= 0.0 - pare[ieHrejmot:ieHrejtot, ip] .= 0.0 - pare[ieHexcess, ip] = heatexcess - end - else - icall = 2 - icool = 1 - ichoke5, ichoke7 = tfcalc!(pari, parg, view(para, :, ip), view(pare, :, ip), ip, icall, icool, initeng) - end end # set cruise-climb climb angle, from fuel burn rate and atmospheric dp/dz - if (!use_NPSS) - TSFC = pare[ieTSFC, ip] - end + TSFC = pare[ieTSFC, ip] V = pare[ieu0, ip] p0 = pare[iep0, ip] ρ0 = pare[ierho0, ip] @@ -685,35 +518,18 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i # ∴γ ≈ ̇mf*p/(ρWV) # For conventional aircraft can represent this as function of TSFC: (see TASOPT docs Eq 438) # gamVcr1 = DoL*p0*TSFC/(ρ0*gee*V - p0*TSFC) - if use_NPSS - gamVcr1 = mdotf * p0 / (ρ0 * W * V) # TSFC not the most useful for turbo-electric systems. Buoyancy weight has been neglected. - else - gamVcr1 = DoL * p0 * TSFC / (ρ0 * gee * V - p0 * TSFC) - end + gamVcr1 = DoL * p0 * TSFC / (ρ0 * gee * V - p0 * TSFC) + para[iagamV, ip] = gamVcr1 - if use_NPSS - Ftotal = pare[ieFe, ip] - else - Ftotal = pare[ieFe, ip] * parg[igneng] - end + Ftotal = pare[ieFe, ip] * parg[igneng] + cosg = cos(gamVcr1) FoW[ip] = Ftotal / (BW * cosg) - DoL - if use_NPSS - FFC[ip] = mdotf * gee / (W * V * cosg) - else - FFC[ip] = Ftotal * TSFC / (W * V * cosg) - end - Vgi[ip] = 1.0 / (V * cosg) + FFC[ip] = Ftotal * TSFC / (W * V * cosg) - if use_NPSS - ΔT = pare[ieT0, ip] - 20.0 - P = pare[iePLH2, ip] - yg = pare[ieyg, ip] - dygdt[ip] = mdotf / ρmix(yg, P) - dPdt[ip] = dPdt_LH2(ΔT, mdotf, P, yg) - end + Vgi[ip] = 1.0 / (V * cosg) #---- set end-of-cruise point "d" using cruise and descent angles, #- and remaining range legs @@ -759,83 +575,28 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i W = para[iafracW, ip] * WMTO BW = W + para[iaWbuoy, ip] Ftotal = BW * (DoL + para[iagamV, ip]) - if !use_NPSS - pare[ieFe, ip] = Ftotal / parg[igneng] - end - - if use_NPSS - Mach = pare[ieM0, ip] - # Run powertrain [TODO] actually should run this to a required thrust not Tt4 - ρ0 = pare[ierho0, ip] - u0 = pare[ieu0, ip] - Φinl = 0.5 * ρ0 * u0^3 * (DAfsurf * fBLIf) / 2.0 - Kinl = 0.5 * ρ0 * u0^3 * (KAfTE * fBLIf) / 2.0 # Assume 2 engines - - if pari[iiengtype] == 0 - NPSS_success, Ftotal, η, P, Hrej, heatexcess, - mdotf, deNOx, EINOx1, EINOx2, FAR, Tt3, OPR, Wc3, Tt41, EGT = NPSS_TEsysOD(NPSS, para[iaalt, ip], Mach, F, pare[ieTt4, ip], Kinl, Φinl, 0.0, 0.0, ifirst, parg, parpt, pare, ip) - - pare[ieOPR, ip] = OPR - pare[ieTt3, ip] = Tt3 - pare[ieWc3, ip] = Wc3 - pare[iedeNOx, ip] = deNOx - pare[ieEINOx1, ip] = EINOx1 - pare[ieEINOx2, ip] = EINOx2 - pare[iemdotf, ip] = mdotf - pare[ieemot:ieethermal, ip] .= η - pare[ieHrejmot:ieHrejtot, ip] .= Hrej - pare[ieHexcess, ip] = heatexcess - else - NPSS_success, Ftotal, heatexcess, - mdotf, EINOx1, FAR, Tt3, OPR, Wc3, Tt41, EGT = NPSS_TFsysOD(NPSS, para[iaalt, ip], Mach, F, pare[ieTt4, ip], ifirst, parg, parpt, pare, ip, 0.0, mofft, Pofft) - - pare[ieOPR, ip] = OPR - pare[ieTt3, ip] = Tt3 - pare[ieWc3, ip] = Wc3 - pare[ieEINOx1, ip] = EINOx1 - pare[iemdotf, ip] = mdotf - pare[ieemot:ieethermal, ip] .= 0.0 - pare[ieHrejmot:ieHrejtot, ip] .= 0.0 - pare[ieHexcess, ip] = heatexcess - end - else - icall = 2 - icool = 1 + pare[ieFe, ip] = Ftotal / parg[igneng] - ichoke5, ichoke7 = tfcalc!(pari, parg, view(para, :, ip), view(pare, :, ip), ip, icall, icool, initeng) + icall = 2 + icool = 1 - end + ichoke5, ichoke7 = tfcalc!(pari, parg, view(para, :, ip), view(pare, :, ip), ip, icall, icool, initeng) + TSFC = pare[ieTSFC, ip] - if (!use_NPSS) - TSFC = pare[ieTSFC, ip] - end V = pare[ieu0, ip] p0 = pare[iep0, ip] ρ0 = pare[ierho0, ip] DoL = para[iaCD, ip] / para[iaCL, ip] - if use_NPSS - gamVcr2 = mdotf * p0 / (ρ0 * W * V) - else - gamVcr2 = DoL * p0 * TSFC / (ρ0 * gee * V - p0 * TSFC) - end + + gamVcr2 = DoL * p0 * TSFC / (ρ0 * gee * V - p0 * TSFC) para[iagamV, ip] = gamVcr2 cosg = cos(gamVcr1) FoW[ip] = Ftotal / (BW * cosg) - DoL - if (!use_NPSS) - FFC[ip] = Ftotal * TSFC / (W * V * cosg) - else - FFC[ip] = mdotf * gee / (W * V * cosg) - end + FFC[ip] = Ftotal * TSFC / (W * V * cosg) Vgi[ip] = 1.0 / (V * cosg) - if use_NPSS - ΔT = pare[ieT0, ip] - 20.0 - P = pare[iePLH2, ip] - dygdt[ip] = mdotf / ρmix(yg, P) - dPdt[ip] = dPdt_LH2(ΔT, mdotf, P, yg) - end ip1 = ipcruise1 ipn = ipcruisen @@ -844,12 +605,6 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i FFCavg = 0.5 * (FFC[ipn] + FFC[ip1]) Vgiavg = 0.5 * (Vgi[ipn] + Vgi[ip1]) - if use_NPSS - dPdtavg = 0.5 * (dPdt[ipn] + dPdt[ip1]) - dygdtavg = 0.5 * (dygdt[ipn] + dygdt[ip1]) - dP = dPdtavg * Vgiavg * dRcruise - dyg = dygdtavg * Vgiavg * dRcruise - end dtcruise = dRcruise * Vgiavg rWcruise = exp(-dRcruise * FFCavg) @@ -858,10 +613,6 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i para[iatime, ipn] = para[iatime, ip1] + dtcruise para[iafracW, ipn] = para[iafracW, ip1] * rWcruise - if use_NPSS - pare[iePLH2, ipn] = pare[iePLH2, ip1] + dP - pare[ieyg, ipn] = pare[ieyg, ip1] + dyg - end # set intermediate points over cruise, if any, just by interpolating for ip = ipcruise1+1:ipcruisen-1 @@ -889,8 +640,6 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i end # Descent - #TODO descent has not been properly implemented yet - below it just performs an approximate scaling based - # on TASOPT to estimate the descent fuel burn ip = ipdescent1 pare[iep0, ip] = pare[iep0, ipcruisen] pare[ieT0, ip] = pare[ieT0, ipcruisen] @@ -909,205 +658,183 @@ function mission!(pari, parg, parm, para, pare, Ldebug, NPSS_PT, NPSS, ipc1)#, i para[iafracW, ip] = para[iafracW, ipcruisen] para[iaWbuoy, ip] = para[iaWbuoy, ipcruisen] + Rd = para[iaRange, ipdescent1] + Re = para[iaRange, ipdescentn] + altd = para[iaalt, ipdescent1] + alte = para[iaalt, ipdescentn] + para[iagamV, ipdescent1] = gamVde1 + para[iagamV, ipdescentn] = gamVden - if use_NPSS - # Mission fuel fractions and weights - fracWa = para[iafracW, ipclimb1] - # fracWe = para[iafracW, ipdescentn] - fracWe = para[iafracW, ipdescent1] * (1 - (1 - 0.993) * ((para[iaalt, ipdescent1] / ft_to_m) / 39858.0))# Pp temp set to 95% of ipdescent1 instead of ipdescentn since descent has not been calculated yet - freserve = parg[igfreserve] - fburn = fracWa - fracWe # Burnt fuel fraction - ffuel = fburn * (1.0 + freserve) - Wfuel = WMTO * ffuel - WTO = Wzero + Wfuel - - parm[imWTO] = WTO - parm[imWfuel] = Wfuel - # printstyled("Wfuel = $Wfuel \n fburn = $fburn \n", color=:red) - - Wburn = WMTO * fburn - parg[igWfburn] = Wburn - parm[imPFEI] = Wburn / gee * parg[igLHVfuel] * 1e6 / (parm[imWpay] * parm[imRange]) - end + for ip = ipdescent1+1:ipdescentn-1 + frac = float(ip - ipdescent1) / float(ipdescentn - ipdescent1) + R = Rd * (1.0 - frac) + Re * frac - if (!use_NPSS) - Rd = para[iaRange, ipdescent1] - Re = para[iaRange, ipdescentn] - altd = para[iaalt, ipdescent1] - alte = para[iaalt, ipdescentn] - para[iagamV, ipdescent1] = gamVde1 - para[iagamV, ipdescentn] = gamVden - - for ip = ipdescent1+1:ipdescentn-1 - frac = float(ip - ipdescent1) / float(ipdescentn - ipdescent1) - R = Rd * (1.0 - frac) + Re * frac - - alt = altd + (Re - Rd) * (gamVde1 * (frac - 0.5 * frac^2) + gamVden * 0.5 * frac^2) - gamVde = gamVde1 * (1.0 - frac) + gamVden * frac - para[iagamV, ip] = gamVde - - altkm = alt / 1000.0 - T0, p0, ρ0, a0, μ0 = atmos(altkm) - pare[iep0, ip] = p0 - pare[ieT0, ip] = T0 - pare[iea0, ip] = a0 - pare[ierho0, ip] = ρ0 - pare[iemu0, ip] = μ0 - - para[iaRange, ip] = R - para[iaalt, ip] = alt - - rhocab = max(parg[igpcabin], p0) / (RSL * TSL) - para[iaWbuoy, ip] = (rhocab - rho0) * gee * parg[igcabVol] - end - para[iaWbuoy, ipdescentn] = 0.0 + alt = altd + (Re - Rd) * (gamVde1 * (frac - 0.5 * frac^2) + gamVden * 0.5 * frac^2) + gamVde = gamVde1 * (1.0 - frac) + gamVden * frac + para[iagamV, ip] = gamVde - # integrate time and weight over descent - for ip = ipdescent1:ipdescentn + altkm = alt / 1000.0 + T0, p0, ρ0, a0, μ0 = atmos(altkm) + pare[iep0, ip] = p0 + pare[ieT0, ip] = T0 + pare[iea0, ip] = a0 + pare[ierho0, ip] = ρ0 + pare[iemu0, ip] = μ0 - # velocity calculation from CL, Weight, altitude - gamVde = para[iagamV, ip] - cosg = cos(gamVde) - W = para[iafracW, ip] * WMTO - BW = W + para[iaWbuoy, ip] - CL = para[iaCL, ip] - rho = pare[ierho0, ip] - V = sqrt(2.0 * BW * cosg / (rho * S * CL)) - Mach = V / pare[iea0, ip] + para[iaRange, ip] = R + para[iaalt, ip] = alt - para[iaMach, ip] = Mach - para[iaReunit, ip] = V * rho / pare[iemu0, ip] + rhocab = max(parg[igpcabin], p0) / (RSL * TSL) + para[iaWbuoy, ip] = (rhocab - rho0) * gee * parg[igcabVol] + end + para[iaWbuoy, ipdescentn] = 0.0 - pare[ieu0, ip] = V - pare[ieM0, ip] = Mach + # integrate time and weight over descent + for ip = ipdescent1:ipdescentn - # set pitch trim by adjusting CLh - Wf = W - Wzero - rfuel = Wf / parg[igWfuel] - itrim = 1 - balance(pari, parg, view(para, :, ip), rfuel, rpay, ξpay, itrim) + # velocity calculation from CL, Weight, altitude + gamVde = para[iagamV, ip] + cosg = cos(gamVde) + W = para[iafracW, ip] * WMTO + BW = W + para[iaWbuoy, ip] + CL = para[iaCL, ip] + rho = pare[ierho0, ip] + V = sqrt(2.0 * BW * cosg / (rho * S * CL)) + Mach = V / pare[iea0, ip] - if (ip == ipdescentn) - # use explicitly specified wing cdf,cdp - icdfun = 0 - else - # use airfoil database for wing cdf,cdp - icdfun = 1 - end - cdsum!(pari, parg, view(para, :, ip), view(pare, :, ip), icdfun) + para[iaMach, ip] = Mach + para[iaReunit, ip] = V * rho / pare[iemu0, ip] - # set up for engine calculation - sing = sin(gamVde) - cosg = cos(gamVde) - DoL = para[iaCD, ip] / para[iaCL, ip] - Fspec = BW * (sing + cosg * DoL) - pare[ieFe, ip] = Fspec / parg[igneng] + pare[ieu0, ip] = V + pare[ieM0, ip] = Mach - if (initeng == 0) - pare[iembf, ip] = pare[iembf, ip-1] - pare[iemblc, ip] = pare[iemblc, ip-1] - pare[iembhc, ip] = pare[iembhc, ip-1] - pare[iepif, ip] = pare[iepif, ip-1] - pare[iepilc, ip] = pare[iepilc, ip-1] - pare[iepihc, ip] = pare[iepihc, ip-1] - inite = 1 + # set pitch trim by adjusting CLh + Wf = W - Wzero + rfuel = Wf / parg[igWfuel] + itrim = 1 + balance(pari, parg, view(para, :, ip), rfuel, rpay, ξpay, itrim) + if (ip == ipdescentn) + # use explicitly specified wing cdf,cdp + icdfun = 0 + else + # use airfoil database for wing cdf,cdp + icdfun = 1 + end + cdsum!(pari, parg, view(para, :, ip), view(pare, :, ip), icdfun) - # make better estimate for new Tt4, adjusted for new ambient T0 - dTburn = pare[ieTt4, ip-1] - pare[ieTt3, ip-1] - OTR = pare[ieTt3, ip-1] / pare[ieTt2, ip-1] - Tt3 = pare[ieT0, ip] * OTR - pare[ieTt4, ip] = Tt3 + dTburn + 50.0 + # set up for engine calculation + sing = sin(gamVde) + cosg = cos(gamVde) + DoL = para[iaCD, ip] / para[iaCL, ip] + Fspec = BW * (sing + cosg * DoL) + pare[ieFe, ip] = Fspec / parg[igneng] - # make better estimate for new pt5, adjusted for new ambient p0 - pare[iept5, ip] = pare[iept5, ip-1] * pare[iep0, ip] / pare[iep0, ip-1] + if (initeng == 0) + pare[iembf, ip] = pare[iembf, ip-1] + pare[iemblc, ip] = pare[iemblc, ip-1] + pare[iembhc, ip] = pare[iembhc, ip-1] + pare[iepif, ip] = pare[iepif, ip-1] + pare[iepilc, ip] = pare[iepilc, ip-1] + pare[iepihc, ip] = pare[iepihc, ip-1] + inite = 1 - else - inite = initeng - end - # use fixed engine geometry, specified Fe - icall = 2 - # use previously-set turbine cooling mass flow - icool = 1 + # make better estimate for new Tt4, adjusted for new ambient T0 + dTburn = pare[ieTt4, ip-1] - pare[ieTt3, ip-1] + OTR = pare[ieTt3, ip-1] / pare[ieTt2, ip-1] + Tt3 = pare[ieT0, ip] * OTR + pare[ieTt4, ip] = Tt3 + dTburn + 50.0 - ichoke5, ichoke7 = tfcalc!(pari, parg, view(para, :, ip), view(pare, :, ip), ip, icall, icool, inite) + # make better estimate for new pt5, adjusted for new ambient p0 + pare[iept5, ip] = pare[iept5, ip-1] * pare[iep0, ip] / pare[iep0, ip-1] - # store effective thrust, effective TSFC - F = pare[ieFe, ip] * parg[igneng] - TSFC = pare[ieTSFC, ip] + else + inite = initeng + end - # store integrands for Range and Weight integration - FoW[ip] = F / (BW * cosg) - DoL - FFC[ip] = F / (W * V * cosg) * TSFC - Vgi[ip] = 1.0 / (V * cosg) + # use fixed engine geometry, specified Fe + icall = 2 + # use previously-set turbine cooling mass flow + icool = 1 - # if F < 0, then TSFC is not valid, so calculate mdot_fuel directly - mfuel = pare[ieff, ip] * pare[iemcore, ip] * parg[igneng] - FFC[ip] = gee * mfuel / (W * cosg * V) + ichoke5, ichoke7 = tfcalc!(pari, parg, view(para, :, ip), view(pare, :, ip), ip, icall, icool, inite) + # store effective thrust, effective TSFC + F = pare[ieFe, ip] * parg[igneng] + TSFC = pare[ieTSFC, ip] - if (ip > ipdescent1) - # corrector integration step, approximate trapezoidal - dh = para[iaalt, ip] - para[iaalt, ip-1] - dVsq = pare[ieu0, ip]^2 - pare[ieu0, ip-1]^2 + # store integrands for Range and Weight integration + FoW[ip] = F / (BW * cosg) - DoL + FFC[ip] = F / (W * V * cosg) * TSFC + Vgi[ip] = 1.0 / (V * cosg) - FoWavg = 0.5 * (FoW[ip] + FoW[ip-1]) - FFCavg = 0.5 * (FFC[ip] + FFC[ip-1]) - Vgiavg = 0.5 * (Vgi[ip] + Vgi[ip-1]) + # if F < 0, then TSFC is not valid, so calculate mdot_fuel directly + mfuel = pare[ieff, ip] * pare[iemcore, ip] * parg[igneng] + FFC[ip] = gee * mfuel / (W * cosg * V) - dR = para[iaRange, ip] - para[iaRange, ip-1] - dt = dR * Vgiavg - rW = exp(-dR * FFCavg) - para[iatime, ip] = para[iatime, ip-1] + dt - para[iafracW, ip] = para[iafracW, ip-1] * rW - end + if (ip > ipdescent1) + # corrector integration step, approximate trapezoidal + dh = para[iaalt, ip] - para[iaalt, ip-1] + dVsq = pare[ieu0, ip]^2 - pare[ieu0, ip-1]^2 - if (ip < ipdescentn) + FoWavg = 0.5 * (FoW[ip] + FoW[ip-1]) + FFCavg = 0.5 * (FFC[ip] + FFC[ip-1]) + Vgiavg = 0.5 * (Vgi[ip] + Vgi[ip-1]) - # predictor integration step, forward Euler - gamVde = para[iagamV, ip+1] - cosg = cos(gamVde) - W = para[iafracW, ip+1] * WMTO + dR = para[iaRange, ip] - para[iaRange, ip-1] + dt = dR * Vgiavg + rW = exp(-dR * FFCavg) - BW = W + para[iaWbuoy, ip+1] + para[iatime, ip] = para[iatime, ip-1] + dt + para[iafracW, ip] = para[iafracW, ip-1] * rW + end - CL = para[iaCL, ip+1] - rho = pare[ierho0, ip+1] - V = sqrt(2 * BW * cosg / (rho * S * CL)) - pare[ieu0, ip+1] = V + if (ip < ipdescentn) - dh = para[iaalt, ip+1] - para[iaalt, ip] - dVsq = pare[ieu0, ip+1]^2 - pare[ieu0, ip]^2 + # predictor integration step, forward Euler + gamVde = para[iagamV, ip+1] + cosg = cos(gamVde) + W = para[iafracW, ip+1] * WMTO - dR = para[iaRange, ip] - para[iaRange, ip-1] - dt = dR * Vgi[ip] - rW = exp(-dR * FFC[ip]) + BW = W + para[iaWbuoy, ip+1] - para[iatime, ip+1] = para[iatime, ip] + dt - para[iafracW, ip+1] = para[iafracW, ip] * rW + CL = para[iaCL, ip+1] + rho = pare[ierho0, ip+1] + V = sqrt(2 * BW * cosg / (rho * S * CL)) + pare[ieu0, ip+1] = V + + dh = para[iaalt, ip+1] - para[iaalt, ip] + dVsq = pare[ieu0, ip+1]^2 - pare[ieu0, ip]^2 + + dR = para[iaRange, ip] - para[iaRange, ip-1] + dt = dR * Vgi[ip] + rW = exp(-dR * FFC[ip]) + + para[iatime, ip+1] = para[iatime, ip] + dt + para[iafracW, ip+1] = para[iafracW, ip] * rW - end end + end - # mission fuel fractions and weights - fracWa = para[iafracW, ipclimb1] - fracWe = para[iafracW, ipdescentn] - freserve = parg[igfreserve] - fburn = fracWa - fracWe - ffuel = fburn * (1.0 + freserve) - Wfuel = WMTO * ffuel - WTO = Wzero + Wfuel + # mission fuel fractions and weights + fracWa = para[iafracW, ipclimb1] + fracWe = para[iafracW, ipdescentn] + freserve = parg[igfreserve] + fburn = fracWa - fracWe + ffuel = fburn * (1.0 + freserve) + Wfuel = WMTO * ffuel + WTO = Wzero + Wfuel - parm[imWTO] = WTO - parm[imWfuel] = Wfuel + parm[imWTO] = WTO + parm[imWfuel] = Wfuel + + # mission PFEI + Wburn = WMTO * fburn + parm[imPFEI] = Wburn/gee * pare[iehfuel, ipcruise1] / (parm[imWpay] * parm[imRange]) - # mission PFEI - Wburn = WMTO * fburn - parm[imPFEI] = Wburn/gee * pare[iehfuel, ipcruise1] / (parm[imWpay] * parm[imRange]) - end return t_prop end diff --git a/src/propsys/PMSM.jl b/src/propsys/PMSM.jl deleted file mode 100644 index b0e85319..00000000 --- a/src/propsys/PMSM.jl +++ /dev/null @@ -1,502 +0,0 @@ -""" -PMSM sizes the electric machine for the motor and generator - of the turbo electric powertrain - -Inputs: - - P : power [W] - - N : Rotational speed [1/s] - - ratSplit : dRot/dStator - - σAg : Air-gap shear stress [N/m²] - - -Outputs: - - Preq : Required input power - - η : Overall efficiency P/Preq - - W : Weight of machine [N] - - l : Length of machine [m] - - d : Diameter of machine [m] - - -Originally based of code by W. Enders -Modified and updated by P. Prashanth -""" -function PMSM(P::Float64, ratAsp::Float64, σAg::Float64, ratSplit::Float64, parte::Array{Float64, 1} ) - - # Setup/ assumptions - Vtip = 200 # [m/s] Rotor tip speed with retaining sleeve - - # Rotor dimensions (See Hendershot p89) - # |T| = |F*r| = (π*D*L)*σ * D/2 - # P = Tω = (π*D*L)*σ * D/2 *V/(D/2) - lRot = sqrt(P/(π*ratAsp*σAg*Vtip)) - dRot = lRot*ratAsp - - parte[ite_lRot] = lRot - # parte[ite_dRot] = dRot - - N = Vtip/(dRot*π) - - # [TODO] From Wilhelm's aircraft sizing code... but need to check why 240 - p = parte[ite_p] - - κM = 240/(2p) - - # Rotational speed - ω = 2π*N - - # Unpack paramters - ratAg = parte[ite_ratAg] - ratM = parte[ite_ratM ] - hRS = parte[ite_hRS ] - ratSp = parte[ite_ratSp] - ratSM = parte[ite_ratSM] - ratSd = parte[ite_ratSd] - ratW = parte[ite_ratW ] - Nshrt = parte[ite_Nshrt] - kpf = parte[ite_kpf ] - - - z = parte[ite_z] - Br = parte[ite_Br] - - μ0 = parte[ite_mu0] - - ρMag = parte[ite_rhoMag ] - ρCu = parte[ite_rhoCu ] - ρFe = parte[ite_rhoFe ] - ρSteel = parte[ite_rhoSteel] - - ratShft = parte[ite_ratShft] - τMax = parte[ite_tauMax ] - - ψ = parte[ite_psi] - # ------------------ - # Machine geometry - # ------------------ - Q = P/ω - Qmax = 1.5*Q - - rRot = dRot/2.0 - - hAg = ratAg*dRot + hRS # Air-gap Thickness - parte[ite_hAg] = hAg - - hM = ratM*hAg - parte[ite_hM] = hM - - # Electric frequency - rRoti = dRot/2 - (hM + hRS) #inner radius of rotor - parte[ite_rRoti] = rRoti - - f = p * N - Ω = f * 2π - Vtip ≤ (rRoti + 0.5*hM)*ω ? println("Warning Vtip test not passed") : - - # Calculate Slot geometry - NS = ratSp * p # Number of slots from slots-polepair ratio ratSp - hS = ratSM * hM # Slot height - hSd = ratSd * hS # Slot depression height - wST = 2π/NS * (rRoti + hM + hAg + hSd + 0.5*hS) # Slot pitch - - wT = wST * ratW # Tooth width - wS = wST - wT - - κS = wS/(rRoti + hAg - hRS + hSd + 0.5*hS) # Angular width of stator gaps see W. Enders thesis - wSi = κS * (rRoti + hAg + hM + hSd) # Arc length of inner part of stator - - δ = wS/wST - - # Stator back Iron height - hSBI = max(rRot/ratSplit - (rRot + hAg + hRS + hSd + hS), hS) #Added this so the min stator back iron height is hS - - # Slots/pole/phases - m = NS/(2*p*z) # Acc to Hanselman p125, this should usually be ≤ 2 - NSz = NS/z # Number of slots per phase - parte[ite_NSz] = NSz - - NSfpc = floor(m*z) - NSc = NSfpc - Nshrt # Actual number of slots per pole - cSp = NSc/NS # Short pitch correction factor - - # End-winding lengths - # Derived from (Ofori-Tenkorang1996), considering short pitch term 1/2p -> Nsct/Ns - l_ax = π * cSp * (rRoti + hAg + hM + hSd + 0.5*hS) * δ/sqrt(1-δ^2) - l_cir = 2π * cSp - - l_ewind = 2. * sqrt(l_ax^2 + (l_cir/2)^2) - # Minimum magnet skew angle to prevent Cogging Torque #[TODO] Check this - seems fishy - κMs = (NSc*2p)/ lcm(Int(NSc*2p), Int(2p)) *1/p * 180.0/π - - # Calculate correction factors - # Calculate pitch factor kp - α = π * NSc/NSfpc # Angular dispalcement of two sides of the coil - kp = sin(α/2) - # Calculation of breadth factor kb - γ = 2π * p /NS - kb = sin(m*γ/2)/(m*sin(γ/2)) - # Winding factor - kw = kp*kb - parte[ite_kw] = kw - - # Skew factor (Jagiela2013) - κMs_rad = κMs * π/180. - ks = sin(p*κMs_rad/2) / (p*κMs_rad/2) - parte[ite_ks] = ks - - # Magentic gap factor - rs = rRoti + hM + hAg - r2 = rRoti + hM - kg = ((rRoti^(p-1))/(rs^(2*p) - rRoti^(2*p))) * - ( (p/(p+1)) * (r2^(p+1) - rRoti^(p+1)) + - (p/(p-1)) * rs^(2*p) * (rRoti^(1-p) - r2^(1-p)) ) - - # Calcualte back EMF - # Accounting for slots, reluctance and leakage - Kc = 1/(1-(1/((wST/wS)*((5*hAg/wS)+1)))); - hAgEff = Kc*hAg # Effective Air-gap - Cphi = (p*κM)/180.0 # Flux concentration factor - K1 = 0.95 # Leakage factor - Kr = 1.05 # Reluctance factor - μrec = 1.05 # Recoil permeability - - PC = hM/(hAgEff * Cphi); # Permeance coefficient - BAg = ((K1 * Cphi)/(1 + (Kr * μrec/PC) ))*Br; # Flux density in the air-gap - parte[ite_BAg] = BAg - - # Calculate magnetic flux and internal voltage - κMrad = κM*(π/180); - - BC1 = (4/π)*BAg*kg*sin(p*κMrad/2); - λ = 2*rs*lRot*NSz*kw*ks*BC1/p - parte[ite_lambda] = λ #Store for off-des - - Erms = Ω*λ/sqrt(2); # RMS back voltage - - # Calculation of inductance and total reactance - # Air-gap inductance [See 6.685 notes] - LAg = (z)*(2/π)*(μ0*NSz^2*kw^2*lRot*rs)/(p^2*(hAg+hM)); - - # Slot leakage inductance - Cperm = μ0*((1/3)*(hS/wSi) + hSd/wSi); - Las = 2*p*lRot*Cperm*(4*(m-Nshrt)+2*Nshrt); - Lam = 2*p*lRot*Nshrt*Cperm; - Ls = Las - Lam; # equ. for 3 phases only - - # End-turn inductance - areaS = wS*hS; # Slot area - Le = 0.25*(μ0*wST*NSz^2)*log(wST*sqrt(π)/sqrt(2*areaS)); - #wST only true, if coils are placed in neighboring slots - - - # Total inductance and reactance per phase - Ltot = LAg + Ls + Le; - Xtot = Ω*Ltot; - - - ## ------------Calculation of machine dimensions and weights--------------- - #Armature - #Total Armature length per phase (assuming two coils per slot) - lArm=2*NSz*(lRot + l_ewind) - parte[ite_lArm] = lArm - - # Armature conductor area - areaArm = 0.5*areaS*kpf - parte[ite_areaArm] = areaArm - # Total mass of armature conductor #mass=pha*l*area*rho - mArm = z*lArm*areaArm*ρCu; - - - #Iron /Stator Core - wSd = parte[ite_wSd] - - # SBI inside radius - rSBIi = rRoti + hM + hAg + hSd + hS; - # SBI outside radius - rSBIo = rSBIi + hSBI; - # Core mass - mSBI = pi*(rSBIo^2 - rSBIi^2)*lRot*ρFe # SBI - mTeeth = (NS*wT*hS + 2*pi*(rRoti+hAg)*hSd - NS*hSd*wSd)*lRot*ρFe # Teeth - mIron = mSBI + mTeeth - - parte[ite_mSBI ] = mSBI - parte[ite_mTeeth] = mTeeth - - # Magnet mass - mM = (p*κMrad)*((rRoti+hM)^2-rRoti^2)*lRot*ρMag; - - # Shaft mass (Hollow shaft) - #tauMax=Qmax/Wt Wt=section modulus (Widerstandsmoment) - #Wt=pi/16*(da^4-di^4)/da, (Dankert,Technische Mechanik, 2018) - lShft = 1.2*lRot #A ssumption to account for bearings, endcaps, etc. - rShfto = 0.5*((16*Qmax)/(τMax*pi*(1-ratShft^4)))^(1/3) # Outer shaft radius - - mShft = (π*rShfto^2*(1-ratShft^2))*lShft*ρSteel # Mass of shaft - tShft = rShfto*(1-ratShft) # thickness of shaft, just for comparison - - - kServ = parte[ite_kServ] - # Total mass - mPMSM = kServ*(mIron + mShft + mM + mArm) - W = mPMSM * gee - parte[ite_Wpmsm] = W - - # Final machine dimensions - lPMSM = lRot + 2l_ax # Total length - dPMSM = 2*rSBIo # Total diameter without housing - - # -------------------------- - # Machine performance - # -------------------------- - ## Calculate design current and Armature resistance - #Armature RMS Current - #Irms=POutD/(z*cos(ψ)*Erms); - Irms = 1/(sqrt(2)*3)*Q/(kw*ks*NSz*BAg*(rRoti+hM)*lRot); # (Lipo2017) - - #Armature resistance per phase - θCu = parte[ite_thetaCu] - σCu = parte[ite_sigCu] - Tarm = parte[ite_Tarm] - - Rarm = lArm*(1+ θCu*(Tarm - 293.15))/(σCu*areaArm) # Pyrhönen2008 - - ##Terminal Voltage - Va = sqrt(Erms^2-((Xtot+Rarm)*Irms*cos(ψ))^2)-(Xtot+Rarm)*Irms*sin(ψ) - - - ## Loss Calculations - #Copper losses - PLCu = z*Irms^2*Rarm # Total Design Copper Losses - - #Iron losses - BSat = parte[ite_BSat] - kst = parte[ite_kst ] - pb0 = parte[ite_pb0 ] - Bb0 = parte[ite_Bb0 ] - fb0 = parte[ite_fb0 ] - ϵb = parte[ite_epsb] - ϵf = parte[ite_epsf] - - Bt = (wT + wS)/wT*BAg; #Tooth Flux Density (eq. 8.77,Lipo2017IntroToACDesign) - parte[ite_Bt] = Bt - - if Bt > BSat - println("ERROR: Bt > Saturation flux density BIronSat") - end - - Bsbi = BAg*rRoti*π/(2*p*kst*hSBI); #BackIron flux density(Hanselman eq 9.7) - parte[ite_Bsbi] = Bsbi - - if Bsbi > BSat - println("ERROR: Bsbi > Saturation flux density BIronSat") - end - - PLsbi = mSBI * pb0 * abs(Bsbi/Bb0)^ϵb * abs(f/fb0)^ϵf; # Design SBI Losses - PLt = mTeeth * pb0 * abs( Bt/Bb0)^ϵb * abs(f/fb0)^ϵf; # Design Teeth Losses - - PLiron = PLsbi + PLt # Total Design Core Losses - - #Windage loss - # [TODO] This uses air as the gas, what if we use H2 vented out of the tank... is it stupidly dangerous? Not if you have a sealed enclosure with no O₂? - Re = N*2π*rRoti*hAg*ρAir/μAir # Reynold's number in air gap - Cf = 0.0725*Re^(-0.2) # Friction coefficient - PLwind = Cf*pi*ρAir*(2π*N)^3*rRoti^4*lRot # Design Windage losses - - # Current density - JarmD = Irms/areaArm; - if JarmD > 3e7 - println("ERROR: JaD > 3e7 [A/m^2]") - end - - - ## Design Efficiency and Current Density Calculation - # Required power and efficiency - PL = PLiron + PLCu + PLwind - Preq = P + PL - η = P/Preq - parte[ite_effdes] = η - - rpm = 60N # Output rotational speed - - SP = P/mPMSM # Specific Power - parte[ite_SPdes] = SP - -return W, Preq, η, rpm, PL, PLiron/PL, PLCu/PL, PLwind/PL, SP - -end - -function PMSM(P::Float64, N::Float64, parte::Array{Float64, 1}) - - λ = parte[ite_lambda] - - kw = parte[ite_kw ] - ks = parte[ite_ks ] - NSz = parte[ite_NSz] - BAg = parte[ite_BAg] - hM = parte[ite_hM ] - lRot = parte[ite_lRot] - hAg = parte[ite_hAg ] - rRoti= parte[ite_rRoti] - - areaArm = parte[ite_areaArm] - lArm = parte[ite_lArm] - - p = parte[ite_p] - z = parte[ite_z] - - mTeeth = parte[ite_mTeeth] - mSBI = parte[ite_mSBI] - - # Calculations - ω = 2π*N - Q = P/ω - - f = p * N - Ω = 2π*f - - Erms = Ω * λ / sqrt(2) - Irms = 1/(3*sqrt(2)) * Q / - (kw * ks * NSz * BAg * (rRoti + hM)*lRot) - - #Armature resistance per phase - θCu = parte[ite_thetaCu] - σCu = parte[ite_sigCu] - Tarm = parte[ite_Tarm ] - - Rarm = lArm*(1+ θCu*(Tarm - 293.15))/(σCu*areaArm) # Pyrhönen2008 - - Jarm = Irms/areaArm - - # Losses - - Bsbi = parte[ite_Bsbi] - Bt = parte[ite_Bt ] - - pb0 = parte[ite_pb0 ] - Bb0 = parte[ite_Bb0 ] - fb0 = parte[ite_fb0 ] - ϵb = parte[ite_epsb] - ϵf = parte[ite_epsf] - - PLsbi = mSBI * pb0 * abs(Bsbi/Bb0)^ϵb * abs(f/fb0)^ϵf - PLt = mTeeth* pb0 * abs( Bt/Bb0)^ϵb * abs(f/fb0)^ϵf - PLiron = PLsbi + PLt - - Re = ω*rRoti*hAg*ρAir/μAir - Cf = 0.0725/Re^0.2 - PLwind = ρAir * π * Cf * ω^3 * rRoti^4 * lRot - - PLcu = z * Irms^2 * Rarm - - PL = PLiron + PLwind + PLcu - - Preq = P + PL - - return PL -end - -""" -inverter - -Simple inverter model that calculates the efficiency and mass of an inverter or rectifier -""" -function inverter(P::Float64, N::Float64, parte::Array{Float64, 1}) - kcf = 20. - SPinv = 19.0e3 #W/kg per GE and Uni Illinois using SiC and GaN switches respectively - https://arc.aiaa.org/doi/pdf/10.2514/6.2017-4701 - - p = parte[ite_p] - - f = p * N - fSwitch = f*kcf - η100 = -2.5e-7*fSwitch + 0.995 - η20 = η100 - 0.0017 - η10 = η100 - 0.0097 - - M = [1.0 0.1 10.0 - 1.0 0.2 5.0 - 1.0 1.0 1.0] - - k1, k2, k3 = M\[η10; η20; η100] - - # η = k1 + k2*(P/Pmax) + k3*(P/Pmax)^-1 but here P = Pmax at design - η = k1 + k2*1 + k3/1 - - parte[ite_k1] = k1 - parte[ite_k2] = k2 - parte[ite_k3] = k3 - - parte[ite_Pinvdes] = P - - Winverter = P/SPinv * gee # Weight in [N] - - return η, Winverter, SPinv - -end -""" -Off design method for inverter -""" -function inverter(P::Float64, parte::Array{Float64, 1}) - - k1 = parte[ite_k1] - k2 = parte[ite_k2] - k3 = parte[ite_k3] - - Pdes = parte[ite_Pinvdes] - - η = k1 + k2*(P/Pdes) + k3*(P/Pdes)^-1 - - return η - -end - -""" -Cable sizing - -#[TODO] See NPSS-PSL which uses real world data for typical wires - -""" -function cable(P, V, lcable, parpt) - - # Conductor: - σ = parpt[ipt_sigcon] #using σ instead of ρ (resistivity) to avoid conflict with density - α = parpt[ipt_alphacon] - Tcon = 273.15 + 50.0 - - σcon = σ/(1+α*(Tcon - 293.15)) - ρcon = parpt[ipt_rhocon] #[kg/m³] - Jmax = parpt[ipt_Jmax] - kpf = parpt[ipt_kpf] - - # Dielectric: - ρins = parpt[ipt_rhoins] - Emax = parpt[ipt_Emax] - tins = V/Emax - - - I = P/V - Acon = I/Jmax - ri = sqrt(Acon/π/kpf) # A = πr²×kpf - ro = ri + tins - Ains = π*(ro^2 - ri^2) - - Rcable = lcable/(Acon*σcon) - parpt[ipt_Rcable] = Rcable - - Lohmic = I^2 * Rcable - Pin = P + Lohmic - η = P/Pin - - parpt[ipt_rholcable] = (Acon*ρcon + Ains*ρins) - mcable = lcable*(Acon*ρcon + Ains*ρins) - Wcable = mcable*gee - - return η, Wcable -end - -function cable(Pin, parpt) - - I = Pin/parpt[ipt_Vcable] - Lohmic = I^2 * parpt[ipt_Rcable] - Pout = Pin - Lohmic - η = Pout/Pin - return η -end \ No newline at end of file diff --git a/src/propsys/propsys.jl b/src/propsys/propsys.jl index c6c8e496..9b6ff7de 100644 --- a/src/propsys/propsys.jl +++ b/src/propsys/propsys.jl @@ -10,7 +10,6 @@ export NPSS_run, startNPSS, endNPSS include("../misc/index.inc") include("../misc/constants.jl") -include("PMSM.jl") include("NPSS_functions.jl") end \ No newline at end of file diff --git a/src/sizing/wsize.jl b/src/sizing/wsize.jl index 11b971c7..29c8c002 100644 --- a/src/sizing/wsize.jl +++ b/src/sizing/wsize.jl @@ -1467,7 +1467,7 @@ function wsize(pari, parg, parm, para, pare, parg[iglnace] = lnace ipc1 = 1 - time_propsys += mission!(pari, parg, parm, para, pare, Ldebug, Nothing, Nothing, ipc1) + time_propsys += mission!(pari, parg, parm, para, pare, Ldebug) parg[igWfuel] = parm[imWfuel] # This is the design mission fuel end