Problem Statement: Problem 2.14 of Riggs (heat conduction). The following equation, which governs heat conduction through a pipe wall, describes how temperature (T) changes with radius (r).
dT d2T
-- + r (---) = 0 Boundary Condition: T(rin)=Tin, T(rout)=Tout
dr dr2
where rin = inner radius
rout = outer radius
Tin = temperature at the inner radius
Tout = temperature at the outer radius
One numerical scheme is to divide the thickness of the wall into
N-1 equidistant intervals (or steps) and calculate the
temperature at N points. The finite difference approximations of
the first derivative dT/dr and the second derivative
d2T/dr2 are:
dT/dr ==> 0.5(Ti+1-Ti-1)/Dr d2T/dr2 ==> (Ti+1-2Ti+Ti-1)/(Dr)2where Dr is the step size.
For T=200°C at rin=0.05 and T=80°C at rout=0.2, set up the equations in the standard "A*x=b" form for N=4. Identify the "A" matrix and the "b" vector.
Solution:
Substituting the given derivative equivalence into the
differential equation gives the following general solution:
for i=1: T1=Tin
for i=N: TN=Tout
for i=2,3,...,N-1 (-0.5+r/dr)Ti-1 - (2r/dr)Ti + (0.5+r/dr)Ti+1 = 0
where r=rin+dr*(i-1)
Substituting the given boundary condition and the given number N=4, we have
dr = (0.2-0.05)/3 = 0.05
for i=1: r=0.05 T1=200
for i=2: r=0.10 (-0.5+0.10/0.05)T1 - (2*0.10/0.05)T2 + (0.5+0.10/0.05)T3 = 0
for i=3: r=0.15 (-0.5+0.15/0.05)T2 - (2*0.15/0.05)T3 + (0.5+0.15/0.05)T4 = 0
for i=4: r=0.20 T4= 80
Thus, we have A*x=b, where each vector/matrix is:
| 1 0 0 0 | | 200 | | T1 |
A=| 1.5 4 2.5 0 | b=| 0 | x=| T2 |
| 0 2.5 6 3.5 | | 0 | | T3 |
| 0 0 0 1 | | 80 | | T4 |
Write a Matlab function cond_fcn that has as its input the boundary conditions, i.e., the following five input variables: rin, rout, Tin, Tout, and N. The function should return two vectors: a column of r and a column of T that corresponds to temperature at different points in the pipe wall.
Solution:
function [r, T] = T (r_in, r_out, T_in, T_out, n)
%-----------------------------------------------------------------------
% Return the temperature profile across the pipe wall due to heat conduction.
%-----------------------------------------------------------------------
% Declare variables ----------------------------------------------------
r = zeros(n+1,1);
A = zeros(n+1,n+1);
b = zeros(n+1,1);
% Assign constants in the linear equation ------------------------------
dr = (r_out-r_in)/n;
for i=2: n
r(i) = r_in + dr*(i-1);
A(i,i-1) = -1./2. + r(i)/dr;
A(i,i) = -2.*r(i)/dr;
A(i,i+1) = 1./2. + r(i)/dr;
end
% Boundary Conditions --------------------------------------------------
r(1) = r_in;
A(1,1) = 1.;
b(1) = T_in;
r(n+1) = r_out;
A(n+1,n+1) = 1.;
b(n+1) = T_out;
% Solve the linear set of equations ------------------------------------
T = A\b;
Write a brief main MATLAB script (program) that calls the above function and plots temperature vs radius. Assume the above function is stored in a file named "cond_fcn.m"
Solution:
At the minimum, the following two lines that are not commented out will do.
More detailed program on heat conduction is here.
% Specify the boundary condition ---------------------------------------
% r_in = 0.05; T_in = 200; % Condition at the inner radius
% r_out = 0.20; T_out= 80; % Condition at the outer radius
% n = 101; % number of points
% Call a routine to find the temperature profile -----------------------
% Usage: [r T]=cond_fcn(r_in, r_out, T_in, T_out, n);
[r T]=cond_fcn(0.05, 0.20, 200, 80, 101);
% Output the answer ----------------------------------------------------
plot(r,T)
Problem Statement: Write a FORTRAN program. You may call a subroutine from a FORTRAN library of your choice, including the Thomas method from Riggs book. Or you may write a program to output the matrix coefficients specific to this problem and pipe the matrix to the general program you wrote in the last problem. With the latter approach, you only have to concentrate on putting a given problem in the standard "Ax=b" form, then make use of the general linear equation solver you have written earlier.
Solution:
Problem Statement: Re-work with MATLAB. Plot the temperature profile.
Solution:
Solution:
Solution:
|
