Help with trapz(): Getting complex values

Question

Shashank 2012-12-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/57193-help-with-trapz-getting-complex-values

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