采用MATLAB，在主应力空间内，绘制mises屈服面，采用方式为：1）沿着静水压力轴线每次递进一步，围绕垂直静水压轴线的面，旋转一周，寻找屈服面上在该静水压下的坐标，，之后再低进一步，重复上一部操作，直到找出屈服面上的所有点，最后将屈服面绘制，并将屈服面划分成四边形

绘制mises屈服面的步骤如下：

1. 确定静水压力轴线，并按照静水压力的递进值，依次计算每个静水压力下的主应力大小和方向。

2. 对于每个静水压力下的主应力方向，以其为旋转轴，将垂直于该方向的面旋转一周，得到该静水压力下的所有应力状态。

3. 对于每个应力状态，计算它的von Mises应力，找出屈服面上的点。

4. 重复步骤2和步骤3，直到找出所有静水压力下的屈服面上的点。

5. 将所有的点连接起来，形成屈服面，并将其划分成四边形。

6. 使用MATLAB中的绘图函数，如scatter和patch，将屈服面绘制出来。

以下是MATLAB代码示例：

% 输入参数
P = [0:10:100]; % 静水压力递进值
sigma1 = 200; % 最大主应力
sigma2 = 150; % 中间主应力
sigma3 = 100; % 最小主应力

% 计算主应力方向
[sigma_x,sigma_y,sigma_z] = sph2cart(deg2rad(45),deg2rad(30),1); % 以x轴为起点，旋转45度，俯角为30度
sigma1_dir = [sigma_x,sigma_y,sigma_z];
sigma2_dir = cross([0,0,1],sigma1_dir);
sigma2_dir = sigma2_dir/norm(sigma2_dir);
sigma3_dir = cross(sigma1_dir,sigma2_dir);

% 计算Mises应力
von_Mises = zeros(length(P),360);
for i = 1:length(P)
    p = P(i);
    for j = 1:360
        theta = deg2rad(j);
        sigma = p*[0,0,1] + sigma1*sigma1_dir*cos(theta)^2 + sigma2*sigma2_dir*sin(theta)*cos(theta) + sigma3*sigma3_dir*sin(theta)^2;
        von_Mises(i,j) = sqrt(sigma(1)^2+sigma(2)^2+sigma(3)^2-sigma(1)*sigma(2)-sigma(1)*sigma(3)-sigma(2)*sigma(3));
    end
end

% 绘制屈服面
figure;
hold on;
for i = 1:length(P)
    p = P(i);
    for j = 1:360
        if von_Mises(i,j) <= von_Mises(1,1)
            scatter3(von_Mises(i,j)*sigma1_dir(1),von_Mises(i,j)*sigma1_dir(2),von_Mises(i,j)*sigma1_dir(3),'k','filled'); % 在该点处绘制黑色点
        end
    end
end

% 划分四边形
for i = 1:length(P)-1
    for j = 1:359
        if von_Mises(i,j) <= von_Mises(1,1) && von_Mises(i+1,j) <= von_Mises(1,1) && von_Mises(i,j+1) <= von_Mises(1,1) && von_Mises(i+1,j+1) <= von_Mises(1,1)
            patch([von_Mises(i,j)*sigma1_dir(1),von_Mises(i+1,j)*sigma1_dir(1),von_Mises(i+1,j+1)*sigma1_dir(1),von_Mises(i,j+1)*sigma1_dir(1)],...
                  [von_Mises(i,j)*sigma1_dir(2),von_Mises(i+1,j)*sigma1_dir(2),von_Mises(i+1,j+1)*sigma1_dir(2),von_Mises(i,j+1)*sigma1_dir(2)],...
                  [von_Mises(i,j)*sigma1_dir(3),von_Mises(i+1,j)*sigma1_dir(3),von_Mises(i+1,j+1)*sigma1_dir(3),von_Mises(i,j+1)*sigma1_dir(3)],'r'); % 绘制红色四边形
        end
    end
end

axis equal;
xlabel('sigma1');
ylabel('sigma2');
zlabel('sigma3');
view(3);