已知作用激光功率为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后材料温度场和应力场

根据高斯光束理论，可以得到高斯光强分布公式：

$I(r)=\frac{2P}{\pi w^2}e^{-2r^2/w^2}$

其中，P为激光功率，w为激光半径，r为距离光轴的距离。

由于题目中给出的激光功率和激光半径都是在平面内的值，因此需要将激光功率单位转换为每个激光束的功率，即：

$P_{beam}=P/\pi w^2=260/(3.14*0.014^2)\approx1.6\times10^7 W/m^2$

根据光吸收率η和激光束功率密度，可以得到每单位质量岩石吸收的热量：

$q=\eta P_{beam}/\rho=0.6\times1.6\times10^7/2=4.8\times10^6 J/kg$

再根据热传导方程和材料参数，可以得到温度场的方程：

$\frac{\partial T}{\partial t}=K\frac{\partial^2 T}{\partial x^2}+\frac{q}{\rho C}$

其中，K为热传导系数，q为吸收的热量，ρ为密度，C为比热容。

应力场的方程为：

$\frac{\partial \sigma_x}{\partial x}+\frac{\partial \sigma_y}{\partial y}+\frac{\partial \sigma_z}{\partial z}=0$

$\sigma_x=\sigma_y=-\frac{\lambda}{2\mu+\lambda}\frac{\partial T}{\partial x},\sigma_z=-\frac{\lambda}{2\mu+\lambda}\frac{\partial T}{\partial z}$

其中，λ和μ为材料的Lamé常数。

使用matlab进行数值求解，采用有限差分方法，代码如下：

clear all;
clc;

% 材料参数
rho=2000;  % 密度，kg/m^3
C=750;  % 比热容，J/(kg.K)
K=4.4;  % 热传导系数，W/(m.K)
eta=0.6;  % 光吸收率
lambda=K*0.7;  % Lamé常数
mu=K/2-lambda/2;

% 计算高斯光强分布
P=260;  % 激光功率，W
w=0.014;  % 激光半径，m
P_beam=P/(pi*w^2);  % 每个激光束功率，W/m^2
[X,Y]=meshgrid(-0.05:0.001:0.05,-0.05:0.001:0.05);
r=sqrt(X.^2+Y.^2);
I=P_beam*exp(-2*r.^2/w^2);

% 初始化温度场、应力场和时间步长
T=300*ones(size(X));  % 初始温度，K
sigma_x=zeros(size(X));
sigma_y=zeros(size(X));
sigma_z=zeros(size(X));
dt=0.01;  % 时间步长，s
t=0;  % 初始时间，s
t_end=3;  % 时间终点，s
dx=0.001;  % 空间步长，m

% 进行数值求解
while t<t_end
    % 计算热传导项
    d2Tdx2=diff(T,2)/dx^2;
    d2Tdx2(:,1)=d2Tdx2(:,2);  % 处理边界
    d2Tdx2(:,end)=d2Tdx2(:,end-1);
    d2Tdz2=diff(T,2,2)/dx^2;
    d2Tdz2(1,:)=d2Tdz2(2,:);
    d2Tdz2(end,:)=d2Tdz2(end-1,:);
    q=eta*P_beam/rho;
    dqdx=diff(q.*I,1)/dx;
    dqdx(:,1)=dqdx(:,2);
    dqdx(:,end)=dqdx(:,end-1);
    dqdz=diff(q.*I,1,2)/dx;
    dqdz(1,:)=dqdz(2,:);
    dqdz(end,:)=dqdz(end-1,:);
    dTdt=K*(d2Tdx2+d2Tdz2)+dqdx/rho/C;
    T=T+dt*dTdt;
    
    % 计算应力项
    dsigmaxdx=diff(sigma_x,1)/dx;
    dsigmaxdx(:,1)=dsigmaxdx(:,2);
    dsigmaxdx(:,end)=dsigmaxdx(:,end-1);
    dsigmadz=diff(sigma_z,1,2)/dx;
    dsigmadz(1,:)=dsigmadz(2,:);
    dsigmadz(end,:)=dsigmadz(end-1,:);
    sigma_x=sigma_x+dt*(-lambda/(2*mu+lambda)*dsigmaxdx);
    sigma_y=sigma_y+dt*(-lambda/(2*mu+lambda)*dsigmaxdx);
    sigma_z=sigma_z+dt*(-lambda/(2*mu+lambda)*dsigmadz);
    
    % 更新时间
    t=t+dt;
end

% 可视化温度场和应力场
figure;
contourf(X,Y,T);
colorbar;
title('Temperature Field');
xlabel('x (m)');
ylabel('z (m)');

figure;
contourf(X,Y,sigma_z);
colorbar;
title('Stress Field');
xlabel('x (m)');
ylabel('z (m)');

运行结果如下：