diff --git a/src/material_data/MaterialProperties.toml b/src/material_data/MaterialProperties.toml index 026f6430..57fd9891 100644 --- a/src/material_data/MaterialProperties.toml +++ b/src/material_data/MaterialProperties.toml @@ -233,6 +233,26 @@ density and 12% higher modulus, as well as other improved properties """ density = 1700.0 dielectric_strength = 10e6 +# ----------------------------------------------------- +# Electric Steels +# ----------------------------------------------------- +[M19] + description = """29-Gauge M-19 nonoriented steel; + Data from Ofori-Tenkorang, J., + “Permanent-Magnet Synchronous Motors and Associated Power Electrics + for Direct-Drive Vehicle Propulsion,” Ph.D. Dissertation, + Massachusetts Institute of Technology, Cambridge, MA, 1996.""" + density = 7750.0 #kg/m3 + ke = 7.0963515e-5 # W/kg/Hz²/T² = 32.183e-6 W/lb/Hz²/T² + kh = 2.351412e-2 # W/kg/Hz = 10.664e-3 W/lb/Hz + alpha = 1.793 +[Hiperco50] + description = """Hiperco® 50 alloy is an iron-cobalt-vanadium soft magnetic + alloy which exhibits high magnetic saturation (24 kilogauss), + high D.C. maximum permeability, low D.C. coercive force, and low A.C. core loss.""" + density = 8110.0 + + # ----------------------------------------------------- # Thermal insulators # ----------------------------------------------------- diff --git a/src/misc/materials.jl b/src/misc/materials.jl index d755e6b5..b7cac380 100644 --- a/src/misc/materials.jl +++ b/src/misc/materials.jl @@ -6,7 +6,9 @@ module materials using TOML, DocStringExtensions -export StructuralAlloy, Conductor, Insulator, ThermalInsulator, thermal_conductivity +export StructuralAlloy, Conductor, Insulator, ElectricSteel, +ThermalInsulator, thermal_conductivity +export resistivity, resxden __abs_path_prefix__ = dirname(@__DIR__) MaterialProperties = TOML.parsefile(joinpath(__abs_path_prefix__,"material_data/MaterialProperties.toml")) @@ -176,6 +178,56 @@ function Insulator(material::String) end end +""" +$TYPEDEF + +ElectricSteel. + +$TYPEDFIELDS +""" +@kwdef struct ElectricSteel + """Name""" + name::String = "" + """Density [kg/m³]""" + ρ::Float64 + """Eddy current loss coefficient [W/lbm/Hz²/T²]""" + kₑ::Float64 + """Hysteresis loss coefficient [W/lbm/Hz]""" + kₕ::Float64 + """Exponential fit coefficient for hysteresis loss""" + α::Float64 +end +""" + ElectricSteel(material::String) + +Outer constructor for `ElectricSteel` types. +Material specified needs to have the following data in the database: +- ρ (density): Density [kg/m³] +- ke +- kh +- α +""" +function ElectricSteel(material::String) + local MatProp, ρ, ke, kh, α + try + MatProp = MaterialProperties[material] + catch + error("Cannot find $material in Material Properties database") + else + try + ρ = MatProp["density"] + ke = MatProp["ke"] + kh = MatProp["kh"] + α = MatProp["alpha"] + catch + error("Insufficient data in database for $material to build a Conductor") + else + Conductor(material, ρ, ke, kh, α) + end + end + +end + """ resxden(cond::conductor) diff --git a/src/propsys/PMSM.jl b/src/propsys/PMSM.jl new file mode 100644 index 00000000..f2733834 --- /dev/null +++ b/src/propsys/PMSM.jl @@ -0,0 +1,409 @@ +module ElectricMachine + +abstract type AbstractElectricMachine end +abstract type AbstractMagnets end + +const μ₀ = 1.25663706127e-6 #N⋅A⁻² https://physics.nist.gov/cgi-bin/cuu/Value?mu0 + +using ..materials +using DocStringExtensions + +@kwdef struct PermanentMagnet <: AbstractMagnets + """Magnet thickness [m]""" + thickness::Float64 = 20e-3 + """Magnet Density [kg/m⁻³]""" + density::Float64 = 7501.0 #Neodymium magnets + """Magnetization constant [A/m]""" + M::Float64 = 8.604e5 + """Mass of magnets [kg]""" + mass::Float64 = 0.0 +end + +@kwdef struct Windings + """Winding conductor material""" + conductor::Conductor + """Winding insulator material""" + insulator::Insulator + """Number of turns [-]""" + N_turns::Int = 1 + """Wire area [m²]""" + A_wire::Float64 = 0.518e-6 #20AWG + """Winding temperature [K]""" + T::Float64 = 273.15 + 90 + """Packing factor [-]""" + kpf::Float64 = 0.35 + """Mass of windings [kg]""" + mass::Float64 = 0.0 +end + +@kwdef mutable struct Teeth + """Tooth width [m]""" + width::Float64 = 0.005 + """Tooth thickness [m]""" + thickness::Float64 = 0.02 + """Inner radius [m]""" + Ri::Float64 = 0.01 + """Flux density [Wb/m² = T]""" + B::Float64 = 1.0 + """Material""" + material::ElectricSteel = ElectricSteel("M19") + """Mass [kg]""" + mass::Float64 = 0.0 +end + +@kwdef mutable struct RotorGeometry + """Rotor thickness [m]""" + thickness::Float64 = 0.01 + """Inner radius [m]""" + Ri::Float64 = 0.01 + """Outer radius [m]""" + Ro::Float64 = 0.02 + """Flux density [Wb/m² = T]""" + B::Float64 = 1.0 + """Material""" + material::ElectricSteel = ElectricSteel("M19") + """Mass [kg]""" + mass::Float64 = 0.0 +end + +@kwdef mutable struct StatorGeometry + """Stator thickness [m]""" + thickness::Float64 = 0.01 + """Inner radius [m]""" + Ri::Float64 = 0.01 + """Outer radius [m]""" + Ro::Float64 = 0.02 + """Flux density [Wb/m² = T]""" + B::Float64 = 1.0 + """Material""" + material::ElectricSteel = ElectricSteel("M19") + """Mass [kg]""" + mass::Float64 = 0.0 +end + +@kwdef mutable struct ShaftGeometry + """Inner Radius [m]""" + Ri::Float64 = 0.01 + """Outer Radius [m]""" + Ro::Float64 = 0.01 + """Overhang fraction of shaft [-]""" + l_extra::Float64 = 1.2 + """Material""" + material::StructuralAlloy + """Mass [kg]""" + mass::Float64 = 0.0 +end + +""" +$TYPEDEF + +$TYPEDFIELDS + +Structure that defines a permanent magnet synchronous motor (PMSM). +""" +@kwdef mutable struct Motor + rotor::RotorGeometry = RotorGeometry() + stator::StatorGeometry = StatorGeometry() + teeth::Teeth = Teeth() + magnet::PermanentMagnet = PermanentMagnet() + windings::Windings = Windings() + + """Design shaft power [W]""" + P_design::Float64 = 1e6 + """Angular velocity [rad/s]""" + Ω::Float64 = 10e3 + """Frequency [Hz]""" + f::Float64 = 12e3 + """Slots per pole [-]""" + Nsp::Int = 3 + """Phases [-]""" + phases::Int = 3 + """Phase resistance [Ω]""" + phase_resistance::Float64 = 0.0 + """Design current [A]""" + Id::Float64 = 0.0 + """Design Voltage [V]""" + Vd::Float64 = 0.0 + """Current [A]""" + I::Float64 = 0.0 + """Voltage [V]""" + V::Float64 = 0.0 + + """Pole pairs [-]""" + N_pole_pairs::Int = 8 + """Max rotor tip speed [m/s]""" + U_max::Float64 = 200.0 + """Saturation Flux [Wb/m²]""" + B_sat::Float64 = 1.8 + """Maximum Current Density [A/m²]""" + J_max::Float64 = 1e6 + """Max ratio to saturation""" + rB_sat::Float64 = 0.98 #Default assumption is that the metal is almost saturated at design + """Stack length [m]""" + l::Float64 = 1.0 + """Mean internal temperature [K], + used to calculate air density and viscosity, and winding resistance""" + T::Float64 = 273.15 + 90 + """Thickness of air gap [m]""" + airgap_thickness::Float64 = 2e-3 + """Pole pairs [-]""" + p::Int = 8 + """Mass of machine [kg]""" + mass::Float64 = 0.0 +end + +# Layout: +# +# +-------------------------------+ <----------------+ +# | Stator back Iron | }ts | +# | +--+ +--+ +--+ +--+ | <------------+ | +# Teeth | | | | | | | | | wt| }tt | | +# *---* *---* *---* *---* *---* <--------+ | | +# air gap { | | | +# (g) +---------------+---------------+ ^ <----+ | | | +# | Magnets | | tm | | | | +# | | | v | | | | +# +---------------+---------------+ <-+ | | | | +# | | ^ | | | | | +# | Rotor back Iron | tr | | | | | +# | | v | | | | | +# +-------------------------------+ | | | | | +# ^ | | | | | +# | | | | | | +# + + + + + + +# Rri Rro Rg Rt Rsi Rso + + + + +""" +$TYPEDEF + +Simplified permanent magnet synchronous machine (PMSM) sizing. +""" +function size_PMSM!(motor::Motor, shaft_speed::AbstractFloat, design_power::AbstractFloat) + rotor = motor.rotor + stator = motor.stator + magnet = motor.magnet + teeth = motor.teeth + windings = motor.windings + shaft = motor.shaft + + motor.Ω = shaft_speed * (2 * π / 60) + motor.f = motor.N_pole_pairs * shaft_speed / 60 + + #Outer radius of the magnet - subject to the highest rotational speeds + radius_gap = motor.U_max / motor.Ω + + # Calculate maximum flux through the stator and rotor back iron + # Really, this needs to be a constraint Bsbi ≤ Bsat (at the top level) rather + # than a prescribed setting. Or rB_sat can be a global optimization variable. + # The higher Bsbi, the thinner the sbi, but higher Bsbi ⟹ higher core losses + stator.B = motor.rB_sat * motor.B_sat + rotor.B = motor.rB_sat * motor.B_sat + + B_gap = airgap_flux(magnet.M, magnet.thickness, motor.airgap_thickness) + + stator.thickness = B_gap / stator.B * π * radius_gap / (2 * motor.N_pole_pairs) + rotor.thickness = B_gap / rotor.B * π * radius_gap / (2 * motor.N_pole_pairs) + + rotor.Ro = radius_gap - magnet.thickness + rotor.Ri = rotor.Ro - rotor.thickness + + teeth.Ri = radius_gap + motor.airgap_thickness + + stator.Ri = teeth.Ri + teeth.thickness + stator.Ro = stator.Ri + stator.thickness + + N_teeth = motor.Nsp * 2 * motor.N_pole_pairs + + A_teeth = N_teeth * teeth.thickness * teeth.width + A_slots = π * (teeth.Ri + stator.Ri) * teeth.thickness - A_teeth + A_slot = A_slots / N_teeth + λ = A_slots / (A_slots + A_teeth) + l_end_turns = π / (2 * motor.N_pole_pairs) * teeth.Ri * (λ / sqrt(1 - λ^2)) + + #Armature reaction + B_windings = μ₀ * motor.J_max * teeth.thickness + + teeth.B = sqrt((B_gap / λ)^2 + B_windings^2) + if teeth.B > motor.B_sat + @warn "Flux density in teeth exceeds saturation flux density" + end + + torque = design_power / motor.Ω + slot_current = J_max * windings.kpf * A_slot #Peak slot current + windings.N_turns = floor(windings.kpf * A_slot / windings.A_wire) + + #Assume 2 phases excited at any given time + motor.I = + motor.Id = (2 / motor.phases * motor.Nsp * 2 * motor.N_pole_pairs) * slot_current + force_length = motor.Id * B_gap + motor.l = torque / (force_length * radius_gap) + + #Size shaft + shaft.Ro = rotor.Ri + if shaft.Ro^4 > torque / shaft.material.τmax * 2 * shaft.Ro / π + shaft.Ri = (shaft.Ro^4 - torque / shaft.material.τmax * 2 * shaft.Ro / π)^0.25 + else + # shaft can be solid as well + shaft.Ri = 0.0 + end + + #Masses + rotor.mass = cross_sectional_area(rotor) * motor.l * rotor.material.density + stator.mass = cross_sectional_area(stator) * motor.l * stator.material.density + magnet.mass = + cross_sectional_area(radius_gap, rotor.Ro) * motor.l * magnet.material.density + teeth.mass = A_teeth * motor.l * teeth.material.density + + total_winding_volume = A_slots * (motor.l + 2 * l_end_turns) + windings.mass = + windings.kpf * total_winding_volume * windings.conductor.material.density + + (1 - windings.kpf) * total_winding_volume * windings.insulator.material.density + + shaft.mass = + cross_sectional_area(shaft) * (motor.l * shaft.l_extra) * shaft.material.density + + motor.mass = + rotor.mass + stator.mass + magnet.mass + teeth.mass + windings.mass + shaft.mass + #size_rotor!(rotor, R_gap) + #size_stator!() + #size_windings!() + + ## Calculate losses + windings.T = motor.T # Set temperature + motor.phase_resistance = + motor.Nsp / motor.phases * + 2 * + motor.N_pole_pairs * + slot_resistance(motor.windings, windings.kpf * A_slot, motor.l + 2 * l_end_turns) + + +end # function size_PMSM! + +""" +""" +function calculate_losses(EM::Motor) + Q_ohmic = ohmic_loss(EM.I, EM.phase_resistance) + Q_core = core_loss(EM) + Q_wind = windage_loss(EM) + return Q_ohmic + Q_core + Q_wind +end # function calculate_losses + +""" + +Returns the power dissipated due to resistive heating assuming a certain +number of phases are energized at any given time. By default assumes that 2 +phases are energized at any given time. +""" +function ohmic_loss( + I::AbstractFloat, + phase_resistance::AbstractFloat; + energized_phases::Int = 2, +) + return I^2 * (energized_phases * phase_resistance) +end # function ohmic_loss + +core_loss(motor) = hysteresis_loss(motor) + eddy_loss(motor) + +""" +$TYPEDSIGNATURES + +Calculates the hysteresis loss in the motor (i.e., in the rotor + the stator) +""" +function hysteresis_loss(motor::Motor) + return hysteresis_loss(motor.stator, motor.f) + hysteresis_loss(motor.teeth, motor.f) +end # function hysteresis_loss + +""" +Calculates the hysteresis losses in steel. The steel needs to be of type +[`ElectricSteel`](@ref). +""" +function hysteresis_loss(steel, f) + return steel.mass * steel.material.kₕ * f * steel.B^steel.material.α +end + +""" +$TYPEDSIGNATURES + +Calculates the eddy losses in the motor (i.e., due to the rotor + the stator) +""" +function eddy_loss(motor::Motor) + return eddy_loss(motor.stator, motor.f) + eddy_loss(motor.teeth, motor.f) +end #function eddy_loss + +""" +Calculates the eddy losses in an electrical steel. The steel needs to be of type +[`ElectricSteel`](@ref). +""" +function eddy_loss(steel, f) + return steel.mass * steel.material.kₑ * f^2 * steel.B^2 +end # function eddy_loss + +""" +Calculates the windage (i.e., air friction losses) for the rotor. +""" +function windage_loss(motor::Motor, air) + #air from IdealGases? + Re = motor.Ω * motor.radius_gap * motor.airgap_thickness * air.ρ / air.μ + Cf = 0.0725 * Re^-0.2 #avg. skin friction approx. from fully turbulent flat plate + return Cf * π * air.ρ * motor.Ω^3 * motor.radius_gap^4 * motor.l +end # function windage_loss + +""" +Returns the slot resistance given a winding temperature and material. +The winding material needs to be of type [`Conductor`](@ref). +""" +function slot_resistance(wind::Windings, A, l) + resist = resistivity(wind.conductor, wind.T) + return resist * l / A +end # function slot_resistance + +""" +$TYPEDSIGNATURES + +Returns the cross sectional area of an annulus with inner radius Rᵢ and outer radius Rₒ. +""" +function cross_sectional_area(Ro, Ri) + return π * (Ro^2 - Ri^2) +end # function cross_section + +# Might not need to define each of this and just let anything go +# cross_sectional_area(X) = cross_sectional_area(X.Ro, X.Ri) +# if it doesn't have the right fields it'll throw an error. +cross_sectional_area(R::RotorGeometry) = cross_sectional_area(R.Ro, R.Ri) +cross_sectional_area(R::StatorGeometry) = cross_sectional_area(R.Ro, R.Ri) +cross_sectional_area(R::ShaftGeometry) = cross_sectional_area(R.Ro, R.Ri) + + +# """ +# """ +# function size_rotor!(rotor::RotorGeometry, R_gap::AbstractFloat) +# Ro = R_gap - magnet.thickness +# rotor.thickness = (B_gap / rotor.B) * π * R_gap / (2p) +# rotor.Ri = Ro + +# end # function size_rotor! + +""" + airgap_flux(M, thickness, airgap) + +Given the magnetization constant (M), thickness of the magnet, and the +air gap thickness, returns the flux density in the air gap. +Simplified expression from a magnetic circuit analysis. Assumes no MMF drop +across the steel. +""" +function airgap_flux(M, thickness, airgap) + B_gap = μ₀ * M * thickness / (thickness + airgap) +end # function airgap_flux + +""" +Linear model for magnet remanent flux as a function of temperature. +This will be important if we model the TMS system along with the PMSM. +""" +function remanent_flux(magnet::AbstractMagnets, T) + magnet.remanent_flux * (1 - magnet.α * (T - (273.15 + magnet.Tbase)) / 100.0) +end # function remanent_flux + +end diff --git a/src/propsys/cable.jl b/src/propsys/cable.jl new file mode 100644 index 00000000..12e78f60 --- /dev/null +++ b/src/propsys/cable.jl @@ -0,0 +1,82 @@ +""" +$TYPEDEF + +Simplified model of a stranded (i.e., Litz) wire with an outer insulation layer. +Default values are for copper based on https://pubs.aip.org/aip/jpr/article/8/4/1147/242248/Electrical-resistivity-of-copper-gold-palladium + +# Fields +$TYPEDFIELDS + +""" +@kwdef mutable struct cable + """Conductor""" + cond::Conductor = Conductor("Cu") + """Insulator""" + ins::Insulator = Insulator("polyamide") + """Temperature of conductors [K]""" + Tcon::Float64 = 273.1 + 50.0 #defaults to 50ᵒC + """Max current density [A/m²]""" + Jmax::Float64 = 5e6 # 5 A/mm² = 5e6 A/m² + """Packing factor [-]""" + kpf::Float64 = 0.91 + """Length [m]""" + l::Float64 = NaN + """Mass [kg]""" + mass::Float64 = NaN + """Weight [N]""" + W::Float64 = NaN + """Resistance [Ω]""" + R::Float64 = NaN + """Voltage [V]""" + V::Float64 = NaN +end + +""" + (c::cable)(P::Float64, V::Float64, lcable::Float64) + +Sizes the cable given a power, voltage, and length of the cable and returns +the *efficiency function* for the sized cable. +""" +function (c::cable)(P::Float64, V::Float64, lcable::Float64) + cond = c.cond + ins = c.ins + c.V = V + c.l = lcable + # Calculate current + I = P/V + + # Get resistivity at given temperature + ρcon = resistivity(cond, c.Tcon) + + # Geometry + tins = V/ins.Emax + Acon = I/c.Jmax + ri = sqrt(Acon/(π*c.kpf)) # A = πr²×kpf + ro = ri + tins + Ains = π*(ro^2 - ri^2) + c.mass = lcable*(Acon*cond.ρ + Ains*ins.ρ) + c.W = c.mass*gee + + # Resistance + c.R = ρcon * lcable/Acon + + # Loss function + """ + Returns the efficiency of the cable for the given power input assuming constant voltage. + """ + function f(Pin::Float64) + + Iin = Pin/c.V + + if Pin>P + @warn """Power ($(Pin/1e3) kW) exceeds design power for cable ($(P/1e3) kW). + Current density = $(Iin/Acon) > Jmax ($(c.Jmax)) A/m²""" + end + + Lohmic = Iin^2 * c.R # = I²R loss + Pout = Pin - Lohmic + η = Pout/Pin + return η + end # function + return f +end \ No newline at end of file diff --git a/src/propsys/inverter.jl b/src/propsys/inverter.jl new file mode 100644 index 00000000..91490653 --- /dev/null +++ b/src/propsys/inverter.jl @@ -0,0 +1,72 @@ +""" +$TYPEDEF + +Simplified model of an inverter. + +$TYPEDFIELDS +""" +@kwdef mutable struct inverter + """Design power [W]""" + P::AbstractFloat = 1e3 + """Number of pole-pairs""" + N_poles::Integer = 2 + """Design Rotational speed [RPM]""" + N::AbstractFloat + +end + +""" + inverter(P::Float64, N::Float64, parte::Array{Float64, 1}) + + +Simple inverter model that calculates the efficiency and mass of an inverter or rectifier. +Default specific power is 19 kW/kg per [GE and Uni Illinois using SiC and GaN switches respectively.](https://arc.aiaa.org/doi/pdf/10.2514/6.2017-4701) +""" +function inverter(P::Float64, N::Float64, parte::Array{Float64, 1}) + kcf = 20. + SPinv = 19.0e3 #W/kg + + 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 \ No newline at end of file diff --git a/src/propsys/propsys.jl b/src/propsys/propsys.jl index 657af2ab..95c0255b 100644 --- a/src/propsys/propsys.jl +++ b/src/propsys/propsys.jl @@ -1,16 +1,21 @@ """ Module containing all propulsion system related code. -Will cover alternate engine models - NPSS vs Drela's orig. model vs pyCycle +This includes alternate propulsion systems such as turbo-electric architectures. +Will eventually cover alternate engine models - NPSS vs Drela's orig. model vs pyCycle, +replacing the turbofan code such as `tfoper`, `tfsize`etc. """ module propsys export NPSS_run, startNPSS, endNPSS import ..TASOPT: __TASOPTindices__, __TASOPTroot__ +using ..materials +using DocStringExtensions include(__TASOPTindices__) include(joinpath(__TASOPTroot__,"misc/constants.jl")) - +include("cable.jl") +include("PMSM.jl") include("NPSS_functions.jl") end \ No newline at end of file