MATLAB: For loop within fzero in 1D heat transfer problem

1d heat transferfinite elementfor loopfzeroheat transferimplicit function

Hi, I am trying to plot the temperature along a 1D length scale using finite element method on nodes, i.e. the temperature of node i-1 is the boundary condition for node i-1. fzero worked to give the first value before I introduced the for loop. The error 'Index exceeds matrix dimensions' appeared when I started using the for loop. Is there anyway to fix this?

%constants: 
h=6.47;   
Lc=0.0168; %distance between nodes  
L=1.274; %total length  
Ta=20+273;   
Tb=38+273;   
As=pi*0.0418^2;  
dt=10;  
Q=-10;  
X=0:L/Lc;  
x=length(X);  
T=zeros(x);  
T(1)=39+273.15; %initial condition  
for i=2:1:x  
fcn=@(T) 0.39*As*(T(i)-T(i-1))/Lc+2*pi*0.0418*Lc*(h*(T(i)-Ta)+0.99*5.67e-8*(T(i).^4-Ta^4))+As*Lc*1069*3650*3.28e-4*(Tb-T(i))+398*As*Lc+398*(2.^(0.1*(T(i)-303.35))-1)*As*Lc-3057*As*Lc*1052*(T(i)-T(i-1))/dt+Q; 
T=fzero(fcn,T(i-1));    <--Error: Index exceeds matrix dimensions  
end

Thank you very much!

Best Answer

%constants: 
h=6.47;   
Lc=0.0168; %distance between nodes  
L=1.274; %total length  
Ta=20+273;   
Tb=38+273;   
As=pi*0.0418^2;  
dt=10;  
Q=-10;  
X=0:L/Lc;  
x=length(X);  
T=zeros(1,x);  
T(1)=39+273.15; %initial condition  
for i=2:1:x  
fcn=@(t) 0.39*As*(t-T(i-1))/Lc+2*pi*0.0418*Lc*(h*(t-Ta)+0.99*5.67e-8*(t.^4-Ta^4))+As*Lc*1069*3650*3.28e-4*(Tb-t)+398*As*Lc+398*(2.^(0.1*(t-303.35))-1)*As*Lc-3057*As*Lc*1052*(t-T(i-1))/dt+Q; 
T(i)=fzero(fcn,T(i-1));  
end

Related Solutions

MATLAB: It gives Invalid use of operator error

Try this:


clc; clear all ;
close all;
clear;
clc;
umatrix = zeros(200, 1);
Cmatrix = zeros(200, 1);
j = 1;
% (1) Get t and di
Pp = 15.5; % Primary pressure
do = 0.01748; % Outer tube diameter
S = 160.65; % Allowable stress
y = 0.4; % Coeffcient given by ASME boiler code
t = Pp * do/(2 * S + 2 * y * Pp); % Wall thickness
di = do - 2 * t; % Inner tube diameter
% (2) With guessed up, get the total number of tube Nsg
mp = (1.38 * 10^4) / 3; % Mass flow rate into one SG
rhop = 706.917; % Primary density, T = 309.055 C and P = 155 bar
for up = 0.1:0.1:20 % To find the optimal velocity
    Nsg = 1.2 * mp / (rhop * up * di^2 * pi / 4); % Number of SG tubes
    % (3) Get Tpout, Rpmf, q''(z=0), q''(z=Lsg)
    
    % To calculate Rpmf
    % Calcuation Rp
    kp = 0.54471; % Thermal conductivity of primary system
    Ap = pi * di^2 / 4; % Flow area of one tube of primary system
    mup = 84.9498 * 10^(-6); % Viscosity of primary system
    cp = 5.71027 * 10^3; % Heat capacity of primary system
    
    hp = (kp / di) * 0.023 * ((mp * di) / (Ap * Nsg * mup))^(0.8) * (cp * mup / kp)^(0.33); % HTC of primary system
    Rp = (do / di) * (1 / hp); % Primary heat transfer resistance
    
    % Calculation Rw
    kw = 12.552; % Wall thermal conductivity of Inconel 690
    Rw = (do * log(do / di)) / (2 * kw); % Tube wall resistance
    
    % Calculation Rf
    Rf = 1.057 * 10^(-5); % It is given

    
    % Calculation Rs
    hf = 1249.38 * 1000; % Enthalpy of fluid
    hg = 2777.06 * 1000; % Enthalpy of gas
    hfg = hg - hf; % Unit : J/kg
    Prf = 0.855144 % Prandtl number of fluid
    cf = 5.33358 * 1000; % Heat capacity of fluid, unit : J/kg
    Cws = 0.006; % It is given
    muf = 92.5798 * 10^(-6); % Dynamic liquid viscosity
    sigma = 0.0184373; % Surface tension
    g = 9.80665; % Gravitational acceleration
    rhof = 746.02; % Density of liquid
    rhog = 34.4984; % Density of gas
    kf = 0.577425; % Thermal conductivity of the liquid
    
    F = (hfg / cf) * Cws * (1 / (muf * hfg) * ((sigma / (g * (rhof - rhog))^(1/2)))^(1/3)) * ((cf * muf / kf)^(1.7)); % Rohsenow correlation
    Tp1 = 326.61; % Tin
    Tp2 = 291.50; % Tout
    Ts = 282.379; % Saturation temperature at P = 6.65MPa
    Rc = Rp + Rw + Rf;
    
    % q''(z=0), q''(z=Lsg)
    syms Q1 Q2
    Q1 = solve((Tp1 - Ts) == Rc * Q1 + F * Q1^(1/3), Q1); % Q1 = q''(z=0)
    Q2 = solve((Tp2 - Ts) == Rc * Q2 + F * Q2^(1/3), Q2); % Q2 = q''(z=Lsg)
    % (4) Get Asg and Lsg
    
    Asg = mp * cp * (Rc * log(Q1 / Q2) + (1/2) * F * (Q2^(-2/3) - Q1^(-2/3))); % SG wall area
    Lsg = (Asg / (Nsg * do * pi)); % Length of a tube
    % (5) Get pumping power
    epsilon = 1.01 * 10^(-5); % SG tube surface roughness
    K = 3; % Total form loss coefficient
    mp = Ap * Nsg * rhop * up; % mp = rhoAg
    f = 0.0055 * (1 + ((2 * 10^4 * epsilon / di) + (10^6 * Ap * mup / ((mp / Nsg) * di)))^(0.333)); % Rough condition
    Ploss = (((K + ((f * Lsg) / di)) * up^2 * rhop) / 2);
    Q = mp / rhop;
    E = Ploss * Q / Asg; % Pumping power per unit area
    % (6) Get yearly cost during lifetime
    
    Cao = 945.86; % Cost per unit area(based on do) of reference SG
    Asgo = 4.645 * 10^3; % Heat exchanger area(based on do) of reference SG
    Cp = 5.33 * 10^(-9); % Cost per unit energy
    Ca = Cao * (Asg / Asgo)^(0.6); % Six-tenth rule, cost per unit area
    i = 0.05; % Interest rate
    n = 30; % SG lifetime
    R = i * (1 + i)^(n) / ((1 + i)^(n) - 1); % Capital recovery factor
    CF = 0.88; % Capacity factor
    theta = 365 * 24 * 3600 * CF; % Pump operation time
    
    C = Asg * (Ca * R + E * theta * Cp); % yearly cost during lifetime
    
    umatrix(j, 1) = up;
    Cmatrix(j, 1) = C;
    
    j = j + 1;
end
plot(umatrix, Cmatrix);
[mincost, I] = min(Cmatrix);
disp(mincost);
disp(umatrix(I, 1));

MATLAB: How to include nonexplicit constants in the solve function

b=solve('b-65','b');
b=subs(b)                    
eq1=['-es-2.69*Ta+' num2str(b)];
eq2='-es+6.11*exp((2.492E6/461.7)*(1/273-1/(Ta+273)))';
[Ta,es]=solve(eq1,eq2,'es','Ta')

Best Answer

Related Solutions

MATLAB: It gives Invalid use of operator error

MATLAB: How to include nonexplicit constants in the solve function

Related Question