Explicit Heat Conduction in 3d

L=1;
t=1;
k=.001;
nx=11;
ny=11;
nz=11;
nt=500;
dx=L/n;
dt=.002;
alpha=k*dt/dx^2;
T0=400*ones(1,n);
T1=300*ones(1,n);
T0(1) = 300;
T0(end) = 300;
for j=2:ny-1
    for i=2:nx-1
        for k=2:nz-1
            for j=2
                T1(i,j,k)=alpha*dt*((T0(i+1,j,k)-2*T0(i,j,k)+T0(i-1,j,k))+(T0(i,j+1,k)-2*T0(i,j,k)+T0(i,j-1,k))+(T0(i,j,k+1)-2*T0(i,j,k)+T0(i,j,k-1))
            end
            T0=T1;
        end
        plot(T1)

[Walter Roberson]

What is the question?

[Yichen Yang]

Consider a typical egg with approximately its semi major(a) and semi minor(b) axes as 4.0 x 3.0 x 3.0 cm, and consider it as a ellipsoid whose equation is given as (x/a)**2 + (y/b)**2 + (z/b) **2 = 1. Discretize the rectangular parallelopiped of size 2a * 2b* 2b by a uniform Cartesian grid of size 40 x 30 x 30 discrete elements. Blank and skip the cells outside the ellipsoid. Now, boil the egg by starting the initial temperature as 300K, and impose a temperature of 350K on the surface. Go to any search engine and find the specific heat, conductivity, and density of the egg yoke. Consider the yellow and white to have the same properties and calculate the evolution of the temperature. Plot temperature contours for a plane that cuts the center of the egg, along the major and one of the minor axes (i. e. x and y coordinates; the egg is a body of revolution). Calculate the time it will take the center (x =0, y=0, z=0) (where egg goes from -a to a, -b to b and -b to b to reach a temperature of 330 K. When can you consider an egg to be soft boiled? And when can you consider it a hard boiled? Look at some literature.

[Yichen Yang]

alpha=1e-4;

t=1600;

nt=160;

xelta_t=t/nt;

xlength=1;

nx=15;

delta_x=xlength/nx;

T1=300;

T2=300;

T3=300;

T4=300;

T5=350;

Tmax=max([T1,T2,T3,T4,T5]);

n=((xlength/delta_x)+1)^2;

m=sqrt(n);

r=t/delta_t+1;

d=alpha*delta_t/delta_x^2;

T=zeros(m,m,r)

for k=1:m

for i=1:m

for j=1:m

if (i==1)

if(i==1)&&(j==1)

T(i,j,k)=(T4+T1)/2;

else if(i==1&&(j>1&&j<m))

T(i,j,k)=T1;

elsei f(i==m&&(j>1&&j<m))

T(i,j,k)=T3;

elseif(j==1&&(i>1&&i<m))

T(i,j,k)=T4;

elseif(j==m&&(i>1&&i<m))

T(i,j,k)=T2;

else

T(i,j,k)=T5;

elseif(i==1)&&(j==m)

T(i,j,k)=(T4+T2)/2;

elseif(i==m)&&(j==m)

T(i,j,k)=(T3+T2)/2;

elseif(i==m)&&(j==1)

T(i,j,k)=(T3+T4)/2;

end

end

end

end

for k=1:r-1

for i=2:m-1

for j=2:m-1

T(i,j,k+1)=T(i,j,k)+d*(T(i+1,j,k)-2*T(i,j,k)+T(i-1,j,k))+(T(i,j+1,k)-2*T(i,j,k)+T(i,j-1,k))+(T(i,j,k+1)-2*T(i,j,k)+T(i,j,k-1))

end

end

end

T(:,:,1);

T(:,:,r);

T;

[Yichen Yang]

this is what I have now

[Yichen Yang]

can u also show me how to plot?

Answer 1

[William Rose]

@Yichen Yang,

Please format your code as code in the window. This will make it easier for others to read and understand. You will also be able to run your code in the window, and the results (plots, error messages, etc.) will appear in the window. Then other people will know what happens when you run your code.

I will do it below.

alpha=1e-4;
t=1600;
nt=160;
xelta_t=t/nt; %This looks like an error. Do you mean delta_t=?
xlength=1;
nx=15;
delta_x=xlength/nx;
T1=300;
T2=300;
T3=300;
T4=300;
T5=350;
Tmax=max([T1,T2,T3,T4,T5]);
n=((xlength/delta_x)+1)^2;
m=sqrt(n);
r=t/delta_t+1;
d=alpha*delta_t/delta_x^2;
T=zeros(m,m,r)
for k=1:m
    for i=1:m
        for j=1:m
            if (i==1)
                if(i==1)&&(j==1)
                    T(i,j,k)=(T4+T1)/2;
                else if(i==1&&(j>1&&j<m))
                        T(i,j,k)=T1;
                elseif (i==m&&(j>1&&j<m))   % I corrected an error. It was "elsei f(..."
                        T(i,j,k)=T3;
                elseif(j==1&&(i>1&&i<m))
                    T(i,j,k)=T4;
                elseif(j==m&&(i>1&&i<m))
                    T(i,j,k)=T2;
                else
                    T(i,j,k)=T5;
                elseif(i==1)&&(j==m)
                    T(i,j,k)=(T4+T2)/2;
                elseif(i==m)&&(j==m)
                    T(i,j,k)=(T3+T2)/2;
                elseif(i==m)&&(j==1)
                    T(i,j,k)=(T3+T4)/2;
                	
                end
                end
            end
        end
        for k=1:r-1
            for i=2:m-1
                for j=2:m-1
                    T(i,j,k+1)=T(i,j,k)+d*(T(i+1,j,k)-2*T(i,j,k)+T(i-1,j,k))+(T(i,j+1,k)-2*T(i,j,k)+T(i,j-1,k))+(T(i,j,k+1)-2*T(i,j,k)+T(i,j,k-1))
                end
            end
        end
        T(:,:,1);
        T(:,:,r);
        T;

OK, now I see what happens when it runs. You have some mixed up else and if and elseif statements, and for statements without matching end statements. I corrected one typographical error where you had "elsei f(...", but this was not enough to fix the problems. Please check your code carefully to reilve the unmatched for and end statments. Make sure all you if, elseif, else statements are organized correctly. In one place you have

            if (i==1)
                if(i==1)&&(j==1)
                    T(i,j,k)=(T4+T1)/2;
                else if(i==1&&(j>1&&j<m)) %did you intend "elseif"?
                        T(i,j,k)=T1;
                elseif (i==m&&(j>1&&j<m))
                        T(i,j,k)=T3;
                elseif(j==1&&(i>1&&i<m))
                    T(i,j,k)=T4;
                elseif(j==m&&(i>1&&i<m))
                    T(i,j,k)=T2;
                else
                    T(i,j,k)=T5;
                elseif(i==1)&&(j==m)...  %This generates an error. You may not use elseif after else.

Fix these issues first.

I suggest you create an array T=temperature with four dimensions, for x, y, z, and time. Your array T has only three dimensions.

Initialize the system with 3 nested for loops over x,y,z. On the inside of the innermost loop, you have an if statement to decide if the current location is inside or outside. If inside, let T=300, otherwise, T=350.

Now run the simulation by using four nested for loops. The outer for loop is over time. The inner loops are over x, y, z. At the center, evaluate whether the current spatial location is inside the egg or not. If the locaiton is outside, let T=350. Otherwise, update the temperature using the discretized 3-D diffusion equation. It will be simpler than the code you have shown above.

[Yichen Yang]

can u also show me how to plot? thanks

Answer 2

[William Rose]

@Yichen Yang,

To plot results from your 4-D array T:

a=2.0; b=1.5; c=1.5;            %egg semi-diameters along x,y,z axes [cm]
dx=0.1; dy=dx; dz=dx;           %voxel dimensions [cm]
nx=42; ny=32; nz=32;            %number of voxels along each axis
%Number of voxels is chosen to assure there is at least one voxel outside every egg voxel.
%This assures that the discrete diffusion equation will not go out of bounds,
%for all voxels inside the egg.
x=(1:nx)*dx-(nx+1)*dx/2;        %x-coordinates of the voxels
y=(1:ny)*dy-(ny+1)*dy/2;        %y-coordinates of the voxels
z=(1:nz)*dz-(nz+1)*dz/2;        %z-coordinates of the voxels

dt=1.0;                         %time step [s]
tmax=600;                       %max time [s]
nt=tmax/dt+1;                   %number of time steps
t=(0:nt-1)*dt;                  %time vector

T0=300;                         %initial inside temperature [K]
Text=350;                       %outside temperature [K]
T=Text*ones(nt,nx,ny,nz);       %allocate array T=temperature

Add your own code to compute the correct values for T(m,i,j,k) at the start and at all subsequent times and locations.

Then you can plot the results:

%Plot temperature on slice though the center, at time t=0
[X,Y]=meshgrid(x,y);
figure; surf(X,Y,squeeze(T(1,:,:,nz/2))','FaceColor','interp','EdgeColor','none');
c=colorbar; c.Label.String='Temp. (\circK)';
view(0,90); axis equal; hold on;
plot3(a*cos(2*pi*(0:100)/100),b*sin(2*pi*(0:100)/100),Text*ones(1,101),'-k');
title('Temperature at t=0');
%Plot temperature on slice though the center, at halfway time
figure; surf(X,Y,squeeze(T((nt+1)/2,:,:,nz/2))','FaceColor','interp','EdgeColor','none');
caxis([T0,Text]); c=colorbar; c.Label.String='Temp. (\circK)';
view(0,90); axis equal; hold on;
plot3(a*cos(2*pi*(0:100)/100),b*sin(2*pi*(0:100)/100),Text*ones(1,101),'-k');
titlestr=sprintf('Temperature at t=%.0f',t((nt+1)/2)); title(titlestr);
%Plot temperature on slice though the center, at time t=tmax
figure; surf(X,Y,squeeze(T(end,:,:,nz/2))','FaceColor','interp','EdgeColor','none');
caxis([T0,Text]); c=colorbar; c.Label.String='Temp. (\circK)';
view(0,90); axis equal; hold on;
plot3(a*cos(2*pi*(0:100)/100),b*sin(2*pi*(0:100)/100),Text*ones(1,101),'-k');
titlestr=sprintf('Temperature at t=%.0f',t(end)); title(titlestr);
%Plot temperature versus time at edge, halfway in, and center.
figure; plot(t,T(:,2,ny/2,nz/2),'-rx',t,T(:,round(nx/4),ny/2,nz/2),'-go',t,T(:,nx/2,ny/2,nz/2),'-b+');
xlabel('Time (s)'); ylabel('Temperature (\circK)'); grid on
title('Temperature versus Time: Edge, Halfway, Center'); ylim([T0 Text]);
ls1=sprintf('x=%.2f',x(2)); ls2=sprintf('x=%.2f',x(round(nx/4))); 
ls3=sprintf('x=%.2f',x(nx/2)); legend(ls1,ls2,ls3);

Two of the four plots are shown below.