I have written a code to reproduce fig 2 in the attached paper but I don't seem to get it right.

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;

[Dyuman Joshi]

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.

[Emmanuel Sarpong]

The paper uses the 1D differential heat conduction equation that's equation 1 and the parameter, I(X,t) is given by equation 2. They apply the finite element method (FEM) to do the computations. The quantities, density (P) in g/cm^3, heat capacity(Cp) in Jg^-1K^-1, diffusivity (D) in cm^2/s, absorption coefficient (alpha) at 532 nm given in cm^-1. The thermal conductivity is the product of diffusivity, density and the specific heat.

Answer 1

[Arnav]

Hi @Emmanuel Sarpong,

I understand that you are trying to solve the partial differential equation

where

Instead of implementing a custom algorithm for solving the partial differential equation using finite element method, it is advisable to use the pdepe function in MATLAB to perform the same for you.

To use pdepe you first need to set up the equation.

function [c, f, s] = pde_function(x, t, T, dTdx) 
    C_p = 0.07827 + 0.00223 * T - 8.66702e-7 * T^2 + 9.63112e-11 * T^3; 
    alpha = k / (q * C_p); 
    c = 1; 
    f = alpha * dTdx; 
    I_t = I0_peak * exp(-((t - 0)^2 / (pulse_duration / (2 * sqrt(2 * log(2))))^2)); 
    s = (1 - R_532) * I_t * exp(-a_532 * x) / (q * C_p); 
end 

After setting up the equation and necessary boundary and initial conditions. You can use the pdepe function as follows:

sol = pdepe(0, @pde_function, initial_condition, boundary_conditions, x, t); 

You can learn more about modelling different PDEs and setting up boundary and initial conditions from the documentation page of pdepe: https://www.mathworks.com/help/matlab/ref/pdepe.html

Regards,

Arnav