Optmization: Need to optimize Tg in the following code w.r.t. dtc.

Ti = 251.9 + 273.15; %Fluid in degK
T0 = 386.1 + 273.15; %Fluid/gas out degK
cpl = 0.32; %Cp fluid KJ/kgdegK
LHV = 41.05; %Fuel value KJ/kg
flowL = 500; %flowrate fluid kg/s
%Air
ea = 0.15; %excess air
cpair = 0.171; %Kcal/kgdegK
Tairin = 114;
Tref = 41; %Ambient
%Fuel
THD = 129674.5; %Heat released by fuel Kw 
Ef = 0.9320; %Fuel efficiency
atfr = 1.344; %in notebook

%Dimensions
ctc = 0.3048; %ctc in m
Lh = 19.236; %Exposed tube length- horizontal m
Lv = 18.431; %Exposed tube length- vertical m
L = Lh + Lv;
d0 = 0.168; %Tube diameter m
pi = 22/7;
Tw = Tref + (Ti+T0)/L; %Tube wall temperature
Tstack = 1212; %Stack exit temperature
qfuel = THD/Ef;

mfuel = qfuel/LHV;
mair = mfuel*atfr*(1+ea);
qair = mair*cpair*(Tairin - Tref);
qloss = 0.05*qfuel;

mstack = mair + mfuel;
cpsteam = 4.1855;
cpco2 = 0.918;
cpn2 = 1.044;
xco2 = 0.054;
xsteam = 0.034; %Total should be 8.8
xn2 = 0.66;
cpgas = xsteam*cpsteam + xco2*cpco2 + xn2*cpn2;
qout = mstack*cpgas*(Tg - Tstack);

qr = qfuel + qair - (qout + qloss);
%Cold plane area
Ar = pi*L*dtc;
%absorptivity
hp = 6.626*10^-34; %m2kg/s
c = 3*10^8; %velocity of light
alpha = 4*pi*qr/hp*c;
%Exchange factor
x = pi*d0/2; %ctc
F = (sqrt(x^2 - 4) - x + 2*sinh(2/x))/2*pi;
%Heat transfer coefficient
k = 0.25*10^-3; %Thermal conductivity
rof = 1320; %density
Gm = 920880; %Dtc. Do I need to opmass flowrate
mu = 130; %Viscosity Pa-s
h0 = (do/k)*(0.023*((Gm*d0/mu*Ai)^0.8)*((cpl*mu/k))^0.3);
sig = 5.67*10^-8; %Stefan-Boltzmann constant
qr = sig*aplha*Ar*F*(Tg^4 - Tw^4) + Ar*h0*(Tg - Tw);
qfuel + qair - (mstack*cpgas*(Tg - Tstack) + qloss) = (pi*L*dtc*h0*(Tg - Tw))/(1-(sig*4*(pi^2)*dtc*L*F*((Tg^4)-(Tw^4)))/hp*c);

I need to optimize Tg w.r.t. dtc or vice versa. How do I do the same without rearranging? (Since, we cannot do that in this case)

[Matt J]

Your code doesn't show how Tg is a function of dtc.

Answer 1

[Walter Roberson]

Is

qfuel + qair - (mstack*cpgas*(Tg - Tstack) + qloss) = (pi*L*dtc*h0*(Tg - Tw))/(1-(sig*4*(pi^2)*dtc*L*F*((Tg^4)-(Tw^4)))/hp*c);

an equation to be solved for Tg ? If so, then

Tg = roots([4*F*L*Pi^2*c*cpgas*dtc*mstack*sig, -4*F*L*Pi^2*c*dtc*sig*(Tstack*cpgas*mstack+qair+qfuel-qloss), 0, 0, -(4*(F*L*Pi^2*Tw^4*c*dtc*sig+(1/4)*hp))*cpgas*mstack-L*Pi*dtc*h0*hp, (4*(F*c*Tw^3*sig*(Tstack*cpgas*mstack+qair+qfuel-qloss)*Pi+(1/4)*h0*hp))*Pi*Tw*L*dtc+hp*(Tstack*cpgas*mstack+qair+qfuel-qloss)])

It is a degree 5 polynomial, so there will be 5 outputs, and it is likely that some of them will be complex.

[Devdatt Thengdi]

Thank You Walter.

[Devdatt Thengdi]

Could you explain what you did there?

[Walter Roberson]

I did a symbolic solve

solve(qfuel+qair-mstack*cpgas*(Tg-Tstack)-qloss == pi*L*dtc*h0*(Tg-Tw)/(1-4*sig*pi^2*dtc*L*F*(Tg^4-Tw^4)*c/hp), Tg)

which returned a symbolic place holder for the roots of a quintic. I then extracted the coefficients from the place holder and put them into a vector and put roots() around it.

[Devdatt Thengdi]

No need. I got it. Now, could you please tell me how to plot this wrt to dtc(0.1:10). It says, 'Vectors must of the same length.' No I get that I am trying to plot one value of dtc against 5 values of Tg. How do I plot it for individual plots? Thanks again.

[Walter Roberson]

Unfortunately roots() does not accept an array of coefficients, so you would have to loop over all of the dtc values and roots() for each one, storing the result. Store the 5 results as 5 different columns of an array, with the array having the same number of rows as dtc has entries. Next, post-process that to drop out imaginary results,

vals(abs(imag(vals)) > 1e-10) = nan;  %arbitrarily say that small enough imaginary components are just round-off error
vals = real(vals);
plot(dtc, vals)
legend({'root 1', 'root 2', 'root 3', 'root 4', 'root 5'})

This would produce 5 different lines, some of which might turn out to have no visible content (because the root had imaginary components).

It is entirely possible that at some point you might appear to get a discontinuity where two of the lines "switch places". That is because the ordering of the roots is not well defined, other than that the first root is going to be real if there are any real roots (and even that is not promised anywhere.)

[Torsten]

For each value of dtc, you will have to figure out which of the 5 values supplied for Tg is the physical correct one.

Best wishes

Torsten.

Answer 2

[Devdatt Thengdi]

So, this is what i did:

        function [dtc,Tg] = Lobo(~)
            ctc = 0.3048; %ctc in m
            dtc= ctc:20; %min=equilateral triangle
            GAST = zeros(length(dtc),5)
            Ti = 251.9 + 273.15; %Fluid in degK
            T0 = 386.1 + 273.15; %Fluid/gas out degK
            cpl = 0.32; %Cp fluid KJ/kgdegK
            LHV = 41.05; %Fuel value KJ/kg
            flowL = 500; %flowrate fluid kg/s
            %Air
            ea = 0.15; %excess air
            cpair = 0.171; %Kcal/kgdegK
            Tairin = 114;
            Tref = 41; %Ambient
            %Fuel
            THD = 129674.5; %Heat released by fuel Kw 
            Ef = 0.9320; %Fuel efficiency
            atfr = 1.344; %in notebook

    %Dimensions
    ctc;
    disp(ctc);
    Lh = 19.236; %Exposed tube length- horizontal m
    Lv = 18.431; %Exposed tube length- vertical m
    L = Lh + Lv;
    d0 = 0.168; %Tube diameter m
    pi = 22/7;
    Tw = Tref + (Ti+T0)/L; %Tube wall temperature
    Tstack = 1212; %Stack exit temperature
    qfuel = THD/Ef;

    mfuel = qfuel/LHV;
    mair = mfuel*atfr*(1+ea);
    qair = mair*cpair*(Tairin - Tref);
    qloss = 0.05*qfuel;

    mstack = mair + mfuel;
    cpsteam = 4.1855;
    cpco2 = 0.918;
    cpn2 = 1.044;
    xco2 = 0.054;
    xsteam = 0.034; %Total should be 8.8
    xn2 = 0.66;
    cpgas = xsteam*cpsteam + xco2*cpco2 + xn2*cpn2;
    hp = 6.626*10^-34; %m2kg/s
    c = 3*10^8; %velocity of light
    %Exchange factor
    x = pi*d0/2; %ctc
    F = (sqrt(x^2 - 4) - x + 2*sinh(2/x))/2*pi;
    %Heat transfer coefficient
    k = 0.25*10^-3; %Thermal conductivity
    rof = 1320; %density
    Gm = 920880; %mass flowrate
    mu = 130; %Viscosity Pa-s
    A0 = pi*d0*d0/4;
    h0 = (d0/k)*(0.023*((Gm*d0/mu*A0)^0.8)*((cpl*mu/k))^0.3);
    sig = 5.67*10^-8; %Stefan-Boltzmann constant
    %qr = sig*alpha*Ar*F*(Tg^4 - Tw^4) + Ar*h0*(Tg - Tw);

    %Ar = pi*L*dtc; %Cold plane area
    %qout = mstack*cpgas*(Tg - Tstack); %Heat going out
    %qr = qfuel + qair - (qout + qloss);
    %alpha = 4*pi*qr/hp*c; %absorptivity

    for i = 1:length(dtc)
    Tg = roots([4.*F*L.*pi^2.*c.*cpgas.*dtc.*mstack.*sig, -4.*F.*L.*pi^2.*c.*dtc.*sig.*(Tstack.*cpgas.*mstack+qair+qfuel-qloss), 0, 0, -(4.*(F.*L.*pi^2.*Tw^4.*c.*dtc.*sig+(1/4).*hp)).*cpgas.*mstack-L.*pi.*dtc.*h0.*hp, (4.*(F.*c.*Tw^3*sig.*(Tstack.*cpgas.*mstack+qair+qfuel-qloss).*pi+(1/4)*h0*hp)).*pi.*Tw.*L*dtc+hp.*(Tstack.*cpgas.*mstack+qair+qfuel-qloss)]);
    GAST(i) = Tg;
    end

    vals = Tg;
    vals(abs(imag(vals)) > 1e-10) = nan; %excluding imaginary valuses as round-off errors 
    vals = real(vals);

    Tg;
    disp(Tg);

    plot (dtc,Tg)
    legend ({'root1' 'root2' 'root3' 'root4' 'root5'})
    xlabel('Tube circle diameter (d0)');
    ylabel('Gas temperarture (Tg)');
    end

This is the result:

>> Lobo

GAST =

       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0
       0     0     0     0     0

      0.3048

In an assignment  A(I) = B, the number of elements in B and I must be the
same.

Error in Lobo (line 68)
GAST(i) = Tg; 
Also, I am getting a blank plot.