I have written a code to reproduce fig 2 in the attached paper but I don't seem to get it right.
7 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I am trying a code to reproduce fig 2 in the attached paper, but I don't seem to get it right. Any help would be apreciated. Thanks
% Material properties
rho = 1.35e3; % Density in kg/m^3
Cp = @(T) (-0.07827 + 0.00223*T - 8.66702e-7*T.^2 + 9.63112e-11*T.^3) * 1e3; % Specific heat in J/g°C converted to J/kgK
alpha = 3.5e3; % Absorption coefficient in m^-1
R = 0.1; % Reflectivity
pulseWidth = 6e-9; % Pulse width in seconds
fluence = 790e-3; % Laser fluence in J/cm^2
I0 = fluence / pulseWidth; % Peak laser intensity
% Spatial domain
L = 1e-3; % Length of the sample in meters (1 mm)
num_elements = 100; % Number of spatial elements
x = linspace(0, L, num_elements); % Spatial grid
% Time domain
total_time = 50e-9; % Total time in seconds
num_time_steps = 500; % Number of time steps
t = linspace(0, total_time, num_time_steps); % Time vector
% Initial temperature distribution
T0 = 20; % Initial temperature in °C
T = T0 * ones(num_elements, 1); % Initial temperature vector
% Thermal diffusivity (D = k / (rho * Cp))
D = 0.4e-4; % Thermal diffusivity in m^2/s
k = D * rho * Cp(T0); % Thermal conductivity in W/m°C
% FEM setup
dx = L / (num_elements - 1); % Spatial step size
dt = total_time / num_time_steps; % Time step size
% FEM matrices
A = zeros(num_elements, num_elements);
B = zeros(num_elements, 1);
% Assembly of FEM matrices
for i = 2:num_elements-1
A(i, i-1) = k / dx^2;
A(i, i) = -2 * k / dx^2;
A(i, i+1) = k / dx^2;
end
% Boundary conditions
A(1, 1) = 1;
A(1, 2) = 0;
A(end, end) = 1;
A(end, end-1) = 0;
% Time-stepping solution
T_history = zeros(num_elements, num_time_steps);
T_history(:, 1) = T;
for n = 1:num_time_steps-1
B(2:end-1) = alpha * I0 * exp(-alpha * x(2:end-1)) * (1 - R);
B(end) = T0; % Bottom boundary condition (fixed temperature)
T = T + dt * (A * T + B);
T(1) = T0; % surface boundary condition
T_history(:, n+1) = T;
end
% Plotting the results
figure;
hold on;
plot(t, T_history(1, :), 'r', 'DisplayName', 'Surface');
plot(t, T_history(find(x >= 1e-6, 1), :), 'g', 'DisplayName', '1 \mum');
plot(t, T_history(find(x >= 5e-6, 1), :), 'b', 'DisplayName', '5 \mum');
xlabel('Time (s)');
ylabel('Temperature (°C)');
legend;
title('Temperature profiles at different depths');
hold off;
2 commentaires
Dyuman Joshi
le 28 Mai 2024
Please provide information and context as to what the figure is and what formulae/methods are being used to obtain the figure.
The paper does not have any required details about your question or your code, it seems.
Réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!