clc;
clear;
close all;
tic
ML = 1 ;
hmax = 20*10^3;
amax = 10;
SM = 2;
rhos = 2700 ;
rhop = 1772 ;
sigmas = 60*10^6 ;
N = 3 ;
g = 9.81 ;
theta = 0;
ct = cosd(theta);
T = 298 ;
Pa = 101.325*10^3 ;
gamma = 1.4 ;
R_gc = 287 ;
a_ss = sqrt(gamma*R_gc*T) ;
Rmax = 1+amax ;
Weq = sqrt((hmax*g)./((log(Rmax)/2).* ...
(log(Rmax)-2)+((Rmax-1)./Rmax))) ;
tb = (Rmax-1)*Weq/(g*Rmax) ;
Meq = Weq/a_ss ;
P0_Pa = (1+Meq^2*((gamma-1)/2))^(gamma/(gamma-1)) ;
P0 = P0_Pa*Pa ;
D = 0.0922;
delta = D*P0/(2*sigmas) ;
lambdamax = 0;
while true
syms Li X_cp X_cg
t=D*P0/(2*sigmas);
Mn = delta*rhos*pi*D*(D+(sqrt(D^2+(D/2)^2))) ;
Vf = 0.5*D*D*delta ;
Mf = 3*(rhos*Vf) ;
Mfb = pi*D*rhos*D*delta ;
Mcyl = (2*pi*(D/2)*(Li+D)+2*pi*(D^2/4))*rhos*delta ;
Mfin = (3/2)*D*rhos*delta ;
Mcone = pi*(D/2)*((D/2)+sqrt(D^2 + (D^2/4)))*rhos*delta ;
Ms = Mcyl+Mfin+Mcone ;
R = Rmax ;
Mp = ((R - 1)*(Ms + ML)) ;
Lp = Mp/(pi*D^2*rhop/4) ;
Ckn = (4*N*(4/3))/(1+sqrt(6)) ;
Xcp = (1.33*D + Ckn*(D+Li+(D/3)))/(2+Ckn) ;
Xcgnum = Mcone*(2/3)*D + Mcyl*((Li+D)/2)+Mp*(2*D+Li-Lp+(Lp/2))...
+ Mfin*(D+Li+(1/3)*D);
Xcgden = Mcone+Mcyl+Mfin+Mp ;
Xcg = Xcgnum/Xcgden ;
eqn1 = Xcp-Xcg-D*SM == 0 ;
L1 =vpasolve(eqn1,Li);
L=L1(real(L1)>0);
Mn = delta*rhos*pi*D*(D+(sqrt(D^2+(D/2)^2))) ;
Vf = 0.5*D*D*delta ;
Mf = 3*(rhos*Vf) ;
Mfb = pi*D*rhos*D*delta ;
Mcyl = (2*pi*(D/2)*(L+D)+2*pi*(D^2/4))*rhos*delta ;
Mfin = (3/2)*D*rhos*delta ;
Mcone = pi*(D/2)*((D/2)+sqrt(D^2 + (D^2/4)))*rhos*delta ;
Ms = Mcyl+Mfin+Mcone ;
R = Rmax ;
Mp = vpa(((R - 1)*(Ms + ML))) ;
Lp = Mp/(pi*D^2*rhop/4) ;
Ckn = (4*N*(4/3))/(1+sqrt(6)) ;
Xcpi = vpasolve((1.3333*D + (Ckn*(D+L+D./2)))/(2+Ckn) - X_cp==0);
Xcp = Xcpi(real(Xcpi)>0) ;
Xcgnum = Mcone*(2/3)*D + Mcyl*((L+D)/2)+Mp*(2*D+L-Lp+(Lp/2))...
+ Mfin*(D+L+(1/3)*D);
Xcgden = Mcone+Mcyl+Mfin+Mp ;
Xcgi = vpasolve((Xcgnum/Xcgden)-X_cg) ;
Xcg = Xcgi(real(Xcgi)>0) ;
Lpf = L+D ;
lambda = ML/(Ms+Mp) ;
M0 = ML+Ms+Mp ;
if (L>D && Xcp>Xcg && Lp<Lpf)
if (lambda >= lambdamax)
Lnew = L ;
Lpnew = Lp ;
Dnew = D ;
tnew = t ;
Xcpnew = Xcp ;
Xcgnew = Xcg ;
Msnew = Ms ;
M0new = M0 ;
Mpnew = Mp ;
lambdamax_new = lambdamax ;
lambdamax = lambda ;
D = D+0.001 ;
else
Lp = Lpnew ;
t = tnew ;
Xcp = Xcpnew ;
Xcg = Xcgnew ;
L = Lnew ;
D = Dnew ;
lambdamax = lambdamax_new ;
Ms = Msnew ;
M0 = M0new ;
Mp = Mpnew ;
break ;
end
else
Lp = Lpnew ;
t = tnew ;
Xcp = Xcpnew ;
Xcg = Xcgnew ;
L = Lnew ;
D = Dnew ;
lambdamax = lambdamax_new ;
Ms = Msnew ;
M0 = M0new ;
Mp = Mpnew ;
break ;
end
end
Impluse = Weq/g;
m_dot = Mp/tb;
Thrust = m_dot*Impluse*g
pressure_R=P0/Pa
tD_r= t/D
t
D
L
Lp
ld_r =L/D
Xcg
Xcp
Mp
Ms
M0
lambdamax
Eplison = lambda_m*M_s/(M_l)
toc