已知激光功率激光功率为P=600w,半径为w=1cm的基模高斯激光；已知岩石样品的密度为ρ=2g/cm3,比热容为C=0.75J/(g.K),热传导系数为K=4.4W/(m.K),假设岩石对光吸收率为η=0.6，岩石为长 x=10cm，宽y=10cm，高z=15cm的长方体。初始条件：温度T0＝300K，激光从（0，0，0）点开始照射，沿x轴具有v=0.13m/s的速度。岩石下表面和侧面设置为热绝缘边界，上表面为激光照射面，激光热流密度作为上表面的边界条件，利用matlab根据有限差分法计算岩石在被激光照射3s后岩石上表面沿着x轴的温度场和应力场

这是一个比较复杂的热传导问题，需要使用有限差分法进行数值计算。以下是具体的步骤和代码实现：

1. 确定计算区域和网格大小

根据题目中给出的岩石尺寸和激光照射方向，可以确定计算区域的大小为 10cm × 10cm × 15cm，可以将其分成网格进行计算。为了保证计算精度，需要选择足够小的网格大小，例如 0.1cm × 0.1cm × 0.1cm。

2. 确定时间步长和总时间

在有限差分法中，时间步长的大小对计算结果有很大影响。根据题目中给出的激光照射时间为 3s，可以选择一个合适的时间步长，例如 0.01s。则总时间步数为 300。

3. 确定初始条件和边界条件

根据题目中给出的初始条件和边界条件，可以设置初始温度场、边界条件和材料参数。

% 初始温度场
T0 = 300; % 初始温度
T = T0 * ones(nx, ny, nz); % 初始化温度场

% 边界条件
q = @(t) 0.6 * 600 * pi * w^2 * exp(-((t-1.5)/0.5)^2); % 激光热流密度
T(:, :, 1) = q(0) / (rho * c) + T0; % 上表面边界
T(:, :, end) = T0; % 下表面边界和侧面边界

% 材料参数
rho = 2; % 岩石密度
c = 0.75; % 岩石比热容
k = 4.4; % 岩石热传导系数

4. 根据有限差分法进行计算

根据热传导方程和有限差分法的离散形式，可以编写以下代码进行计算。

for n = 1 : nt
    T_old = T;
    for i = 2 : nx-1
        for j = 2 : ny-1
            for k = 2 : nz-1
                T(i, j, k) = T_old(i, j, k) + dt * (k/(rho*c)) * ((T_old(i+1, j, k) - 2*T_old(i, j, k) + T_old(i-1, j, k))/(dx^2) + (T_old(i, j+1, k) - 2*T_old(i, j, k) + T_old(i, j-1, k))/(dy^2) + (T_old(i, j, k+1) - 2*T_old(i, j, k) + T_old(i, j, k-1))/(dz^2));
            end
        end
    end
    % 边界条件
    T(:, :, 1) = q(n*dt) / (rho * c) + T0;
    T(:, :, end) = T0;
end

5. 计算应力场

根据热应力理论，可以计算出岩石的应力场。由于岩石是各向同性材料，可以使用简化的 Hooke 定律进行计算。

alpha = k / (rho * c); % 热扩散系数
E = 73e9; % 岩石杨氏模量
nu = 0.25; % 岩石泊松比
alpha_T = 1 / (3*k) * (3*alpha - 4*E*nu/(1-2*nu)); % 热膨胀系数
sigma_x = zeros(nx, ny, nz); % 初始化应力场
for i = 2 : nx-1
    for j = 2 : ny-1
        for k = 2 : nz-1
            delta_T = T(i, j, k) - T0;
            sigma_x(i, j, k) = E / (1 - nu^2) * (delta_T - alpha_T * T0);
        end
    end
end

6. 绘制温度场和应力场

可以使用 MATLAB 的图形界面工具箱或自行编写代码绘制温度场和应力场的三维图像。

参考代码如下：

% 绘制温度场
[x, y, z] = meshgrid(1:nx, 1:ny, 1:nz);
figure;
xslice = [1, round(nx/2), nx];
yslice = [1, round(ny/2), ny];
zslice = [1, round(nz/2), nz];
slice(x, y, z, T, xslice, yslice, zslice);
xlabel('x (cm)');
ylabel('y (cm)');
zlabel('z (cm)');
colorbar;
title('Temperature Field');

% 绘制应力场
figure;
slice(x, y, z, sigma_x, xslice, yslice, zslice);
xlabel('x (cm)');
ylabel('y (cm)');
zlabel('z (cm)');
colorbar;
title('Stress Field');