MATLAB Answers

Help with trapz(): Getting complex values

7 views (last 30 days)
Shashank
Shashank on 21 Dec 2012
Here is my code: I am getting complex values for "T_integral" and that is messing up my temperature profile TSil. Any help will be appreciated.
%%Initializing Variables
Lz = 0.00038; % thickness of the silicon wafer
Ft = 1; % simulation time
Nz = 10; % number of grid points
Nt = 100; % number of time steps
dz = Lz/Nz; % grid cell size
dt = Ft/Nt; % time step size
z = linspace(0,Lz,Nz+1);
%%Thermal Properties
alpha_Sil = 0.00008; % alpha for silicon (m^2/sec)
k_Sil = 149; % thermal conductivity for silicon (W/m K)
d = alpha_Sil*dt/(dz^2);
Ta = 298; % temperature of ambient air (K)
ha = 15; % heat transfer coefficient of air at Ta (W/m^2 K)
%%Optical Properties
a_Sil = 52.6; % absorption coefficient (m^-1)
r_SilAir = 0.34; % refelctivity of the Si-air interface (experimentally measured)
t_Sil = 0.9802; % transmissivity of the silicon layer due to absorption alone
e_Sil = 1 - r_SilAir; % emissivity of silicon surface
%%Constants
sigma = 5.67037E-08; % Stefan-Boltzmann constant
F_12 = 0.017517; % energy emitted by a black body contained within the two wavelength intervals of the IR camera
Ec = 10.3439; % energy measured by the IR camera
%%Finite Difference Scheme (Crank Nicholson in space & F.E in time)
TSil = 329*ones(1,Nz+1); % initial condition (arbitrary temperature profile)
TSilnew = TSil;
t = 0;
for kt = 1:Nt+1
t = t + dt;
for kz = 1:Nz+1
Y(kz) = (TSil(1,kz)^4).*exp(-a_Sil*z(kz));
end
T_integral = trapz(z,Y(:));
TSilnew(1) = (((Ec - (r_SilAir*F_12*sigma*(Ta^4)) - (e_Sil*a_Sil*sigma*F_12*T_integral))/(e_Sil*t_Sil*F_12*sigma)).^0.25);
for kz = 2:Nz
TSilnew(kz) = (d/(2*(1+d)))*(TSilnew(kz-1) + TSilnew(kz+1) + TSil(kz-1) + TSil(kz+1)) + ((1-d)/(1+d))*TSil(kz);
end
TSilnew(kz+1) = (d/(1+d-(dz*d*ha/k_Sil)))*(TSilnew(Nz) + TSil(Nz)) + ((1-d-(dz*d*ha/k_Sil))/(1+d-(dz*d*ha/k_Sil)))*TSil(Nz+1) - ((2*dz*ha*d/k_Sil)/(1+d-(dz*d*ha/k_Sil)))*Ta;
TSil = TSilnew;
end

Answers (2)

per isakson
per isakson on 21 Dec 2012
Edited: per isakson on 21 Dec 2012
  • Put the code in a function
  • Set a conditional, not(isreal(Y)), breakpoint on the line
T_integral = trapz(z,Y(:));
  • Run the function and you will see that eventually Y is complex
  • Set a conditional, not(isreal(TSilnew)), breakpoint on the line
TSilnew(1) = (((Ec - (r_SilAir*F_12*sigma*(Ta^4)).....
  • Run the function and you'll see that TSilnew(1) is set to a complex value for kt=56. You have a negative value raised to 1/4.

Shashank
Shashank on 21 Dec 2012
Per, How do I make that value real then?! And why would a simple integration result in a complex value at all?
  3 Comments
Shashank
Shashank on 21 Dec 2012
I realized that. Thanks for helping me out Per.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!