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

Question

Emmanuel Sarpong 2024-5-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2123451-i-have-written-a-code-to-reproduce-fig-2-in-the-attached-paper-but-i-don-t-seem-to-get-it-right

回答： Arnav 2024-9-5

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 个评论
显示无隐藏无

Dyuman Joshi 2024-5-28

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 2024-5-28

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 2024-9-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2123451-i-have-written-a-code-to-reproduce-fig-2-in-the-attached-paper-but-i-don-t-seem-to-get-it-right#answer_1511229

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

2 个评论
显示无隐藏无

回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

2 个评论 显示 无隐藏 无

回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论