Heat Distribution in Circular Cylindrical Rod

Steady-State Solution

First, compute the steady-state solution. If the final time in the transient analysis is sufficiently large, the transient solution at the final time must be close to the steady state solution. By comparing these two results, you can check the accuracy of the transient analysis.

The 2-D model is a rectangular strip whose x-dimension extends from the axis of symmetry to the outer surface and y-dimension extends over the actual length of the rod (from -1.5 m to 1.5 m). Create the geometry by specifying the coordinates of its four corners.

g = decsg([3 4 0 0 .2 .2 -1.5 1.5 1.5 -1.5]');

Plot the geometry with the edge labels.

figure
pdegplot(g,EdgeLabels="on")
axis equal

Create an femodel for steady-state thermal analysis and include the geometry in the model.

model = femodel(AnalysisType="thermalSteady", ...
                Geometry=g);

Specify that the problem is axisymmetric.

model.PlanarType = "axisymmetric";

The rod is composed of a material with these thermal properties.

k = 40; % Thermal conductivity, W/(m*C)
rho = 7800; % Density, kg/m^3
cp = 500; % Specific heat, W*s/(kg*C)
q = 20000; % Heat source, W/m^3

For a steady-state analysis, specify the thermal conductivity of the material.

model.MaterialProperties = ...
    materialProperties(ThermalConductivity=k);

Specify the internal heat source.

model.FaceLoad = faceLoad(Heat=q);

Define the boundary conditions. There is no heat transferred in the direction normal to the axis of symmetry (edge 1). You can use the default boundary condition for this edge. Edge 2 is kept at a constant temperature T = 100 °C.

model.EdgeBC(2)=edgeBC(Temperature=100);

Specify the convection boundary condition on the outer boundary (edge 3). The surrounding temperature at the outer boundary is 100 °C, and the heat transfer coefficient is $50\mathrm{W}/\left(\mathrm{m}{\cdot }^{\circ }\mathrm{C}\right)$.

model.EdgeLoad(3) = ...
    edgeLoad(ConvectionCoefficient=50,...
             AmbientTemperature=100);

The heat flux at the bottom of the rod (edge 4) is $5000\mathrm{W}/{\mathrm{m}}^2$.

model.EdgeLoad(4) = edgeLoad(Heat=5000);

Generate and plot the mesh.

model = generateMesh(model);
figure
pdemesh(model)
axis equal

Solve the model and plot the result.

result = solve(model);
T = result.Temperature;
figure
pdeplot(result.Mesh,XYData=T,Contour="on")
axis equal
title("Steady-State Temperature")

Transient Solution

Switch the analysis type of the model to transient thermal analysis.

model.AnalysisType = "thermalTransient";

Specify the thermal conductivity, mass density, and specific heat of the material.

model.MaterialProperties = ...
    materialProperties(ThermalConductivity=k,...
                       MassDensity=rho,...
                       SpecificHeat=cp);

Specify that the initial temperature of 0 °C in the rod.

model.FaceIC = faceIC(Temperature=0);

Compute the transient solution for solution times from t = 0 to t = 50000 seconds.

tfinal = 50000;
tlist = 0:100:tfinal;
result = solve(model,tlist);

Plot the temperature distribution at t = 50000 seconds.

T = result.Temperature;

figure 
pdeplot(result.Mesh, ...
        XYData=T(:,end), ...
        Contour="on")
axis equal
title(sprintf(['Transient Temperature' ...
               ' at Final Time (%g seconds)'],tfinal))

Find the temperature at the bottom surface of the rod: first, at the center axis and then on the outer surface.

Tcenter = interpolateTemperature(result,[0.0;-1.5],1:numel(tlist));
Touter = interpolateTemperature(result,[0.2;-1.5],1:numel(tlist));

Plot the temperature at the left end of the rod as a function of time. The outer surface of the rod is exposed to the environment with a constant temperature of 100 °C. When the surface temperature of the rod is less than 100 °C, the environment heats the rod. The outer surface is slightly warmer than the inner axis. When the surface temperature is greater than 100 °C, the environment cools the rod. The outer surface becomes cooler than the interior of the rod.

figure
plot(tlist,Tcenter)
hold on
plot(tlist,Touter,"--")
title("Temperature at the Bottom as a Function of Time")
xlabel("Time, s")
ylabel("Temperature, C")
grid on
legend("Center Axis","Outer Surface","Location","SouthEast")