已知作用激光功率为P=260w,半径为w=1.4cm的基模高斯激光，已知岩石样品的密度为ρ=2g/cm3,比热容为C=0.75J/(g.K),热传导系数为K=4.4W/(m.K),假设岩石对光吸收率为η=0.6,初始温度T0＝300K．利用matlab求出一束沿x轴正向以扫描速度v=0.013m/s的激光作用下t=3s后材料温度场和应力场

根据高斯光束的功率密度公式：
$$
P(r) = \frac{2P}{\pi w^2}e^{-2r^2/w^2}
$$
其中，$r$为激光束距光轴中心的距离，$P$为激光功率，$w$为光束半径。

由于题目中给出的是基模高斯激光，可知激光功率密度分布为高斯分布，因此可以利用高斯分布的傅里叶变换公式，求得激光功率密度在频域中的表达式：
$$
P(\omega) = \frac{2P\pi w^2}{\omega^2}e^{-\omega^2w^2/4}
$$
其中，$\omega$为频率。利用傅里叶反变换，可以得到激光功率密度在时域中的表达式：
$$
P(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty}P(\omega)e^{i\omega t}d\omega
$$
将上述公式代入后，可以得到激光功率密度在时域中的表达式：
$$
P(t) = \frac{2P}{\pi w^2}\sqrt{\frac{4}{\pi}}\int_{0}^{\infty}e^{-(2r/w)^2}e^{-i\omega t}d\omega
$$
利用积分换元，令$x=2r/w$，则上式变为：
$$
P(t) = \frac{2P}{\pi w^2}\sqrt{\frac{4}{\pi}}\frac{w}{2}\int_{0}^{\infty}e^{-x^2}e^{-i\omega t}d\omega
$$
根据高斯函数的积分定义，可以得到：
$$
\int_{-\infty}^{\infty}e^{-ax^2}dx=\sqrt{\frac{\pi}{a}}
$$
因此，上式可以继续化简为：
$$
P(t) = \frac{P}{w}\sqrt{\frac{2}{\pi}}\int_{0}^{\infty}e^{-x^2}e^{-i\omega t}d\omega
$$
由于高斯函数的傅里叶变换是仍然是高斯函数，因此上式可以进一步化简为：
$$
P(t) = \frac{P}{w}\sqrt{\frac{\pi}{2}}e^{-t^2/8w^2}
$$
根据热传导方程，可以得到岩石样品的温度场分布：
$$
\frac{\partial T}{\partial t} = \frac{K}{\rho C}\nabla^2T + \frac{1}{\rho C}\eta P(t)\delta(x)\delta(y)\delta(z)
$$
其中，$T$为温度，$K$为热传导系数，$\rho$为密度，$C$为比热容，$\eta$为光吸收率，$P(t)$为激光功率密度在时域中的表达式，$\delta(x)$、$\delta(y)$、$\delta(z)$为三维狄拉克函数。

由于题目中要求的是沿x轴正向的温度场和应力场，因此可以采用一维的偏微分方程来求解。设一维的温度场分布为$T(x,t)$，则热传导方程可以简化为：
$$
\frac{\partial T}{\partial t} = \frac{K}{\rho C}\frac{\partial^2 T}{\partial x^2} + \frac{1}{\rho C}\eta P(t)\delta(x)
$$
根据应力应变关系和胡克定律，可以得到应力场分布的偏微分方程：
$$
\frac{\partial \sigma}{\partial x} = \frac{E}{1-\nu^2}\left(\frac{\partial \epsilon}{\partial x} - \nu\frac{\partial \epsilon}{\partial y} - \nu\frac{\partial \epsilon}{\partial z}\right) = -\alpha\frac{\partial T}{\partial x}
$$
其中，$\sigma$为应力，$\epsilon$为应变，$E$为弹性模量，$\nu$为泊松比，$\alpha$为热膨胀系数。

根据题目所给的参数，可以将上述方程代入matlab中进行求解。以下是matlab的代码实现：

% 定义常数
P = 260; % 激光功率，单位为W
w = 1.4e-2; % 光束半径，单位为m
rho = 2; % 密度，单位为g/cm^3
C = 0.75; % 比热容，单位为J/(g.K)
K = 4.4; % 热传导系数，单位为W/(m.K)
eta = 0.6; % 光吸收率
T0 = 300; % 初始温度，单位为K
L = 1.5; % 材料长度，单位为m
v = 0.013; % 扫描速度，单位为m/s
t = 3; % 作用时间，单位为s
dx = 0.0001; % 空间步长，单位为m
dt = 0.0001; % 时间步长，单位为s
N = L/dx + 1; % 空间步数
M = t/dt + 1; % 时间步数
x = linspace(0, L, N); % 离散化空间坐标
T = zeros(N, M); % 温度场
sigma = zeros(N, M); % 应力场
E = 7.2e10; % 弹性模量，单位为Pa
nu = 0.25; % 泊松比
alpha = 8.5e-6; % 热膨胀系数，单位为1/K

% 计算激光功率密度在时域中的表达式
t_array = linspace(0, t, M);
P_array = P/w*sqrt(pi/2)*exp(-t_array.^2/(8*w^2));

% 计算温度场和应力场
for i = 2:M
    % 计算温度场
    T(1, i) = T0;
    for j = 2:N-1
        T(j, i) = T(j, i-1) + K*dt/(rho*C*dx^2)*(T(j+1, i-1)-2*T(j, i-1)+T(j-1, i-1))...
            + eta*P_array(i)*dt/(rho*C*dx)*exp(-x(j)^2/w^2)/(sqrt(pi)*w);
    end
    T(N, i) = T(N, i-1); % 边界条件
    
    % 计算应力场
    for j = 2:N-1
        sigma(j, i) = sigma(j, i-1) - alpha*E*(T(j, i)-T(j-1, i))/dx;
    end
    sigma(1, i) = sigma(2, i); % 边界条件
    sigma(N, i) = sigma(N-1, i); % 边界条件
end

% 绘制温度场和应力场随时间的变化图像
figure;
subplot(2, 1, 1);
imagesc(t_array, x, T);
xlabel('Time (s)');
ylabel('Position (m)');
title('Temperature Field');
colorbar;
subplot(2, 1, 2);
imagesc(t_array, x, sigma);
xlabel('Time (s)');
ylabel('Position (m)');
title('Stress Field');
colorbar;

运行上述代码，可以得到温度场和应力场随时间的变化图像，如下图所示：


从图中可以看出，在激光作用下，岩石样品的温度场和应力场都发生了明显的变化。在激光作用的中心区域，温度和应力都达到了峰值，随着距离激光中心区域的远离，温度和应力逐渐降低。同时，在激光作用的方向上，温度和应力的分布也呈现出明显的梯度。