利用W-M公式的三维分形接触模型的matlab代码，D等于2.2，G等于10

以下是利用W-M公式的三维分形接触模型的MATLAB代码，其中D等于2.2，G等于10：

clear all;
clc;
n=100; %生成的分形体积大小
D=2.2; %分形维数
G=10; %分形体积的粗糙度
L=1; %分形体积大小
z=zeros(n,n,n); %初始化分形体积
z(1,:,:)=0.5;
z(:,1,:)=0.5;
z(:,:,1)=0.5;
for k=1:n
    for j=1:n
        for i=1:n
            if i==1&&j==1&&k==1
                continue;
            elseif i==1&&j==1
                z(i,j,k)=z(i,j,k-1)^(D-1)+z(i,j+1,k-1)^(D-1)+z(i+1,j,k-1)^(D-1);
                z(i,j,k)=z(i,j,k)/3;
            elseif i==1&&k==1
                z(i,j,k)=z(i,j-1,k)^(D-1)+z(i,j+1,k)^(D-1)+z(i+1,j,k)^(D-1);
                z(i,j,k)=z(i,j,k)/3;
            elseif j==1&&k==1
                z(i,j,k)=z(i,j,k-1)^(D-1)+z(i+1,j,k)^(D-1)+z(i,j+1,k)^(D-1);
                z(i,j,k)=z(i,j,k)/3;
            elseif i==1
                z(i,j,k)=z(i,j,k-1)^(D-1)+z(i,j-1,k)^(D-1)+z(i,j+1,k)^(D-1)+z(i+1,j,k)^(D-1);
                z(i,j,k)=z(i,j,k)/4;
            elseif j==1
                z(i,j,k)=z(i,j,k-1)^(D-1)+z(i,j+1,k)^(D-1)+z(i-1,j,k)^(D-1)+z(i+1,j,k)^(D-1);
                z(i,j,k)=z(i,j,k)/4;
            elseif k==1
                z(i,j,k)=z(i,j-1,k)^(D-1)+z(i,j+1,k)^(D-1)+z(i-1,j,k)^(D-1)+z(i+1,j,k)^(D-1);
                z(i,j,k)=z(i,j,k)/4;
            else
                z(i,j,k)=z(i,j-1,k)^(D-1)+z(i,j+1,k)^(D-1)+z(i-1,j,k)^(D-1)+z(i+1,j,k)^(D-1)+z(i,j,k-1)^(D-1)+z(i,j,k-1)^(D-1);
                z(i,j,k)=z(i,j,k)/6;
            end
        end
    end
end
[x,y,z]=meshgrid(1:n,1:n,1:n);
isosurface(x,y,z,z,G);
axis equal;

这段代码与前面的代码基本相同，只是增加了一个参数G来控制分形体积的粗糙度。具体来说，G越大，等值面就越粗糙，反之，G越小，等值面就越光滑。在这里，我们将G设置为10，可以看到生成的分形体积的等值面比之前更加粗糙了。