init

lab6/lab6b_me.m

@@ -0,0 +1,3 @@
+%% still copying this from somewhere; not making it myself AHHhhhh
+[t, y] = ode45(@(t, y) 2*t, [0 5], 0)
+plot(t, y)

proj/plotMidpointTemperature.m

@@ -0,0 +1,16 @@
+function plotMidpointTemperature(time_series, temp_at_midpoint)
+% PLOTMIDPOINTTEMPERATURE Plots temperature at the midpoint over time.
+%
+% Inputs:
+%   time_series     - Time points corresponding to temp_at_midpoint (1 x Nt vector)
+%   temp_at_midpoint - Temperature values at the midpoint (1 x Nt vector)
+
+    plot(time_series, temp_at_midpoint, 'LineWidth', 2, 'Color', [0.8500 0.3250 0.0980]);
+    grid on;
+    box on;
+    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 max(time_series)]); % Set x-limits to the simulation time
+end

proj/plotProfiles.m

@@ -0,0 +1,26 @@
+function plotProfiles(x_coords, profiles_data, labels)
+% PLOTPROFILES Plots temperature profiles at different time points.
+%
+% Inputs:
+%   x_coords    - Spatial grid points (Nx x 1 vector)
+%   profiles_data - Cell array of temperature profiles (each cell is 1 x Nx vector)
+%   labels      - Cell array of strings for plot legends (e.g., {'t=0s', 't=10s'})
+
+    hold on;
+    grid on;
+    box on;
+    
+    colors = lines(length(profiles_data)); % Generate distinct colors
+    for k = 1:length(profiles_data)
+        if ~isempty(profiles_data{k}) % Ensure the profile data exists
+            plot(x_coords, profiles_data{k}, 'DisplayName', labels{k}, 'LineWidth', 1.5, 'Color', colors(k,:));
+        end
+    end
+    
+    xlabel('Position along bar, x (m)');
+    ylabel('Temperature, T (^oC)');
+    title('Temperature Profiles Over Time');
+    legend('Location', 'southeast');
+    ylim([0 100]); % Set y-limits to reasonable temperature range
+    hold off;
+end

proj/proj.md

@@ -0,0 +1,26 @@
+# Project
+
+## **Q1: Time Derivative ($\frac{\partial T}{\partial t}$)**
+
+- This is approximated using a **forward finite difference** formula: $$\frac{T_i^{n+1} - T_i^n}{\Delta t}$$
+
+## **Q2: Second Spatial Derivative ($\frac{\partial^2 T}{\partial x^2}$)**
+
+- This is approximated using a **central finite difference** formula: $$\frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{(\Delta x)^2}$$
+
+## **Q3: Reformulate the Equation**
+
+- Substitute the finite difference approximations into the 1D Heat Equation $ \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} $:
+        $$\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}$$
+
+## **Q4: Express $ T_i^{n+1} $**
+
+- Rearranging the reformulated equation to explicitly solve for temperature at the next time step:
+  $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 can be simplified by defining $r = \alpha \frac{\Delta t}{(\Delta x)^2}$:
+  $T_i^{n+1} = T_i^n + r (T_{i+1}^n - 2T_i^n + T_{i-1}^n)$.
+
+## **Q5: Boundary and Initial Conditions**
+
+- **Initial Condition ($T_0$):** Defines the temperature of the _entire bar (internal nodes)_ at the starting time ($t=0$). For question 7 of this project, it's $0^{\circ}C$.
+- **Boundary Conditions ($T_g$, $T_d$):** Define the fixed temperatures at the _ends of the bar_ ($x=0$ and $x=L$) for all times ($t>0$). These are $T_g=100^{\circ}C$ and $T_d=50^{\circ}C$ in question 7 of this project. These values are applied directly to the boundary nodes and are not updated by the finite difference equation.

proj/runHeatProject.m

@@ -0,0 +1,19 @@
+material.rho = 2700; material.Cp = 900; material.lambda = 237; material.L = 0.4;
+thermal_cond.Tg = 100; thermal_cond.Td = 50; thermal_cond.T0 = 0;
+steps.dx = 0.01; steps.dt = 0.5;
+
+t_final_s_sim = 2000; % Total simulation time
+plot_points = [0, 10, 50, 100, 200, 500, 1000, 2000]; % Specific times to capture profiles
+
+% Call the solver function
+[x_coords, final_T_profile, midpoint_T_series, time_vec, all_temp_profiles, profile_labels] = ...
+    solveHeatEquation(material, thermal_cond, steps, t_final_s_sim, plot_points);
+
+% Create and populate the figure with subplots
+figure('Color', 'white', 'Position', [100, 100, 1000, 400]);
+
+subplot(1, 2, 1);
+plotProfiles(x_coords, all_temp_profiles, profile_labels);
+
+subplot(1, 2, 2);
+plotMidpointTemperature(time_vec, midpoint_T_series);

proj/solveHeatEquation.m

@@ -0,0 +1,106 @@
+% solveHeatEquation.m
+function [x, T_final_profile, T_time_series, time_vector, all_profiles, all_profile_times] = ...
+    solveHeatEquation(material, thermal_cond, steps, t_final_s, plot_time_points)
+% SOLVEHEATEQUATION Solves the 1D Heat Equation using Explicit Finite Differences.
+%
+% Inputs:
+%   material         - Struct with fields: rho, Cp, lambda, L
+%   thermal_cond     - Struct with fields: Tg, Td, T0
+%   steps            - Struct with fields: dx, dt
+%   t_final_s        - Final simulation time in seconds (u32 or f32)
+%   plot_time_points - Array of time points (e.g., [0 10 50 100 500]) to save full profiles.
+%
+% Outputs:
+%   x                - Spatial grid points (Nx x 1 vector)
+%   T_final_profile  - Temperature profile at t_final_s (1 x Nx vector)
+%   T_time_series    - Temperature at midpoint over time (1 x Nt vector)
+%   time_vector      - Time points corresponding to T_time_series (1 x Nt vector)
+%   all_profiles     - Cell array of temperature profiles at plot_time_points
+%   all_profile_times- Corresponding time labels for all_profiles
+
+    % Extract parameters for readability
+    rho = material.rho; Cp = material.Cp; lambda = material.lambda; L = material.L;
+    Tg = thermal_cond.Tg; Td = thermal_cond.Td; T0 = thermal_cond.T0;
+    dx = steps.dx; dt = steps.dt;
+
+    % Calculate alpha and r (Fourier Number)
+    alpha = lambda / (rho * Cp);
+    r = alpha * dt / (dx^2);
+
+    if r > 0.5
+        warning('The explicit method is unstable for r > 0.5. Consider reducing dt or increasing dx.');
+    end
+
+    % Spatial Discretization
+    x = 0:dx:L;
+    Nx = length(x);
+    mid_index = round(Nx / 2); % Index for midpoint
+
+    % Initial Conditions
+    T_current = zeros(1, Nx);
+    T_current(2:Nx-1) = T0;
+    T_current(1) = Tg;
+    T_current(Nx) = Td;
+
+    % Time Discretization
+    time_vector = 0:dt:t_final_s;
+    Nt = length(time_vector);
+
+    T_time_series = zeros(1, Nt);
+    T_time_series(1) = T_current(mid_index);
+
+        % Store profiles for plotting (initialization, before the main loop)
+    all_profiles = cell(size(plot_time_points)); % Pre-allocate
+    all_profile_times = cell(size(plot_time_points));
+    
+    % Track if a profile for a given plot_time_point has been captured
+    % This ensures each target time is captured only once and correctly.
+    captured_flags = false(size(plot_time_points));
+    
+    % Capture initial profile (t=0) - this should always be the first plot point
+    % Assuming plot_time_points(1) is 0
+    all_profiles{1} = T_current;
+    all_profile_times{1} = sprintf('t = %.0f s', plot_time_points(1));
+    captured_flags(1) = true; % Mark t=0 as captured
+
+    % Main Time Loop
+    for n = 1:Nt-1
+        T_next = T_current; % Initialize T_next. Boundary conditions are retained.
+
+        % Apply finite difference equation for internal nodes
+        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 for next iteration
+
+        % Store midpoint temperature
+        T_time_series(n+1) = T_current(mid_index);
+
+    % Store full profiles at specified plot times
+    % Iterate through plot_time_points to find matches
+    for p_idx = 1:length(plot_time_points)
+        if ~captured_flags(p_idx) % Only try to capture if not already captured
+            % Check if current simulation time (time_vector(n+1)) is at or just past the target plot time
+            % Use a small tolerance relative to dt to account for floating point
+            if time_vector(n+1) >= plot_time_points(p_idx) - dt/2 && time_vector(n+1) < plot_time_points(p_idx) + dt/2
+                all_profiles{p_idx} = T_current;
+                all_profile_times{p_idx} = sprintf('t = %.0f s', plot_time_points(p_idx)); % Use the *target time* for label
+                captured_flags(p_idx) = true; % Mark as captured
+            end
+        end
+    end
+
+    % Add a check at the very end of the simulation to ensure the final plot_time_point is captured
+    % if t_final_s matches one of the plot_time_points and it wasn't captured exactly during the loop
+    [is_final_plot_time, idx_final_plot] = ismembertol(t_final_s, plot_time_points, dt/2);
+    if is_final_plot_time && ~captured_flags(idx_final_plot)
+        all_profiles{idx_final_plot} = T_current;
+        all_profile_times{idx_final_plot} = sprintf('t = %.0f s', t_final_s);
+        captured_flags(idx_final_plot) = true;
+    end
+
+    end
+    
+    T_final_profile = T_current; % The profile at the very end of the simulation
+end