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;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; 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; finite element method MATLAB Answers — New Questions

​