Answer Happy

Convert Pascal Coding to Matlab
Pascal Code:
Program Prob3_M; {solution to Problem 3_M}
uses crt, printer;
const N = 11; {radial nodes}
M = 21; {axial nodes}
k = 40.0; {thermal conductivity (W/mK)}
alpha = 1e–5; {thermal diffusivity (m^2/s)}
Rs = 0.025; {disk radius (m)}
Zs = 0.005; {disk thickness (m)}
qmax = 3.0e6; {peak flux (W/m^2)}
Ro = 0.05; {parameter in flux equation (m)}
Tinit = 20.0; {initial temperature (C)}
Tmax = 300.0; {maximum desired temperature at top of axis
(C)}
var dr, dz, rhoC, dtMax, dt, time:real;
i, j : integer;
C,R,L,U,D : array [1..N,1..M] of real;
q,V,Aco,Aci,Af : array [1..N] of real;
Toid,Tnew : array [1..N+1,1..M+1] of real;
begin
{calculate size of the control volumes}
dr : = Rs/(N – 1);
dz : = Zs/(M – 1);
rhoC: = k/alpha;
{calculate control volume surface areas and volumes}
Af[1] : = pi*dr*dr/4.0;
Af[N] : = pi*dr*(Rs – dr/4.0);
for i : = 2 to N – 1 do Af : = 2.0*pi*dr*dr*(i – 1);
for i : = 1 to N do V : = Af*dz;
Aco[N] : = 2.0*pi*Rs*dz;
for i : = 1 to N – 1 do Aco : = 2.0*pi*dr*dz*(i – 0.5);
Aci[1] : = 0.0;
for i : = 2 to N do Aci : = 2.0*pi*dr*dz*(i – 1.5);
{calculate absorbed flux as function of i}
q [1] : = pi/4.0*qmax*dr*dr*(1.0 – 0.9/8.0*sqr(dr/Ro));
q [N] : = pi*qmax*dr*dr*(N – 1.25–0.9/8.0*sqr(dr/Ro)*(2.0*N*N*N
–7.5*N*N
+ 9.5*N – 65.0/16.0));
for i : = 2 to N – 1 do
q : = pi*qmax*dr*dr*(2.0*(1 – 1.0) –
0.9/4.0*sqr(dr/Ro)*(4.0*i*i*i – 12.0*i*
i + 13.0*i – 5.0));
{fill in the coefficient matrices}
for i : = 1 to N do {first, zero all of them out}
for j : = 1 to M do
begin
R[i, j] : = 0.0;
L[i, j] : = 0.0;
U[i, j] : = 0.0;
D[i, j] : = 0.0;
C[i, j] : = 0.0;
end
for i : = 2 to N – 1 do
for j : = 2 to M – 1 do
begin
R[i, j] : = Aco/dr;
L[i, j] : = Aci/dr;
U[i, j] : = Af/dz;
D[i, j] : = Af/dz;
end;
i : = 1;
for j : = 2 to M – 1 do
begin
R[i, j] : = Aco[i]/dr;
U[i, j] : = Af[i]/dz;
D[i, j] : = Af[i]/dz;
end;
i : = 1;
j : = M;
R[i, j] : = Aco[i]/dr;
D[i, j] : = Af[i]/dz;
i : = 1;
j : = 1;
R[i, j] : = Aco[i]/dr;
U[i, j] : = Af[i]/dz;
i : = N;
for j : = 2 to M – 1 do
begin
L[i, j] : = Aci[i]/dr;
U[i, j] : = Af[i]/dz;
D[i, j] : = Af[i]/dz;
end;
i : = N;
j : = M;
L[i, j] : = Aci[i]/dr;
D[i, j] : = Af[i]/dz;
i : = N;
j : = 1;
L[i, j] : = Aci[i]/dr;
U[i, j] : = Af[i]/dz;
j : = M;
for I : = 2 to N – 1 do
begin
R[i, j] : = Aco[i]/dr;
L[i, j] : = Aci[i]/dr;
D [i, j] : = Af[i]/dz;
end;
j : = 1;
for I : = 2 to N – 1 do
begin
R [i, j] : = Aco[i]/dr;
L [i, j] : = Aci[i]/dr;
U [i, j] : = Af[i]/dz;
end;
{find maximum permissible dt}
dtMax : = 0.0;
for i : = 1 to N do
for j : = 1 to M do
begin
dt : = V[i]/alpha/(R[i, j] + L[i, j] + U[i, j] + D[i, j]);
if = dt > dtMax then dtMax : = dt;
end;
dt : = 0.5*dtMax; {actual value to be used in solution}
{fill in the cij matrix}
for i : = 1 to N do C[i, M] : = dt*q [i]/rhoC/V[i];
{establish the initial conditions}
for i : = 1 to N do
for j : = 1 to M do
Told [i, j] : = Tinit;
{carry out the solution}
time : = 0.0;
repeat
time: = time + dt;
writeln (time: 10:5);
for i: = 1 to N do
for j: = 1 to M do
Tnew [i,j] : = Told [i,j]*(1.0 – alpha*dt/V[i]*(R[i,j]+ L[i,j]+
U[i,j]+
D[i,j))
+ alpha*dt/V[i]*(R[i, j]*Told [i + 1, j] + L[i, j]*Told [i – 1, j]
+ U[i, j]*Told
[i, j + 1] + D[i, j]*Told [i, j – 1]) + C[i, j];
if Tnew[1, m] > Tmax then {print out distribution and
quit}
begin
writeln(1st, time : 8 : 4, ‘sec dt = ‘, dt : 15 : 10);
write(1st,’ ’);
for i : = 1 to N do write (1st, 1 : 10);
writeln(1st);
for j : = M downto 1 do
begin
write(1st, j : 4);
for i : = 1 to N do write(1st, Tnew[i, j]; 10 : 5);
writeln(1st);
end;
writeln(1st);
halt;
end
(otherwise, keep going}
for i : = 1 to N do
for j : = 1 to M do
Told[i, j] ; = Tnew[i, j];
until time < - 1.0;
end.