lda/nmu-m
2
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

misc fix

Browse Source
This commit is contained in:
lda committed 2025-06-26 17:22:10 +00:00
1 parent 7daf75bbb3
commit 3db31075e3
4 files changed
+285 -3
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

proj/proj.m
+278
View File
@@ -0,0 +1,278 @@
%% Project: 1D Heat Transfer Simulation (Numerical Solution of Heat Equation)
% This Live Script provides a complete solution to the MATLAB Project focusing
% on determining heat transfer through a 1D bar by conduction using the explicit
% finite difference method. It addresses all questions from the project document
% (Questions 1-7).
%% Introduction to the Heat Equation
% The general 3D Heat Equation describes temperature distribution over time
% and space:
%
% $$\rho c_p \frac{\partial T}{\partial t} = \lambda \left( \frac{\partial^2
% T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial
% z^2} \right)$$
%
% Where:
%
% - $T$: Temperature ($^{\circ}C$)
%
% - $\rho$: Mass density ($kg/m^3$)
%
% - $c_p$: Specific heat capacity ($J/(kg \cdot ^{\circ}C)$)
%
% - $\lambda$: Thermal conductivity ($W/(m \cdot ^{\circ}C)$)
%
% - $x, y, z$: Spatial dimensions (m) - $t$: Time (s)
%
% For a 1D bar, assuming homogeneous and isotropic material and insulated sides
% (except for the ends), the equation simplifies to:
%
% $$\rho c_p \frac{\partial T}{\partial t} = \lambda \frac{\partial^2 T}{\partial
% x^2}$$
%
% This can be rewritten using thermal diffusivity $\alpha = \frac{\lambda}{\rho
% c_p}$:
%
% $$\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}$$
%% Define Physical Parameters and Material Properties
% Bar Dimensions and Discretization
L = 0.4; % Length of the bar (m)
dx = 0.01; % Spatial step size (m)
dt = 0.5; % Time step size (s)
% Thermal Conditions
Tg = 100; % Temperature at left end (x=0) in Celsius
Td = 50; % Temperature at right end (x=L) in Celsius
T0 = 0; % Initial temperature of the bar (internal nodes) in Celsius
% Material Properties (Aluminum properties as hinted in the document)
rho = 2700; % Density (kg/m^3)
Cp = 900; % Specific heat capacity (J/kg/C)
lambda = 237; % Thermal conductivity (W/m/C)
% Simulation End Time (chosen to observe steady-state behavior)
t_final_s_sim = 2000; % Total simulation time in seconds
% Time points at which to capture and plot full temperature profiles
% This allows for a denser visualization early on, and sparser later as equilibrium is approached.
specific_early_points = [0, 5, 10, 20];
mid_range_points = 50:50:750;
late_range_points = 750:250:t_final_s_sim;
plot_time_points = unique([specific_early_points, mid_range_points, late_range_points]);
%% Question 1: Express the left term of the equation (*) using numerical differentiation.
% The left term is the first partial derivative of temperature with respect
% to time, $\frac{\partial T}{\partial t}$. Using a forward finite difference
% approximation for the time derivative at node $i$ and time step $n$:
%
% $$\frac{\partial T}{\partial t} \approx \frac{T_i^{n+1} - T_i^n}{\Delta t}$$
%
% This approximates the rate of temperature change at a specific spatial point
% $i$ from time $t_n$ to $t_{n+1}$.
%% Question 2: Express the right term of the equation (*) using numerical differentiation.
% The right term is the second partial derivative of temperature with respect
% to space, $\frac{\partial^2 T}{\partial x^2}$.
%
% Using a central finite difference approximation for the second spatial derivative
% at node $i$ and time step $n$:
%
% $$\frac{\partial^2 T}{\partial x^2} \approx \frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{(\Delta
% x)^2}$$
%
% This relates the curvature of the temperature profile at node $i$ to the temperatures
% of its immediate neighbors at the current time step.
%% Question 3: Reformulate the equation (*) using the results of Question 1 and 2.
% Substitute the finite difference approximations into the 1D Heat Equation:
%
% $$\frac{T_i^{n+1} - T_i^n}{\Delta t} = \alpha \frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{(\Delta
% x)^2}$$
%
% Where $\alpha = \frac{\lambda}{\rho c_p}$ is the thermal diffusivity.
%% Question 4: From the new formula of the equation, express $T_i^{n+1}$.
% Rearranging the reformulated equation to explicitly solve for the temperature
% at the next time step $T_i^{n+1}$:
%
% $$T_i^{n+1} = T_i^n + \alpha \frac{\Delta t}{(\Delta x)^2} (T_{i+1}^n - 2T_i^n
% + T_{i-1}^n)$$
%
% This is the core update equation for each interior node.
%
% For convenience, we define the Fourier number (or numerical stability parameter)
% $r$:
%
% $$r = \alpha \frac{\Delta t}{(\Delta x)^2}$$
%
% So the update equation becomes:
%
% $$T_i^{n+1} = T_i^n + r (T_{i+1}^n - 2T_i^n + T_{i-1}^n)$$
%% Question 5: What are the boundary conditions and initial temperature?
% a. Temperature on the left $x=0$, $T_g$?
%
% The temperature at the left end of the bar ($x=0$) is a fixed boundary condition:
% $T_g = 100^{\circ}C$.
%
% This means $T_1^n = 100^{\circ}C$ for all time steps $n$.
%
% b. Temperature on the right $x=L$, $T_d$?
%
% The temperature at the right end of the bar ($x=L$) is a fixed boundary condition:
% $T_d = 50^{\circ}C$.
%
% This means $T_{N_x}^n = 50^{\circ}C$ for all time steps $n$.
%
% c. Initial temperature $T_0$?
%
% The initial temperature of the *internal* nodes of the bar (at time $t=0$)
% is $T_0 = 0^{\circ}C$.
%
% This means $T_i^0 = 0^{\circ}C$ for all internal nodes $i$.
%% Question 6: Write the algorithm to calculate all temperatures at each point of the bar and at any time.
% This section outlines the step-by-step algorithm implemented in MATLAB.
%
% 1. **Calculate Derived Constants:** Compute thermal diffusivity (alpha) and
% the Fourier number (r). Check the stability criterion: r must be <= 0.5 for
% the explicit method to be stable.
%
% 2. **Spatial Grid Setup:** Create a vector 'x' representing the discrete
% spatial nodes along the bar. Determine the total number of spatial nodes (Nx).
%
% 3. **Initialize Temperature Array:** Create 'T_current' (a 1D array/vector)
% to hold temperatures at the current time step. Set all interior nodes to 'T0'.
% Apply fixed boundary conditions 'Tg' and 'Td' to the first and last nodes.
%
% 4. **Time Loop:** Iterate through discrete time steps from t=0 up to 't_final_s_sim'.
%
% a. For each time step, create a 'T_next' array (initialized to 'T_current').
%
% b. **Apply Finite Difference Equation: For each *interior* node (from index
% 2 to Nx-1), calculate its temperature at the next time step using the derived
% formula from Question 4.
%
% c. **Update Boundary Conditions:** The boundary node temperatures (T_next(1)
% and T_next(Nx)) remain fixed at Tg and Td, respectively. These are carried over
% from T_current or explicitly set to override any updates to ensure they remain
% fixed.
%
% d. **Advance Time Step:** Set 'T_current = T_next' to prepare for the next
% iteration.
%
% e. **Data Storage for Plotting:** Store the full temperature profile (T_current)
% at predefined 'plot_time_points'. Also, record the temperature at the midpoint
% of the bar over all time steps for a separate time-series plot.
%% Question 7: Use MATLAB to solve the equation.
% This section provides the MATLAB implementation based on the algorithm and detailed parameters.
fprintf('--- Solving 1D Heat Equation Project ---\n');
% Calculate Derived Constants
alpha = lambda / (rho * Cp);
r = alpha * dt / (dx^2);
fprintf('Thermal Diffusivity (alpha): %.4e m^2/s\n', alpha);
fprintf('Fourier Number (r): %.4f\n', r);
if r > 0.5
warning('The explicit method is unstable for r > 0.5. Consider reducing dt or increasing dx.');
end
% Spatial Grid Setup
x = 0:dx:L; % Spatial grid points
Nx = length(x); % Number of spatial nodes
% Find the index for the midpoint of the bar for plotting the right pane
mid_index = find(abs(x - L/2) < 1e-9, 1); % Find index closest to L/2
if isempty(mid_index) % Fallback if L/2 isn't an exact node
mid_index = round(Nx / 2);
end
% Initialize Temperature Array
T_current = zeros(1, Nx); % Initialize temperature array with zeros
T_current(2:Nx-1) = T0; % Internal nodes at initial temperature (0 C)
T_current(1) = Tg; % Left boundary condition (100 C)
T_current(Nx) = Td; % Right boundary condition (50 C)
% Time Grid Setup
time_vector = 0:dt:t_final_s_sim;
Nt = length(time_vector); % Number of time steps
% Storage for results
T_midpoint_over_time = zeros(1, Nt); % Stores temp at midpoint over ALL time steps
T_midpoint_over_time(1) = T_current(mid_index); % Initial midpoint temp
all_profiles = cell(size(plot_time_points)); % Stores full profiles at specific times
all_profile_times = cell(size(plot_time_points)); % Stores labels for these profiles
% Capture initial profile (t=0) if it's one of the requested plot_time_points
[is_t0_plot_point, idx_t0_plot_point] = ismembertol(0, plot_time_points, dt/2, 'DataScale', 1);
if is_t0_plot_point
all_profiles{idx_t0_plot_point} = T_current;
all_profile_times{idx_t0_plot_point} = sprintf('t = %.0f s', 0);
end
% Main Time-Stepping Loop
fprintf('\nStarting simulation...\n');
for n = 1:Nt-1 % Loop from t=0 (n=1) to t_final_s (n=Nt)
T_next = T_current; % Initialize T_next. Boundary values from T_current are implicitly carried over.
% Apply finite difference equation for internal nodes (from index 2 to Nx-1)
for i = 2:Nx-1
T_next(i) = T_current(i) + r * (T_current(i+1) - 2*T_current(i) + T_current(i-1));
end
T_current = T_next; % Update temperature for the next iteration
% Store midpoint temperature for time series plot
T_midpoint_over_time(n+1) = T_current(mid_index);
% Store full spatial profiles for plotting at specified plot_time_points
[is_plot_time, idx_to_store] = ismembertol(time_vector(n+1), plot_time_points, dt/2, 'DataScale', 1);
if is_plot_time
all_profiles{idx_to_store} = T_current;
all_profile_times{idx_to_store} = sprintf('t = %.0f s', time_vector(n+1));
end
end
fprintf('Simulation finished.\n');
%% Visualization (Question 7 Continued)
% Create a figure to display the results in two subplots
figure('Color', 'white', 'Position', [100, 100, 1000, 400]); % Adjust figure size for better viewing
% Left Pane: Temperature Profiles Over Time (T vs x)
subplot(1, 2, 1);
hold on;
grid on;
box on; % Add a box around the plot
% Plot each captured profile
colors = lines(length(plot_time_points)); % Generate distinct colors based on number of plot points
for k = 1:length(plot_time_points)
if ~isempty(all_profiles{k}) % Ensure data exists for this plot point
plot(x, all_profiles{k}, 'DisplayName', all_profile_times{k}, 'LineWidth', 1.5, 'Color', colors(k,:));
end
end
xlabel('Position along bar, x (m)');
ylabel('Temperature, T (^oC)');
title('Temperature Profiles Over Time');
% Conditionally display legend if not too many lines
if length(plot_time_points) < 8 % Adjust threshold as needed
legend('Location', 'southeast');
end
hold off;
% Right Pane: Temperature at Midpoint Over Time (T vs t)
subplot(1, 2, 2);
plot(time_vector, T_midpoint_over_time, 'LineWidth', 2, 'Color', [0.8500 0.3250 0.0980]);
grid on;
box on; % Add a box around the plot
xlabel('Time, t (s)');
ylabel('Temperature at Midpoint (^oC)');
title('Temperature at Midpoint (x=L/2) Over Time');
ylim([0 100]); % Set y-limits to reasonable temperature range
xlim([0 t_final_s_sim]); % Set x-limits to the simulation time
% Final Output
fprintf('\nSimulation and plotting complete. Check the generated figure for results.\n');
% To save the figure, use:
% saveas(gcf, 'Figure_1.png');
% (Uncomment and run the above line after the figure appears to save it.)
proj/proj.mlx
BIN
View File
Binary file not shown.
proj/runHeatProject.m
+1 -1
View File
@@ -23,4 +23,4 @@ plotProfiles(x_coords, all_temp_profiles, profile_labels);
subplot(1, 2, 2);
plotMidpointTemperature(time_vec, midpoint_T_series);
exportgraphics(f, 'Figure_1.pdf', 'ContentType', 'vector','Resolution', 150, 'Units', 'centimeters', Width=1920);
exportgraphics(f, 'Figure_1.pdf', 'ContentType', 'vector');
proj/solveHeatEquation.m
+6 -2
View File
@@ -23,9 +23,13 @@ function [x, T_final_profile, T_time_series, time_vector, all_profiles, all_prof
Tg = thermal_cond.Tg; Td = thermal_cond.Td; T0 = thermal_cond.T0;
dx = steps.dx; dt = steps.dt;
% Calculate alpha and r (Fourier Number)
% Calculate Thermal Diffusivity (alpha)
alpha = lambda / (rho * Cp);
fprintf('Thermal Diffusivity (alpha): %.4e m^2/s\n', alpha);
% Calculate the Fourier Number (r) and check stability
r = alpha * dt / (dx^2);
fprintf('Fourier Number (r): %.4f\n', r);
if r > 0.5
warning('The explicit method is unstable for r > 0.5. Consider reducing dt or increasing dx.');
@@ -53,7 +57,7 @@ function [x, T_final_profile, T_time_series, time_vector, all_profiles, all_prof
all_profiles = cell(length(plot_time_points), 1);
all_profile_times = cell(length(plot_time_points), 1);
% Capture initial profile (t=0)
% Capture initial profile (t=0) if it's one of the requested plot_time_points
if ismembertol(0, plot_time_points, dt/2, 'DataScale', 1)
all_profiles{1} = T_current;
all_profile_times{1} = sprintf('t = %.0f s', 0); % Using sprintf for consistent formatting