拓扑优化经典99行程序解读

拓扑优化的99行程序学习报告4月19日2011《结构优化设计》课程学习报告任课教师：李书一、前言：在最近的结构优化设计课程上学习了O.Sigmund的《A 99 line topology optimization code written in Matlab》一文，对拓扑优化的理论原理与实际的计算机程序实现都有了一定的理解，文章主要是通过拓扑优化的原理来实现对简单结构的静力学问题的优化求解，而编写的代码仅有99行，包括36行的主程序，12行的OC优化准则代码，16行的网格过滤代码和35行的有限元分析代码。

自1988 年丹麦学者Bendsoe与美国学者Kikuchi提出基于均匀化方法的结构拓扑优化设计基本理论以来,均匀化方法应用到具有周期性结构的材料分析中,近几年该方法已经成为分析夹杂、纤维增强复合材料、混凝土材料等效模量,以及材料的细观结构拓扑优化常用的手段之一。

其基本思想是在组成拓扑结构的材料中引入微结构，优化过程中以微结构的几何尺寸作为设计变量，以微结构的消长实现其增删，并产生介于由中间尺寸微结构组成的复合材料，从而实现了结构拓扑优化模型与尺寸优化模型的统一。

文章就是通过均匀化的基础，结合拓扑结构优化的工程实际，以计算机模拟的方法将拓扑优化的一般过程呈现出来，有助于初涉拓扑优化的读者对拓扑优化有个基础的认识。

二、拓扑优化问题描述为了简化问题的描述，文中假设设计域是简单的矩形形式，且在进行有限元离散的时候采用正方形单元对其进行离散。

这样不仅便于进行单元离散和单元编号，也利于对结构进行几何外形的描述。

一般说来，基于指数逼近法的拓扑优化最小化的问题可作如下描述：文中采用的对结构材料属性的描述是所谓的“指数逼近法”或者称为SIMP 逼近法，即（Solid Isotropic Material with Penalization带惩罚因子的各项同性材料模型法），该方法是拓扑优化中常用的变密度材料插值模型中最具代表性的一种。

拓扑优化方法拓扑优化方法是一种有效的优化方法，目前被广泛应用于求解复杂优化问题。

本文通过介绍拓扑优化方法的基本原理、典型案例、优势与应用等方面，来深入探讨拓扑优化的相关知识。

一、什么是拓扑优化方法拓扑优化方法（Topology Optimization，简称TO）是一种解决复杂最优化问题的有效优化方法，它是利用拓扑的可变性，用于求解复杂拓扑结构组合优化问题的一种新兴方法。

拓扑优化方法既可以用来求解有限元分析（Finite Element Analysis，简称FEA）中有序结构问题，也可以用来求解无序结构问题。

二、拓扑优化方法的基本原理拓扑优化方法的基本原理是：在设定的最优化目标函数及运算范围内，利用优化技术，使得复杂结构拓扑结构达到最优，从而达到最优化设计目标。

拓扑优化方法的优势主要体现在重量最小化、强度最大化、结构疲劳极限优化等多种反向设计问题上。

此外，由于拓扑优化方法考虑到结构加工、安装、维护等方面，其结构设计更加实用性好。

三、拓扑优化方法的典型案例1、航空外壳优化：目前，航空外壳的拓扑优化设计可以使得外壳的重量减轻50%以上，同时提升外壳的强度和耐久性。

2、机械联轴器优化：拓扑优化方法可以有效的提高机械联轴器长期使用的耐久性，减少其体积和重量，满足高性能要求。

3、结构优化：通过拓扑优化方法，可以有效地减少刚性框架结构的重量，优化结构设计，改善结构性能，大大降低制造成本。

四、拓扑优化方法的优势1、灵活性强：拓扑优化方法允许在设计过程中改变结构形态，可以有效利用具有局部不稳定性的装配元件；2、更容易操作：拓扑优化方法比传统的有序结构模型更容易实现，不需要做过多的运算；3、成本低：拓扑优化方法可以有效降低产品的工艺制造成本，在改进出色性能的同时，可以节省大量人力物力；4、可重复性高：拓扑优化方法可以实现由抽象到具体的可重复的设计，可以实现大量的应用系统。

五、拓扑优化方法的应用拓扑优化方法目前被广泛应用在机械、航空航天、汽车等机械工程领域，具体应用包括但不限于：机械手和夹具的设计优化，汽车机架优化，电器结构优化，机械外壳优化，振动优化，和结构强度优化等等。

拓扑优化的99行程序学习报告4月19日2011《结构优化设计》课程学习报告任课教师：李书一、 前言：在最近的结构优化设计课程上学习了O.Sigmund 的《A 99 line topology optimization code written in Matlab 》一文，对拓扑优化的理论原理与实际的计算机程序实现都有了一定的理解，文章主要是通过拓扑优化的原理来实现对简单结构的静力学问题的优化求解，而编写的代码仅有99行，包括36行的主程序，12行的OC 优化准则代码，16行的网格过滤代码和35行的有限元分析代码。

自1988 年丹麦学者Bendsoe 与美国学者Kikuchi 提出基于均匀化方法的结构拓扑优化设计基本理论以来,均匀化方法应用到具有周期性结构的材料分析中,近几年该方法已经成为分析夹杂、纤维增强复合材料、混凝土材料等效模量,以及材料的细观结构拓扑优化常用的手段之一。

其基本思想是在组成拓扑结构的材料中引入微结构，优化过程中以微结构的几何尺寸作为设计变量，以微结构的消长实现其增删，并产生介于由中间尺寸微结构组成的复合材料，从而实现了结构拓扑优化模型与尺寸优化模型的统一。

文章就是通过均匀化的基础，结合拓扑结构优化的工程实际，以计算机模拟的方法将拓扑优化的一般过程呈现出来，有助于初涉拓扑优化的读者对拓扑优化有个基础的认识。

二、 拓扑优化问题描述为了简化问题的描述，文中假设设计域是简单的矩形形式，且在进行有限元离散的时候采用正方形单元对其进行离散。

这样不仅便于进行单元离散和单元编号，也利于对结构进行几何外形的描述。

一般说来，基于指数逼近法的拓扑优化最小化的问题可作如下描述：01min :()()():0min 1NT p Te e e Xe c X U KU x u k u V X subjecttof V KU FX x =⎫==⎪⎪⎪=⎬⎪⎪=⎪<≤≤⎭∑文中采用的对结构材料属性的描述是所谓的“指数逼近法”或者称为SIMP逼近法，即（Solid Isotropic Material with Penalization 带惩罚因子的各项同性材料模型法），该方法是拓扑优化中常用的变密度材料插值模型中最具代表性的一种。

1.%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY2000 %%%2.%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLESIGMUND %%%3.%%%% 一个由 OLE SIGMUND编写的99行拓扑优化代码，2000年1月 %%%4.%%%% 为加速而修改的代码，2002年9月,由OLE SIGMUND编写 %%%5.function top(nelx,nely,volfrac,penal,rmin);6.% INITIALIZE7.x(1:nely,1:nelx) = volfrac;8.loop = 0;9.change = 1.;10.% START ITERATION11.while change > 0.0112.loop = loop + 1;13.xold = x;14.% FE-ANALYSIS15.[U]=FE(nelx,nely,x,penal);16.% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS17.[KE] = lk;18.c = 0.;19.for ely = 1:nely20.for elx = 1:nelx21.n1 = (nely+1)*(elx-1)+ely;22.n2 = (nely+1)* elx +ely;23.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);24.c = c + x(ely,elx)^penal*Ue'*KE*Ue;25.dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;26.end27.end28.% FILTERING OF SENSITIVITIES29.[dc] = check(nelx,nely,rmin,x,dc);30.% DESIGN UPDATE BY THE OPTIMALITY CRITERIA METHOD31.[x] = OC(nelx,nely,x,volfrac,dc);32.% PRINT RESULTS33.change = max(max(abs(x-xold)));34.disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...35.' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...36.' ch.: ' sprintf('%6.3f',change )])37.% PLOT DENSITIES38.colormap(gray); imagesc(-x); axis equal; axis tight; axisoff;pause(1e-6);39.end40.%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%41.function [xnew]=OC(nelx,nely,x,volfrac,dc)42.l1 = 0; l2 = 100000; move = 0.2;43.while (l2-l1 > 1e-4)44.lmid = 0.5*(l2+l1);45.xnew =max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));46.if sum(sum(xnew)) - volfrac*nelx*nely > 0;47.l1 = lmid;48.else49.l2 = lmid;50.end51.end52.%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%53.function [dcn]=check(nelx,nely,rmin,x,dc)54.dcn=zeros(nely,nelx);55.for i = 1:nelx56.for j = 1:nely57.sum=0.0;58.for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx)59.for l = max(j-floor(rmin),1):min(j+floor(rmin),nely)60.fac = rmin-sqrt((i-k)^2+(j-l)^2);61.sum = sum+max(0,fac);62.dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);63.end64.end65.dcn(j,i) = dcn(j,i)/(x(j,i)*sum);66.end67.end68.%%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%69.function [U]=FE(nelx,nely,x,penal)70.[KE] = lk;71.K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));72.F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);73.for elx = 1:nelx74.for ely = 1:nely75.n1 = (nely+1)*(elx-1)+ely;76.n2 = (nely+1)* elx +ely;77.edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];78.K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;79.end80.end81.% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)82.F(2,1) = -1;83.fixeddofs = union([1:2:2*(nely+1)],[2*(nelx+1)*(nely+1)]);84.alldofs = [1:2*(nely+1)*(nelx+1)];85.freedofs = setdiff(alldofs,fixeddofs);86.% SOLVING87.U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);88.U(fixeddofs,:)= 0;89.%%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%90.function [KE]=lk91.E = 1.;92.nu = 0.3;93.k=[ 1/2-nu/6 1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ...94.-1/4+nu/12 -1/8-nu/8 nu/6 1/8-3*nu/8];95.KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)96.k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)97.k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)98.k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)99.k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)100.k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)101.k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)102.k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];103.% 104.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%105.% This Matlab code was written by Ole Sigmund, Department of Solid % 106.% Mechanics, Technical University of Denmark, DK-2800 Lyngby, Denmark. %107.% Please sent your comments to the author:[email]sigmund@fam.dtu.dk[/email]%108.% %109.% The code is intended for educational purposes and theoretical details %110.% are discussed in the paper %111.% "A 99 line topology optimization code written in Matlab" % 112.% by Ole Sigmund (2001), Structural and Multidisciplinary Optimization, %113.% Vol 21, pp. 120--127. %114.% %115.% The code as well as a postscript version of the paper can be % 116.% downloaded from the web-site:[url]http://www.topopt.dtu.dk[/url] %117.% %118.% Disclaimer: %119.% The author reserves all rights but does not guaranty that the code is %120.% free from errors. Furthermore, he shall not be liable in any event %121.% caused by the use of the program.。

3188-1-1.htmlSigmund教授所编写的top优化经典99行程序，可以说是我们拓扑优化研究的基础；每一个新手入门都会要读懂这个程序，才能去扩展，去创新；99行程序也有好多个版本，用于求解各种问题，如刚度设计、柔顺机构、热耦合问题，但基本思路大同小异；本文拟对其中的一个版本进行解读，愿能对新手有点小小的帮助。

不详之处，还请论坛内高手多指点读懂了该程序，只能说是略懂拓扑优化理论了，我手里就有一些水平集源程序是成千上万行，虽然在99行的基础上成熟了很多，但依然还有很多的发展空间。

源程序如下：%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%function top(nelx,nely,volfrac,penal,rmin);nelx=80;nely=20;volfrac=0.4;penal=3;rmin=2;% INITIALIZEx(1:nely,1:nelx) = volfrac;loop = 0;change = 1.;% START ITERATIONwhile change > 0.01loop = loop + 1;xold = x;% FE-ANAL YSIS[U]=FE(nelx,nely,x,penal);% OBJECTIVE FUNCTION AND SENSITIVITY ANAL YSIS[KE] = lk;c = 0.;for ely = 1:nelyfor elx = 1:nelxn1 = (nely+1)*(elx-1)+ely;n2 = (nely+1)* elx +ely;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);c = c + x(ely,elx)^penal*Ue'*KE*Ue;dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;endend% FILTERING OF SENSITIVITIES[dc] = check(nelx,nely,rmin,x,dc);% DESIGN UPDA TE BY THE OPTIMALITY CRITERIA METHOD[x] = OC(nelx,nely,x,volfrac,dc);% PRINT RESULTSchange = max(max(abs(x-xold)));disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...' ch.: ' sprintf('%6.3f',change )])% PLOT DENSITIEScolormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);end%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [xnew]=OC(nelx,nely,x,volfrac,dc)l1 = 0; l2 = 100000; move = 0.2;while (l2-l1 > 1e-4)lmid = 0.5*(l2+l1);xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));if sum(sum(xnew)) - volfrac*nelx*nely > 0;l1 = lmid;elsel2 = lmid;endend%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dcn]=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);for i = 1:nelxfor j = 1:nelysum=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)fac = rmin-sqrt((i-k)^2+(j-l)^2);sum = sum+max(0,fac);dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);endenddcn(j,i) = dcn(j,i)/(x(j,i)*sum);endend%%%%%%%%%%FE-ANAL YSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function [U]=FE(nelx,nely,x,penal)[KE] = lk;K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);for elx = 1:nelxfor ely = 1:nelyn1 = (nely+1)*(elx-1)+ely;n2 = (nely+1)* elx +ely;edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;endend% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)F(2*(nelx/2+1)*(nely+1),1) = 1;fixeddofs = [2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1)];alldofs = [1:2*(nely+1)*(nelx+1)];freedofs = setdiff(alldofs,fixeddofs);% SOLVINGU(freedofs, ：）= K(freedofs,freedofs) \ F(freedofs,：）;U(fixeddofs,：）= 0; % 这两行复制后应换成英文字符，我这里为了防止生成QQ表% 情修改了一下格式%%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [KE]=lkE = 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)];% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%程序执行方法：（just for matlab new users）打开matlab，点开new M-file，将上述源程序复制粘贴到M-文件中，修改蓝色部分的格式，保存。

%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%%%%% CODE MODIFIED FOR INCREASED SPEED, September 2002, BY OLE SIGMUND %%%function top(nelx,nely,volfrac,penal,rmin);nelx=80;nely=20;volfrac=0.4;penal=3;rmin=2;% INITIALIZEx(1:nely,1:nelx) = volfrac;loop = 0;change = 1.;% START ITERA TIONwhile change > 0.01loop = loop + 1;xold = x;% FE-ANAL YSIS[U]=FE(nelx,nely,x,penal);% OBJECTIVE FUNCTION AND SENSITIVITY ANAL YSIS[KE] = lk;c = 0.;for ely = 1:nelyfor elx = 1:nelxn1 = (nely+1)*(elx-1)+ely;n2 = (nely+1)* elx +ely;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);c = c + x(ely,elx)^penal*Ue'*KE*Ue;dc(ely,elx) = -penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue;endend% FILTERING OF SENSITIVITIES[dc] = check(nelx,nely,rmin,x,dc);% DESIGN UPDA TE BY THE OPTIMALITY CRITERIA METHOD[x] = OC(nelx,nely,x,volfrac,dc);% PRINT RESULTSchange = max(max(abs(x-xold)));disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...' ch.: ' sprintf('%6.3f',change )])% PLOT DENSITIEScolormap(gray); imagesc(-x); axis equal; axis tight; axis off;pause(1e-6);end%%%%%%%%%% OPTIMALITY CRITERIAUPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [xnew]=OC(nelx,nely,x,volfrac,dc)l1 = 0; l2 = 100000; move = 0.2;while (l2-l1 > 1e-4)lmid = 0.5*(l2+l1);xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));if sum(sum(xnew)) - volfrac*nelx*nely > 0;l1 = lmid;elsel2 = lmid;endend%%%%%%%%%% MESH-INDEPENDENCYFILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dcn]=check(nelx,nely,rmin,x,dc)dcn=zeros(nely,nelx);for i = 1:nelxfor j = 1:nelysum=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)fac = rmin-sqrt((i-k)^2+(j-l)^2);sum = sum+max(0,fac);dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);endenddcn(j,i) = dcn(j,i)/(x(j,i)*sum);endend%%%%%%%%%%FE-ANAL YSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%function [U]=FE(nelx,nely,x,penal)[KE] = lk;K = sparse(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);for elx = 1:nelxfor ely = 1:nelyn1 = (nely+1)*(elx-1)+ely;n2 = (nely+1)* elx +ely;edof = [2*n1-1; 2*n1; 2*n2-1; 2*n2; 2*n2+1; 2*n2+2; 2*n1+1; 2*n1+2];K(edof,edof) = K(edof,edof) + x(ely,elx)^penal*KE;endend% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)F(2*(nelx/2+1)*(nely+1),1) = 1;fixeddofs = [2*(nely/2+1),2*nelx*(nely+1)+2*(nely/2+1)];alldofs = [1:2*(nely+1)*(nelx+1)];freedofs = setdiff(alldofs,fixeddofs);% SOLVINGU(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);U(fixeddofs,:)= 0;%%%%%%%%%% ELEMENT STIFFNESSMATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [KE]=lkE = 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)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%。