The Corresponding Files (Click to Save):算法

Code: top.m
References: A 99 line topology optimization code written in MATLAB

The program can be promoted by line:數組

top(30,10,0.5,3.0,1.5)

代碼註釋

之前學了點皮毛就直接用商業軟件搬磚了，總以爲有點心虛，因此入門第一步，拜讀了Sigmund的99行Matlab拓撲優化程序，參考了其餘大神們的註釋，做爲初學者還有很多看不懂的地方，我本身又補充整理了一遍。我是零基礎小白，因此代碼前面先交代了一些理論背景，我本身能看懂相信全部人都能看懂，哈哈哈。app

%%%%%%%%  A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000  %%%%%%%%
%%%%%%%%  COMMENTED - OUT 1.0 BY HAOTIAN_W,                       JULY 2020  %%%%%%%%
function top(nelx,nely,volfrac,penal,rmin)
% ===================================================================================
% nelx    : 水平方向上的離散單元數;
% nely    : 豎直方向上的離散單元數;
%
% volfrac : 容積率，材料體積與設計域體積之比，對應的工程問題就是"將結構減重到百分之多少";
%
% penal   : 懲罰因子，SIMP方法是在0-1離散模型中引入連續變量x、係數p及中間密度單元，從而將離
%           散型優化問題轉換成連續型優化問題，而且令0≤x≤1，p爲懲罰因子，經過設定p>1對中間密
%           度單元進行有限度的懲罰，儘可能減小中間密度單元數目，使單元密度儘量趨於0或1;
% 
%           合理選擇懲罰因子的取值，能夠消除多孔材料，從而獲得理想的拓撲優化結果:
%               當penal<=2時     存在大量多孔材料，計算結果沒有可製造性;
%               當penal>=3.5時   最終拓撲結果沒有大的改變;
%               當penal>=4時     結構整體柔度的變化很是緩慢，迭代步數增長，計算時間延長;
%            
% rmin    : 敏度過濾半徑，防止出現棋盤格現象;
% ===================================================================================
% 結構優化的數學模型經常使用名詞:
%   1) 設計變量   在設計中可調整的、變化的基本參數，本算例是單元密度;
%   2) 目標函數   設計變量的函數，優化設計的目標，本算例是柔度最小;
%   3) 約束條件   幾何、應力、位移等，難點在創建約束方程，本算例是幾何體積約束;
%   4) 終止準則   結束迭代的條件，本算例是目標函數變化量<=0.010;
%   5) 載荷工況   定義結構全部可能的受力狀況，本算例是單一載荷;
% ===================================================================================
% 基於SIMP理論的優化準則法迭代分析流程:
%   1) 定義設計域，選擇合適的設計變量、目標函數以及約束函數等其餘邊界條件;
%   2) 結構離散爲有限元網格，計算優化前的單元剛度矩陣;
%   3) 初始化單元設計變量，即給定設計域內的每一個單元一個初始單元相對密度;
%   4) 計算各離散單元的材料特性參數，計算單元剛度矩陣，組裝結構總剛度矩陣，計算結點位移;
%   5) 計算整體結構的柔度值及其敏度值，求解拉格朗日乘子;
%   6) 用優化準則方法進行設計變量更新;
%   7) 檢查結果的收斂性，如未收斂則轉4)，循環迭代，如收斂則轉8);
%   8) 輸出目標函數值及設計變量值，結束計算。
% ===================================================================================

% x是設計變量，給設計域內的每一個單元一個初始相對密度，值爲volfrac
x(1:nely,1:nelx) = volfrac; 

% loop儲存迭代次數
loop = 0; 

% change儲存每次迭代以後目標函數的改變值，用以判斷是否收斂
change = 1.;

% 當目標函數改變量<=0.01時說明收斂，結束迭代
while change > 0.01
    loop = loop + 1; 
    
    % 將前一次的設計變量賦值給xold，x用來儲存這一次的結果，以後還要比較它們以判斷是否收斂
    xold = x;  

    % 每次迭代都進行一次有限元分析，計算結點位移，並儲存在全局位移數組U中
    [U] = FE(nelx,nely,x,penal);   
    
    % 計算單元剛度矩陣
    [KE] = lk;
    
    % c是用來儲存目標函數的變量
    c = 0.;
    
    % 遍歷設計域矩陣元素，從左上角第一個元素開始，一列一列
    for ely = 1:nely
        for elx = 1:nelx
            
            % 節點位移儲存在U中，若是想得到節點位移必須先知道節點編號，進行索引;
            % 節點編號能夠根據當前單元在設計域中的位置算出;
            % n1是左上角節點編號，n2是右上角節點編號;
            n1 = (nely+1)*(elx-1)+ely; 
            n2 = (nely+1)*elx+ely;
           
            % 局部位移數組Ue儲存4個節點共8個自由度位移，每一個節點分別有x、y兩個自由度; 
            % 由於是矩形單元，因此根據n一、n2兩個節點的編號能夠推演出單元全部節點的自由度編號;
            % 順序是：[左上x；左上y；右上x；右上y；右下x；右下y；左下x；左下y];
            % 只適用於矩形單元劃分網格;
            Ue = U([2*n1-1;2*n1; 2*n2-1;2*n2; 2*n2+1;2*n2+2; 2*n1+1;2*n1+2],1);
         
            % SIMP模型，將設計變量x從離散型變成指函數型，指數就是懲罰因子：x(ely,elx)^penal
            % 計算整體結構的柔度值，這裏目標函數是柔度最小，參見論文中公式（1）
            c = c + x(ely,elx)^penal*Ue'*KE*Ue;
            
            % 計算整體結構的敏度值，實際上dc就是c對Xe的梯度，參見論文中公式（4）
            dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;
        end
    end
    
    % 無關網格敏度過濾
    [dc] = check(nelx,nely,rmin,x,dc); 
    
    % 採用優化準則法（OC）求解當前模型，得出知足體積約束的結果，更新設計變量
    [x] = OC(nelx,nely,x,volfrac,dc);
    
    % 更新目標函數改變值
    change = max(max(abs(x-xold)));
    
    % 打印迭代信息: It.迭代次數，Obj.目標函數，Vol.材料體積比，ch.迭代改變量
    disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
       ' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...
        ' ch.: ' sprintf('%6.3f',change )])
    
    % 當前迭代結果圖形顯示
    colormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);
end


% 採用優化準則法（OC）迭代子程序 -------------------------------------------------------
% 數學模型主要求解算法有：優化準則法(OC)、序列線性規劃法(SLP)、移動漸進線法(MMA); 
% OC適用於單約束優化問題求解，好比這裏的"體積約束下的柔度最小化"問題，當求解複雜的多約束拓撲
% 優化問題時，採用SLP和MMA一般更方便;
% 參見論文中公式（2）（3）
function [xnew] = OC(nelx,nely,x,volfrac,dc)
% Input：  水平單元數nelx, 豎直單元數nely, 設計變量x, 材料體積比volfrac, 目標函數靈敏度dc;
% Output： 更新後的設計變量xnew;

% 定義一個取值區間，二分法，獲得知足體積約束的拉格朗日算子
l1 = 0; l2 = 100000;

% 正向最大位移
move = 0.2;

while (l2-l1 > 1e-4)
    
    % 二分法，取區間中點
    lmid = 0.5*(l2+l1);
  
    % 參見論文中的公式（2）
    % sqrt(-dc./lmid)對應公式中Be^eta（eta=1/2），eta阻尼係數是爲了確保計算的收斂性
    xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));
    
    % sum(sum(xnew))是更新後的材料體積, volfrac*nelx*nely是優化目標，用它們的差值判斷是否收斂
    % sum()能夠去看看help，注意一下行、列問題，不看也行，反正知道sum(sum())是全部元素求和就行
    if sum(sum(xnew)) - volfrac*nelx*nely > 0
        l1 = lmid;
    else
        l2 = lmid;
    end
end


% 無關網格敏度過濾子程序 --------------------------------------------------------------
% 參見論文中公式（5）（6）
function [dcn] = check(nelx,nely,rmin,x,dc)
% Input：  水平單元數nelx, 豎直單元數nely, 敏度過濾半徑rmin, 設計變量x, 整體結構敏度dc;
% Output： 過濾後的目標函數敏度dcn;

% dcn清零，用來保存更新的目標函數靈敏度
dcn = zeros(nely,nelx);

for i = 1:nelx   
    for j = 1:nely
        sum=0.0; 
        
        % 在過濾半徑定義的範圍內遍歷
        for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx)
            for l = max(j-floor(rmin),1):min(j+floor(rmin),nely)
                
                % 參見論文中公式（6）,fac即公式中卷積算子Hf
                % qrt((i-k)^2+(j-l)^2) 是計算此單元與相鄰單元的距離，即公式中dist(e,f)
                fac = rmin-sqrt((i-k)^2+(j-l)^2);
                sum = sum+max(0,fac);
                
                % 參見論文中公式（5）
                dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);
            end
        end
        dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
    end
end


% 有限元求解子程序 --------------------------------------------------------------------
function [U] = FE(nelx,nely,x,penal)
% Input：  水平單元數nelx, 豎直單元數nely, 設計變量x, 懲罰因子penal;
% Output： 全局節點位移U;

% 計算單元剛度矩陣
[KE] = lk; 

% 整體剛度矩陣的稀疏矩陣
K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));

% 力矩陣的稀疏矩陣
F = sparse(2*(nely+1)*(nelx+1),1); 

% U清零，用來保存更新的全局節點位移
U = zeros(2*(nely+1)*(nelx+1),1);

for elx = 1:nelx
    for ely = 1:nely
        
        % 計算單元左上角、右上角節點編號
        n1 = (nely+1)*(elx-1)+ely; 
        n2 = (nely+1)* elx   +ely;
        
        % 同上主程序，計算單元4個節點8個自由度
        edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];
        
        % 將單元剛度矩陣KE 組裝成 整體剛度矩陣K
        K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;
    end
end

% 施加載荷，本算例應用了一個在左上角的垂直單元力
F(2,1) = -1;

% 施加約束，消除線性方程中的固定自由度來實現支承結構，本算例左邊第一列和右下角固定
fixeddofs = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);

% 剩下的不加約束的節點自由度，setdiff()從..中除去..
alldofs   = [1:2*(nely+1)*(nelx+1)];
freedofs  = setdiff(alldofs,fixeddofs);

% 求解線性方程組，獲得各節點自由度的位移值儲存在U中
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:); 

% 受約束節點固定自由度位移值爲0
U(fixeddofs,:)= 0;


% 求解單元剛度矩陣子程序 --------------------------------------------------------------
% 有限元方法計算的一個重要的係數矩陣，表徵單元體的受力與變形關係;
% 特色：對稱性、奇異性、主對角元素恆正、奇數（偶數）行和爲0;
% 矩形單元4節點 8*8矩陣;
function [KE] = lk

% 材料楊氏彈性模量
E = 1.; 

% 材料泊松比
nu = 0.3;

k=[ 1/2-nu/6   1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ... 
   -1/4+nu/12 -1/8-nu/8  nu/6       1/8-3*nu/8];
KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
                  k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
                  k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
                  k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
                  k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
                  k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
                  k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
                  k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];

公式彙總

柔度最小目標函數	$\left.\begin{array}{rl} \min _{\mathbf{x}}: & c(\mathbf{x})=\mathbf{U}^{T} \mathbf{K} \mathbf{U}=\sum_{e=1}^{N}\left(x_{e}\right)^{p} \mathbf{u}_{e}^{T} \mathbf{k}_{0} \mathbf{u}_{e} \\ \text { subject to }: & \frac{V(\mathbf{x})}{V_{0}}=f \\ : & \mathbf{K U}=\mathbf{F} \\ : & \mathbf{0}<\mathbf{x}_{\min } \leq \mathbf{x} \leq \mathbf{1} \end{array}\right\}$	（1）
OC準則優化更新	$\begin{aligned} &x_{e}^{\text {new }}=\\ &\left\{\begin{array}{c} \max \left(x_{\min }, x_{e}-m\right) \text { if } x_{e} B_{e}^{\eta} \leq \max \left(x_{\min }, x_{e}-m\right) \\ x_{e} B_{e}^{\eta} \text { if } \max \left(x_{\min }, x_{e}-m\right)<x_{e} B_{e}^{\eta}<\min \left(1, x_{e}+m\right) \\ \min \left(1, x_{e}+m\right) \text { if } \min \left(1, x_{e}+m\right) \leq x_{e} B_{e}^{\eta} \end{array}\right. \end{aligned}$	（2）
	$B_{e}=\frac{-\frac{\partial c}{\partial x_{e}}}{\lambda \frac{\partial V}{\partial x_{e}}}$	（3）
目標函數靈敏度	$\frac{\partial c}{\partial x_{e}}=-p\left(x_{e}\right)^{p-1} \mathbf{u}_{e}^{T} \mathbf{k}_{0} \mathbf{u}_{e}$	（4）
靈敏度濾波	$\frac{\widehat{\partial c}}{\partial x_{e}}=\frac{1}{x_{e} \sum_{f=1}^{N} \hat{H}_{f}} \sum_{f=1}^{N} \hat{H}_{f} x_{f} \frac{\partial c}{\partial x_{f}}$	（5）
卷積算子(加權因子)	$\begin{array}{l} \hat{H}_{f}=r_{\min }-\operatorname{dist}(e, f) \\ \left\{f \in N \mid \operatorname{dist}(e, f) \leq r_{\min }\right\}, \quad e=1, \ldots, N \end{array}$	（6）

參考資料

[1] Sigmund, O. A 99 line topology optimization code written in Matlab. Struct Multidisc Optim 21, 120–127 (2001). ide

[TopOpt] 99行拓撲優化程序徹底註釋

代碼註釋

公式彙總

參考資料