拓扑优化149行均匀化matlab代码注释详解

function CH = homo3D(lx,ly,lz,lambda,mu,voxel)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% lx = Unit cell length in x-direction.

% ly = Unit cell length in y-direction.

% lz = Unit cell length in z-direction.

% lambda = Lame's first parameter for solid materials.

% mu = Lame's second parameter for solid materials.

% voxel = Material indicator matrix. Used to determine nelx/nely/nelz

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% INITIALIZE 初始化

% voxel是输入参数，表示体素模型。这一行代码获取了体素模型的尺寸，即在x、y和z轴上的维度，

% 分别存储在'nelx','nely','nelz'中

[nelx, nely, nelz] = size(voxel); %size of voxel model along x,y and z axis

% Stiffness matrix 计算每个体素的尺寸'dx''dy''dz'，并计算总的体素数量'nel'

dx = lx/nelx; dy = ly/nely; dz = lz/nelz;

nel = nelx*nely*nelz;

% 调用'hexahedron'函数，传递体素尺寸的一半作为参数，以计算单元刚度矩阵和单元载荷向量。

% 这些值将在后续有限元分析中使用。

[keLambda, keMu, feLambda, feMu] = hexahedron(dx/2,dy/2,dz/2);

% Node numbers and element degrees of freedom for full (not periodic) mesh

% (非周期性的)网格节点编号和单元自由度。'nodenrs'是节点编号矩阵，

% 'edofVec'是包含所有节点的全局自由度编号的列向量

nodenrs = reshape(1:(1+nelx)*(1+nely)*(1+nelz),1+nelx,1+nely,1+nelz);

edofVec = reshape(3*nodenrs(1:end-1,1:end-1,1:end-1)+1,nel,1);

% 定义了一些修正系数，用于构建全局节点和自由度的映射关系。

% 'addx'和'addy'中包含了一些修正值，

% 然后使用'repmat'函数将它们重复添加到'edofVec'中，以生产全局自由度矩阵'edof'

addx = [0 1 2 3*nelx+[3 4 5 0 1 2] -3 -2 -1];

addxy = 3*(nely+1)*(nelx+1)+addx;

edof = repmat(edofVec,1,24) + repmat([addx addxy],nel,1); % 全局自由度矩阵，包含每个自由度编号

% 总的来说，这个初始化过程设置了有限元分析中所需的一些基本参数和数据结构，

% 包括体素模型的尺寸、刚度矩阵、全局节点编号和自由度。这些信息将在后续的有限元分析中使用。

%% IMPOOSE PERIODIC BOUNDARY CONDITIONS 施加周期性边界条件

% Use original edofMat to index into list with the periodic dofs

% 'nn'表示总节点数，'nnP'表示唯一节点数，即在非周期情况下的节点数

nn = (nelx+1)*(nely+1)*(nelz+1); % Total number of nodes 总结点数

nnP = (nelx)*(nely)*(nelz); % Total number of unique nodes 唯一节点数

% 将非周期节点编号整理成三维矩阵，以便进行周期性边界条件扩充

nnPArray = reshape(1:nnP, nelx, nely, nelz); % 将单元编号整理成三维矩阵

% 将节点编号在三个方向上分别进行扩充，扩充的方式是通过在边界处添加镜像节点，以模糊周期性边界条件

% Extend with a mirror of the back border 扩充镜像后边界条件

nnPArray(end+1,:,:) = nnPArray(1,:,:);

% Extend with a mirror of the left border 扩充镜像做边界条件

nnPArray(:, end+1, :) = nnPArray(:,1,:);

% Extend with a mirror of the top border 扩充镜像上边界条件

nnPArray(:, :, end+1) = nnPArray(:,:,1);

% Make a vector into which we can index using edofMat: 制作一个向量用于edofMat索引

% 创建一个向量'dofVector'，用于存储周期性边界条件扩充后的节点自由度。

% 通过将'nnPArray'中的节点编号进行索引，生成具有周期性边界条件的节点自由度

dofVector = zeros(3*nn, 1);

dofVector(1:3:end) = 3*nnPArray(:)-2;

dofVector(2:3:end) = 3*nnPArray(:)-1;

dofVector(3:3:end) = 3*nnPArray(:);

% 通过将'dofVector'更新全局自由度矩阵'edof',以考虑周期性边界条件

edof = dofVector(edof);

% 计算考虑周期性边界条件后的总自由度数量。

ndof = 3*nnP;

% 总的来说，这段代码的作用是在有限元分析中考虑周期性边界条件，

% 通过将边界处的节点进行镜像扩充，更新全局自由度矩阵，以模拟系统的周期性行为。

% 这种操作的目的是简化模型，减小计算域的尺寸，从而降低计算的复杂度。

%% ASSEMBLE GLOBAL STIFFNESS MATRIX AND LOAD VECTORS 组装全局刚度矩阵和载荷向量

% Indexing vectors

iK = kron(edof,ones(24,1))';

jK = kron(edof,ones(1,24))';

% Material properties assigned to voxels with materials

% 将材料属性 lambda 和 mu 限制在具有材料的体素上。

% 这是通过将条件 (voxel==1) 应用于 lambda 和 mu 实现的，

% 其中 (voxel==1) 为真表示该体素内有材料，为假表示无材料。

lambda = lambda*(voxel==1); mu = mu*(voxel==1);

% The corresponding stiffness matrix entries

% 计算刚度矩阵 K 中的条目，考虑了材料的Lambda和Mu。: 表示将矩阵展开成列向量，.' 表示转置。

sK = keLambda(:)*lambda(:).' + keMu(:)*mu(:).';

% 使用稀疏矩阵构造函数 sparse 构建全局刚度矩阵 K，并将其对称化

K = sparse(iK(:), jK(:), sK(:), ndof, ndof);

K = 1/2*(K+K');

% Assembly three load cases corresponding to the three strain cases

% 这段代码的目的是为了构建加载向量 F。

% edof 是一个矩阵，表示每个单元的自由度。

% repmat 函数将其重复6次，以生成与加载向量中的6种应变情况相对应的索引。

iF = repmat(edof',6,1);

% 构建一个矩阵 jF，其中包含了对应于每种应变情况的类型。

% ones(24, nel) 生成一个大小为 (24, nel) 的矩阵，每一列都是对应的应变类型。

jF = [ones(24,nel); 2*ones(24,nel); 3*ones(24,nel);...

4*ones(24,nel); 5*ones(24,nel); 6*ones(24,nel);];

% 计算加载向量的值。feLambda(:) 和 feMu(:)

% 将 feLambda 和 feMu 展开成列向量，然后通过点乘 lambda 和 mu 得到对应的值。

sF = feLambda(:)*lambda(:).'+feMu(:)*mu(:).';

% 使用稀疏矩阵构造函数 sparse 构建加载向量 F。

% iF(:)、jF(:)、sF(:) 分别表示索引、类型和值。

% ndof 是总的自由度数，6 表示加载向量有6列，分别对应于六种不同的应变情况xx,yy,zz,xy,xz,yz。

% 最终，加载向量 F 的每个元素代表了在每个节点上的每个自由度上施加的外部力或位移。这是有限元分析中用于描述外部约束和荷载的一种形式。

F = sparse(iF(:), jF(:), sF(:), ndof, 6); % ndof和6将F的大小指定为ndof*6

%% SOLUTION 求解

% solve by PCG method, remember to constrain one node

% 这两行代码的目的是从整个模型的自由度中提取出属于活跃（或有材料）的单元的自由度信息，并将其整理成一个排序的、唯一的自由度数组。

% 具体解释如下：

% edof 是整个有限元模型的自由度矩阵，其中包含了每个单元的自由度信息。这个矩阵的行数对应于单元数，列数是每个单元的自由度数量的总和。

% voxel==1 是一个逻辑数组，其元素是根据 voxel 矩阵中元素是否为 1 来确定的。在有限元模型中，voxel 矩阵通常表示材料的分布，其中 1 表示有材料的位置。

% edof(voxel==1,:) 通过逻辑索引，选择出 voxel 中为 1 的位置对应的单元的自由度信息。

% unique 函数用于去除重复的自由度。这样做的原因是，一个单元可能与其他多个单元共享一些自由度。

% sort 函数用于对提取出的自由度进行排序。

% 综合起来，activedofs 最终包含了属于有材料的单元的自由度信息，这些自由度按照顺序排列，并且没有重复。

% 这个信息将在后续的有限元分析中用于构建子系统的刚度矩阵和载荷向量。

activedofs = edof(voxel==1,:); activedofs = sort(unique(activedofs(:)));

% 初始化位移矩阵X，用于存储每个自由度上的位移

X = zeros(ndof,6);

% cholesky分解

% 通过 ichol 函数计算预条件共轭梯度方法中的预条件矩阵。预条件矩阵的作用是改善共轭梯度法的收敛性。

L = ichol(K(activedofs(4:end),activedofs(4:end)));% 预条件矩阵

for i = 1:6 % 遍历六种不同的应变情况，每一个应变情况对应一个列向量，表示位移场中的一个分量

% X(activedofs(4:end),i)表示将求解得到的位移结构保存到矩阵X中

% activedofs(4:end)是活跃节点的索引，'i'表示当前循环中要求解的应变分量

X(activedofs(4:end),i) = pcg(K(activedofs(4:end),...

activedofs(4:end)),F(activedofs(4:end),i),1e-10,300,L,L');

%pcg(K(activedofs(4:end),

% ...activedofs(4:end)), F(activedofs(4:end),i), 1e-10, 300, L, L'):

% 使用 PCG 方法求解线性系统。

% 其中，K(activedofs(4:end), activedofs(4:end)) 是刚度矩阵的子集，

% F(activedofs(4:end), i) 是载荷矩阵的子集，1e-10 是收敛容差，

% 300 是最大迭代次数，L 和 L' 分别是预条件矩阵及其共轭转置。

end

% X(activedofs(4:end),:) = K(activedofs(4:end),activedofs(4:end))...

% \F(activedofs(4:end),:); % Solving by direct method

% 这是求解的另一种方式，使用直接方法 \ 求解线性方程组。

% 通常，预条件共轭梯度方法在大型有限元模型中的收敛速度更快，因此更为常用。

%% HOMOGENIZATION 均匀化

% 这部分代码实现了对微观层级的均匀化。具体而言，它计算了宏观尺度上的弹性模量'CH'

% 以描述材料整体弹性性质。

% The displacement vectors corresponding to the unit strain cases

% X0存储每个体素的位移向量，维度为(元素数量,单元节点自由度，应变情况数量)

X0 = zeros(nel, 24, 6);

% 初始化一个单元位移数组，存储六个不同应变情况下的元素位移

% The element displacements for the six unit strains

X0_e = zeros(24, 6);

% 合成的单元刚度矩阵ke和力向量fe

%fix degrees of nodes [1 2 3 5 6 12];

ke = keMu + keLambda; % Here the exact ratio does not matter, because

fe = feMu + feLambda; % it is reflected in the load vector

% 通过求解线性系统，计算六个不同应变情况下的元素位移

X0_e([4 7:11 13:24],:) = ke([4 7:11 13:24],[4 7:11 13:24])...

\fe([4 7:11 13:24],:);

% 为X0的每个应变情况赋值

X0(:,:,1) = kron(X0_e(:,1)', ones(nel,1)); % epsilon0_11 = (1,0,0,0,0,0)

X0(:,:,2) = kron(X0_e(:,2)', ones(nel,1)); % epsilon0_22 = (0,1,0,0,0,0)

X0(:,:,3) = kron(X0_e(:,3)', ones(nel,1)); % epsilon0_33 = (0,0,1,0,0,0)

X0(:,:,4) = kron(X0_e(:,4)', ones(nel,1)); % epsilon0_12 = (0,0,0,1,0,0)

X0(:,:,5) = kron(X0_e(:,5)', ones(nel,1)); % epsilon0_23 = (0,0,0,0,1,0)

X0(:,:,6) = kron(X0_e(:,6)', ones(nel,1)); % epsilon0_13 = (0,0,0,0,0,1)

% 初始化均匀化弹性张量

CH = zeros(6);

% 单胞体积

volume = lx*ly*lz;

for i = 1:6

for j = 1:6 % 用两个嵌套for循环遍历弹性张量CH的每个元素

% 计算sum_L和sum_M分别对应于kelamada和kemu部分的总和。

% 这里使用了存储在X中的位移信息，以及弹性张量的相应部分

% X是位移矩阵

sum_L = ((X0(:,:,i) - X(edof+(i-1)*ndof))*keLambda).*...

(X0(:,:,j) - X(edof+(j-1)*ndof));

sum_M = ((X0(:,:,i) - X(edof+(i-1)*ndof))*keMu).*...

(X0(:,:,j) - X(edof+(j-1)*ndof));

% 将形状调整为和体素模型相同的性质

sum_L = reshape(sum(sum_L,2), nelx, nely, nelz);

sum_M = reshape(sum(sum_M,2), nelx, nely, nelz);

% Homogenized elasticity tensor

% 计算总的全局弹性矩阵 论文公式12

CH(i,j) = 1/volume*sum(sum(sum(lambda.*sum_L + mu.*sum_M)));

end

end

end

%% 计算单元刚度矩阵和载荷向量 Ke = 积分 BDB dV，

function [keLambda, keMu, feLambda, feMu] = hexahedron(a, b, c)

% Constitutive matrix contributions 定义弹性材料的弹性矩阵的两个分量，

% 'CMu'表示剪切模量矩阵，'CLamada'拉伸模量矩阵

CMu = diag([2 2 2 1 1 1]); CLambda = zeros(6); CLambda(1:3,1:3) = 1;

% Three Gauss points in both directions 选择三个高斯积分点，用于在每个方向上数值积分。

% 'xx','yy','zz'是积分点的位置，'ww'是相应的权重

xx = [-sqrt(3/5), 0, sqrt(3/5)]; yy = xx; zz = xx;

ww = [5/9, 8/9, 5/9];

% Initialize 初始化用于存储刚度矩阵和载荷向量的矩阵

keLambda = zeros(24,24); keMu = zeros(24,24);

feLambda = zeros(24,6); feMu = zeros(24,6);

for ii = 1:length(xx) % 使用三重循环遍历高斯积分点，计算积分点的位置

for jj = 1:length(yy)

for kk = 1:length(zz)

%integration point 插值点

x = xx(ii); y = yy(jj); z = zz(kk);

%stress strain displacement matrix 计算高斯积分点相关的应变-位移矩阵

qx = [ -((y-1)*(z-1))/8, ((y-1)*(z-1))/8, -((y+1)*(z-1))/8,...

((y+1)*(z-1))/8, ((y-1)*(z+1))/8, -((y-1)*(z+1))/8,...

((y+1)*(z+1))/8, -((y+1)*(z+1))/8];

qy = [ -((x-1)*(z-1))/8, ((x+1)*(z-1))/8, -((x+1)*(z-1))/8,...

((x-1)*(z-1))/8, ((x-1)*(z+1))/8, -((x+1)*(z+1))/8,...

((x+1)*(z+1))/8, -((x-1)*(z+1))/8];

qz = [ -((x-1)*(y-1))/8, ((x+1)*(y-1))/8, -((x+1)*(y+1))/8,...

((x-1)*(y+1))/8, ((x-1)*(y-1))/8, -((x+1)*(y-1))/8,...

((x+1)*(y+1))/8, -((x-1)*(y+1))/8];

% Jacobian 计算雅可比矩阵，并使用它计算全局坐标系下的应变-位移矩阵

J = [qx; qy; qz]*[-a a a -a -a a a -a; -b -b b b -b -b b b;...

-c -c -c -c c c c c]';

qxyz = J\[qx;qy;qz];

% 计算应变-位移矩阵B

B_e = zeros(6,3,8);

for i_B = 1:8

B_e(:,:,i_B) = [qxyz(1,i_B) 0 0;

0 qxyz(2,i_B) 0;

0 0 qxyz(3,i_B);

qxyz(2,i_B) qxyz(1,i_B) 0;

0 qxyz(3,i_B) qxyz(2,i_B);

qxyz(3,i_B) 0 qxyz(1,i_B)];

end

B = [B_e(:,:,1) B_e(:,:,2) B_e(:,:,3) B_e(:,:,4) B_e(:,:,5)...

B_e(:,:,6) B_e(:,:,7) B_e(:,:,8)];

% Weight factor at this point 计算积分点的权重

weight = det(J)*ww(ii) * ww(jj) * ww(kk);

% Element matrices 单元矩阵

% 使用权重和应变-位移矩阵B计算拉伸和剪切刚度矩阵 论文公式8的两项

keLambda = keLambda + weight * B' * CLambda * B;

keMu = keMu + weight * B' * CMu * B;

% Element loads 单元载荷

% 使用权重和应变-位移矩阵计算拉伸和剪切载荷向量 论文公式11的两项

feLambda = feLambda + weight * B' * CLambda;

feMu = feMu + weight * B' * CMu;

end

end

end

end