Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
])
#Horizontal Tail Constraints
with SignomialsEnabled():
constraints.extend([
# Trim condition for each flight segment
TCS([cruise['x_{AC}']/aircraft.wing['mac'] <= aircraft.wing['c_{m_{w}}']/cruise['C_{L}'] + \
cruise['x_{CG}']/aircraft.wing['mac'] + aircraft.HT['V_{ht}']*(cruise['C_{L_{ht}}']/cruise['C_{L}'])]),
TCS([climb['x_{AC}']/aircraft.wing['mac'] <= aircraft.wing['c_{m_{w}}']/climb['C_{L}'] + \
climb['x_{CG}']/aircraft.wing['mac'] + aircraft.HT['V_{ht}']*(climb['C_{L_{ht}}']/climb['C_{L}'])]),
aircraft.HT['L_{ht_{max}}'] == 0.5*aircraft['\\rho_0']*aircraft['V_{ne}']**2*aircraft.HT['S_{ht}']*aircraft.HT['C_{L_{ht,max}}'],
#compute mrat, is a signomial equality
SignomialEquality(aircraft.HT['m_{ratio}']*(1+2/aircraft.wing['AR']), 1 + 2/aircraft.HT['AR_{ht}']),
#tail volume coefficient
aircraft.HT['V_{ht}'] == aircraft.HT['S_{ht}']*aircraft.HT['x_{CG_{ht}}']/(aircraft.wing['S']*aircraft.wing['mac']),
#enforce max tail location is the end of the fuselage
aircraft.HT['x_{CG_{ht}}'] <= aircraft.fuse['l_{fuse}'],
aircraft.HT['l_{ht}'] >= aircraft.HT['x_{CG_{ht}}'] - cruise['x_{CG}'],
aircraft.HT['l_{ht}'] >= aircraft.HT['x_{CG_{ht}}'] - climb['x_{CG}'],
#Stability constraint, is a signomial
TCS([aircraft['SM_{min}'] + aircraft['\\Delta x_{CG}']/aircraft.wing['mac'] + aircraft.wing['c_{m_{w}}']/aircraft.wing['C_{L_{max}}'] <= aircraft.HT['V_{ht}']*aircraft.HT['m_{ratio}'] + aircraft.HT['V_{ht}']*aircraft.HT['C_{L_{ht,max}}']/aircraft.wing['C_{L_{max}}']]),
# TCS([aircraft.wing['x_w'] >= cruise['x_{CG}'] + cruise['\\Delta x_w']]),
# TCS([aircraft.wing['x_w'] >= climb['x_{CG}'] + climb['\\Delta x_w']]),
#Mhaero >= rMh*Lhmax*(xtail-xwing), #[SP]
#Mvaero >= rMv*Lvmax*(xtail-xwing), #[SP]
hfuse == Rfuse, # may want to consider adding deltaRfuse later...
xhbend >= self.wingbox['x_{wing}'],
# Horizontal bending model
# Maximum axial stress is the sum of bending and pressurization stresses
Ihshell <= ((pi+4*thetadb)*Rfuse**2)*Rfuse*tshell + 2/3*hdb**3*tdb, # [SP]
#Ivshell <= (pi*Rfuse**2 + 8*wdb*Rfuse + (2*pi+4*thetadb)*wdb**2)*Rfuse*tshell, #[SP] #Ivshell approximation needs to be improved
sigbend == rE*sigskin,
# Horizontal bending material model
# Calculating xbend, the location where additional bending material is required
SignomialEquality(A0,A2*(xshell2-xhbend)**2 + A1*(xtail-xhbend)), #[SP] #[SPEquality]
A2 >= self.ops['N_{land}']*(Wpay+Wshell+Wwindow+Winsul+Wfloor+Wseat)/(2*lshell*hfuse*sigMh), # Landing loads constant A2
A1 >= (self.ops['N_{land}']*Wtail + rMh*self.htail['L_{h_{max}}'])/(hfuse*sigMh), # Aero loads constant A1
A0 == (Ihshell/(rE*hfuse**2)), # Shell inertia constant A0
Ahbendf >= A2*(xshell2-self.wingbox['x_f'])**2 + A1*(xtail-self.wingbox['x_f']) - A0, #[SP] # Bending area forward of wingbox
Ahbendb >= A2*(xshell2-self.wingbox['x_b'])**2 + A1*(xtail-self.wingbox['x_b']) - A0, #[SP] # Bending area behind wingbox
Vhbendf >= A2/3*((xshell2-self.wingbox['x_f'])**3 - (xshell2-xhbend)**3) \
+ A1/2*((xtail-self.wingbox['x_f'])**2 - (xtail - xhbend)**2) \
+ A0*(xhbend-self.wingbox['x_f']), #[SP]
Vhbendb >= A2/3*((xshell2-self.wingbox['x_b'])**3 - (xshell2-xhbend)**3) \
+ A1/2*((xtail-self.wingbox['x_b'])**2 - (xtail - xhbend)**2) \
+ A0*(xhbend-self.wingbox['x_b']), #[SP]
Vhbendc >= .5*(Ahbendf + Ahbendb)*self.wingbox['c_0']*self.wingbox['\\bar_w'],
Vhbend >= Vhbendc + Vhbendf + Vhbendb,
Whbend >= g*rhobend*Vhbend,
#loose constraints on speed needed to prevent N from sliding out
#to zero or infinity
N1 >= .1,
N1 <= 2,
#note residuals 2 and 3 differ from TASOPT, by replacing mhc with mlc
#in residual 4 I was able to remove the LPC/HPC mass flow equality
#in residual 6 which allows for convergence
#residual 2 HPT mass flow
TCS([mhtD == (fp1)*mhc*(Pt25/Pt41)*(Tt41/Tt25)**.5]),
#residual 3 LPT mass flow
TCS([(fp1)*mlc*(Pt18/Pt45)*(Tt45/Tt18)**.5 == mltD]),
#residual 4
SignomialEquality(u7**2 +2*Cpfanex*T7, 2*Cpfanex*Tt7),
rho7 == P7/(R*T7),
TCS([mf*(Pt2/Pref)*(Tref/Tt2)**.5 == rho7*A7*u7]),
#residual 5 core nozzle mass flow
SignomialEquality(u5**2 +2*Cptex*T5, 2*Cptex*Tt5),
rho5 == P5/(R*T5),
#compute core mass flux
mCore == rho5 * A5 * u5/(fp1),
#compute fan mas flow
mFan == rho7*A7*u7,
#residual 6 LPC/HPC mass flow constraint
TCS([mlc*(Pt18/Pref)*(Tref/Tt18)**.5 == mCore]),
#heat of combustion of jet fuel
hf = Variable('h_f', 42.8, 'MJ/kg', 'Heat of Combustion of Jet Fuel') #http://hypeRbook.com/facts/2003/EvelynGofman.shtml...prob need a better source
#engine weight
W_engine = Variable('W_{engine}', 'N', 'Weight of a Single Turbofan Engine')
pilc = Variable('\pi_{lc}', '-', 'LPC Pressure Ratio')
pihc = Variable('\pi_{hc}', '-', 'HPC Pressure Ratio')
alpha = Variable('alpha', '-', 'By Pass Ratio')
with SignomialsEnabled():
constraints = [
#making f+1 GP compatible --> needed for convergence
SignomialEquality(fp1,f+1),
#mass flow sizing
mCore == Fd/(Fsp*a0*(alphap1)), #B.194
#component area sizing
#fan area
P2 == Pt2*(hold2)**(-3.512),
T2 == Tt2 * hold2**-1,
h2 == Cp1 * T2,
rho2 == P2/(R * T2), #B.196
u2 == M2*(Cp1*R*T2/(781*units('J/kg/K')))**.5, #B.197
A2 == (alphap1)*mCore/(rho2*u2), #B.198
#HPC area
P25 == Pt25*(hold25)**(-3.824857),
T25 == Tt25 * hold25**-1,
sigth == dPover * Rfuse / tskin,
# Floor loading
lfloor >= lshell + 2 * Rfuse,
Pfloor >= Nland * (Wpay + Wseat),
Mfloor == 9. / 256. * Pfloor * wfloor,
Afloor >= 2. * Mfloor / (sigfloor * hfloor) + 1.5 * Sfloor / taufloor,
Vfloor == 2 * wfloor * Afloor,
Wfloor >= rhofloor * g * Vfloor + 2 * wfloor * lfloor * Wppfloor,
Sfloor == (5. / 16.) * Pfloor,
hfloor <= 0.1 * Rfuse,
# Tail cone sizing
taucone == sigskin,
Wcone >= rhocone * g * Vcone * (1 + fstring + fframe),
SignomialEquality(xtail, lnose + lshell + .5 * lcone), #[SP] #[SPEquality]
# BENDING MODEL
# Maximum axial stress is the sum of bending and pressurization
# stresses
Ihshell <= ((pi + 4 * thetadb) * Rfuse**2 + 8.*(1-thetadb**2/2) * (dRfuse/2.)*Rfuse + \
(2*pi + 4 * thetadb)*(dRfuse/2)**2) * Rfuse * tshell + \
2 / 3 * (hdb + dRfuse/2.)**3 * tdb, # [SP]
Ivshell <= (pi*Rfuse**2 + 8*wdb*Rfuse + (2*pi+4*thetadb)*wdb**2)*Rfuse*tshell, #[SP] #Ivshell
# approximation needs to be improved
sigbend == rE * sigskin,
# Horizontal bending material model
# Calculating xhbend, the location where additional bending
# material is required
xhbend >= xwing,
SignomialEquality(A0h, A2h * (xshell2 - xhbend) ** 2 + A1h * (xtail - xhbend)), # [SP] #[SPEquality]
Vhbendf >= A2h / 3 * ((xshell2 - xf)**3 - (xshell2 - xhbend)**3) \
+ A1h / 2 * ((xtail - xf)**2 - (xtail - xhbend)**2) \
+ A0h * (xhbend - xf), # [SP]
Vhbendb >= A2h / 3 * ((xshell2 - xb)**3 - (xshell2 - xhbend)**3) \
+ A1h / 2 * ((xtail - xb)**2 - (xtail - xhbend)**2) \
+ A0h * (xhbend - xb), # [SP]
Vhbendc >= .5 * (Ahbendf + Ahbendb) * c0 * w,
Vhbend >= Vhbendc + Vhbendf + Vhbendb,
Whbend >= g * rhobend * Vhbend,
# Wing variable substitutions
SignomialEquality(xf, xwing + .5 * c0 * \
w), # [SP] [SPEquality]
SignomialEquality(xb, xwing - .5 * c0 * \
w), # [SP] [SPEquality]
sigMh <= sigbend - rE * dPover / 2 * Rfuse / tshell,
SignomialEquality(xwing, lnose + 0.6 * lshell), # TODO remove
# Volume relations
Vcyl == Askin * lshell,
Vnose == Snose * tskin,
Vbulk == Sbulk * tskin,
Vdb == Adb * lshell,
# TODO Revert to posynomial after debugging
# [SP] #[SPEquality]
SignomialEquality(
Vcabin, Afuse * (lshell + 0.67 * lnose + 0.67 * Rfuse)),
Pcabin = Variable('P_{cabin}','Pa','Cabin Air Pressure')
W_buoy = Variable('W_{buoy}','lbf','Buoyancy Weight')
Tcabin = Variable('T_{cabin}','K','Cabin Air Temperature')
rhocabin = Variable('\\rho_{cabin}','kg/m^3','Cabin Air Density')
constraints = []
with SignomialsEnabled():
constraints.extend([
# Cabin Air properties
rhocabin == Pcabin/(state['R']*Tcabin),
Pcabin == 75000*units('Pa'),
Tcabin == 297*units('K'),
# Buoyancy weight #TODO relax the equality
SignomialEquality(W_buoy,(rhocabin - state['\\rho'])*g*aircraft['V_{cabin}']), #[SP] #[SPEquality]
# W_buoy >= (rhocabin - state['\\rho'])*g*aircraft['V_{cabin}'], # [SP]
# speed must be greater than stall speed
state['V'] >= Vstall,
# compute the drag
D >= self.wingP['D_{wing}'] + self.fuseP['D_{fuse}'] + self.VTP['D_{vt}'] + self.HTP['D_{ht}'],
LoD == W_avg/D,
# Wing looading
WLoad == .5 * self.wingP['C_{L}'] * self.aircraft['S'] * state.atm['\\rho'] * state['V']**2 / self.aircraft.wing['S'],
# Geometric average of start and end weights of flight segment
W_avg >= (W_start * W_end)**.5 + W_buoy, # Buoyancy weight included in Breguet Range
# Maximum wing loading constraint
Pcabin = Variable('P_{cabin}','Pa','Cabin Air Pressure')
W_buoy = Variable('W_{buoy}','lbf','Buoyancy Weight')
Tcabin = Variable('T_{cabin}','K','Cabin Air Temperature')
rhocabin = Variable('\\rho_{cabin}','kg/m^3','Cabin Air Density')
constraints = []
with SignomialsEnabled():
constraints.extend([
# Cabin Air properties
rhocabin == Pcabin/(state['R']*Tcabin),
Pcabin == 75000*units('Pa'),
Tcabin == 297*units('K'),
# Buoyancy weight #TODO relax the equality
SignomialEquality(W_buoy,(rhocabin - state['\\rho'])*g*aircraft['V_{cabin}']), #[SP] #[SPEquality]
# W_buoy >= (rhocabin - state['\\rho'])*g*aircraft['V_{cabin}'], # [SP]
# speed must be greater than stall speed
state['V'] >= Vstall,
# compute the drag
D >= self.wingP['D_{wing}'] + self.fuseP['D_{fuse}'] + self.VTP['D_{vt}'] + self.HTP['D_{ht}'],
LoD == W_avg/D,
# Wing looading
WLoad == .5 * self.wingP['C_{L}'] * self.aircraft['S'] * state.atm['\\rho'] * state['V']**2 / self.aircraft.wing['S'],
# Geometric average of start and end weights of flight segment
W_avg >= (W_start * W_end)**.5 + W_buoy, # Buoyancy weight included in Breguet Range
# Maximum wing loading constraint