基于matlab有限体积法的传热结构拓扑优化程序

设计域

设计结果

代码

主代码（备注：程序没有经过优化，所以看起来比较多，请指正）

clc;
clear;
Lx = 0.2; % length x-direction
Ly = 0.2; % length y -direction
cv_x=100;cv_y=100;
ni = cv_x+1; % grid points x-direction
nj = cv_y+1; % grid points y-direction
ei = cv_x+2; ej = cv_y+2;   %虚拟单元
Dx = Ly/cv_x; Dy = Ly/cv_y;     
dx = Dy; dy = Dy;
xcord = 0:Dx:Lx;
ycord = 0:Dy:Ly;
DV =Dx*Dy;
Ex = 0-dx/2:dx:(ei-1)*dx;
Ey = 0-dx/2:dy:(ej-1)*dy;
bT0=[1 5];bT1 =[200,200];
%% 
lameda0 = 900;    %900 W/m^2*K
cp = 690;    %690 J/kg*K
rho = 2220;  %(kg/m^3)
Tinit =293.15;
qB1 = 5E6;
maxtime = 100;
p=1.5;
T(1:ni,1:nj) = Tinit;
theta(1:ei,1:ej) = 0;
theta(2:ei-1,2:ej-1) =0.8;
volfrac=0.50;
r=0.005;
for loop = 1:300
    thetaold = theta;
    Elameda=lameda0*theta.^p;
    T = thermalFVMout(T,Elameda,ni,nj,Dx,Dy,Tinit,qB1,rho,cp,maxtime,bT0,bT1);
    c = sum(sum(T))/(ni*nj);
    dc = calcdc(ni,nj,ei,ej,Dy,Dx,lameda0,T,theta,p);
    dc= SensitivityFilter(Ex,Ey,ei,ej,r,dc,theta);
    theta(2:ei-1,2:ej-1) = MMA(ei-2,ej-2,theta(2:ei-1,2:ej-1),volfrac,-dc(2:ei-1,2:ej-1));
    change = max(max(abs(theta(2:ei-1,2:ej-1)-thetaold(2:ei-1,2:ej-1))));
    disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
        ' Vol.: ' sprintf('%6.3f',sum(sum(theta(2:ni-1,2:nj-1)))/((ni-2)*(nj-2))) ...
        ' ch.: ' sprintf('%6.3f',change)])
    cStore(loop) =c;
    % PLOT DENSITIES
%     thetaAll(:,:,loop)=theta;
%     TAll(:,:,loop)=T;
    figure(1);
    pcolor(T(1:ni,1:nj));xlabel('x');ylabel('y');title('Temperature distribution');colorbar;
    filename = ['C:\Users\DELL018\Desktop\Thermal30C01\T\T_',num2str(loop),];
    saveas(gcf,filename,'bmp');
    figure(2);
    pcolor(theta(2:ei-1,2:ej-1));shading interp;xlabel('x');ylabel('y');title('Theta distribution');colorbar;
    filename = ['C:\Users\DELL018\Desktop\Thermal30C01\Theta\T_',num2str(loop),];
    saveas(gcf,filename,'bmp');
    figure(3); colormap(flipud(gray)); imagesc(theta(2:ei-1,2:ej-1)); axis equal; axis tight;colorbar;
    filename = ['C:\Users\DELL018\Desktop\Thermal30C01\Theta02\T_',num2str(loop),];
    saveas(gcf,filename,'bmp');
    figure(4);
    plot(cStore);
end
%     figure(1);
%     pcolor(T(2:ni-1,2:nj-1));xlabel('x');ylabel('y');title('Theta distribution');colorbar;
%     figure(2);
%     pcolor(theta(3:ei-2,3:ej-2));shading interp;xlabel('x');ylabel('y');title('Theta distribution');colorbar;
 save('dataAll','T','theta','cStore');

有限体积求解

function T = thermalFVMout(T,Elameda,ni,nj,Dx,Dy,Tinit,qB1,rho,cp,maxtime,bT0,bT1)
DV =Dx*Dy;
dx = Dx;
dy = Dy;
dt = 0.01;%传热时间
volfrac = 0.5;
lamedaEW(1:ni+1,1:nj)=0;
lamedaNS(1:ni,1:nj+1)=0;
for i=1:ni+1
    for j = 1:nj
        if j == 1
            lamedaEW(i,j) = Elameda(i,j+1);
        elseif j==nj
            lamedaEW(i,j) = Elameda(i,j);
        else
            lamedaEW(i,j) = (Elameda(i,j)+Elameda(i,j+1))*0.5;
        end
    end
end

for j=1:nj+1
    for i = 1:ni
        if i == 1
            lamedaNS(i,j) = Elameda(i+1,j);
        elseif i==ni
            lamedaNS(i,j) = Elameda(i,j);
        else
            lamedaNS(i,j) = (Elameda(i,j)+Elameda(i+1,j))*0.5;
        end
    end
end
for i=2:ni-1
    for j=2:nj-1
        aE(i,j) = Dy/(dx/lamedaEW(i+1,j));
        aW(i,j) = Dy/(dx/lamedaEW(i,j));
        aN(i,j) = Dx/(dy/lamedaNS(i,j+1));
        aS(i,j) = Dx/(dy/lamedaNS(i,j));
        aP0 = rho*cp*Dx*Dy/dt;
        SP(i,j) = 0;
        SC(i,j) = 0;
        b(i,j) = 0;
        SC(i,j) = qB1;
        aP(i,j) = aE(i,j)+aW(i,j)+aN(i,j)+aS(i,j)+aP0-SP(i,j)*Dx*Dy;
    end
end

T(1:ni,1:nj) = Tinit;
maxres = 1.0e-6;
maxit = 100;
time=0;
omega=1;
while (time < (maxtime+dt/2))
    Told=T;
    sumres = 1;
    counter = 0;
    while (sumres>maxres&counter<maxit)
        sumres = 0;
        for i = 2:ni-1
            for j = 2:nj-1
                T(i,j)=omega*(aE(i,j)*T(i+1,j)+aW(i,j)*T(i-1,j)+aN(i,j)*T(i,j+1)...
                    +aS(i,j)*T(i,j-1)+SC(i,j)*DV+aP0*Told(i,j))/aP(i,j)+(1-omega)*T(i,j);
                res =  abs(aP(i,j)*T(i,j)-(aE(i,j)*T(i+1,j)+aW(i,j)*T(i-1,j)+...
                    aN(i,j)*T(i,j+1)+aS(i,j)*T(i,j-1)+SC(i,j)*DV+aP0*Told(i,j)));
                sumres=sumres+res;
            end
        end
        %% 边界处理
        i = 1;
        for j=2:nj-1
            if (j>=bT0(1))&&(j<=bT0(2))
                T(i,j) = 293.15;
            else
                T(i,j) = 1/4*(2*T(i+1,j)+T(i,j-1)+T(i,j+1)+qB1*dx*dy/lamedaEW(i+1,j)/2);
            end
        end
        i = ni;
        for j=2:nj-1
            if(j>bT1(1))&&(j<=bT1(2))
                T(i,j) = 293.15;
            else
            T(i,j) = 1/4*(2*T(i-1,j)+T(i,j-1)+T(i,j+1)+qB1*dx*dy/lamedaEW(i,j)/2);
            end
        end
        
        j = 1;
        for i=2:ni-1
                if (i>=bT0(1))&&(i<=bT0(2))
                T(i,j) = 293.15;
                else
                T(i,j) = 1/4*(2*T(i,j+1)+T(i-1,j)+T(i+1,j)+qB1*dx*dy/lamedaNS(i,j+1)/2);
                end
        end
        j = nj;
        for i=2:ni-1
            if (i>=bT1(1))&&(i<=bT1(2))
                T(i,j) = 293.15;
            else
            T(i,j) = 1/4*(2*T(i,j-1)+T(i-1,j)+T(i+1,j)+qB1*dx*dy/lamedaNS(i,j)/2);
            end
        end
        %% 四个角点处理
        T(1,1) = 1/2*(T(2,1)+T(1,2)+qB1*dy*dx/lamedaEW(i+1,j)/4);
        T(ni,1) = 1/2*(T(ni-1,1)+T(ni,2)+qB1*dy*dx/lamedaEW(i,j)/4);
        T(1,nj) = 1/2*(T(2,nj)+T(1,nj-1)+qB1*dy*dx/lamedaEW(i+1,j)/4);
        T(ni,nj) = 1/2*(T(ni-1,nj)+T(ni,nj-1)+qB1*dy*dx/lamedaEW(i,j)/4);
        counter = counter + 1;
    end
    time = time + dt;
end

敏度求解

function dc = calcdc(ni,nj,ei,ej,Dy,Dx,lameda0,T,theta,p)

Elameda=lameda0*theta.^p;
Ap(1:nj,1:nj) = 0;
lamedaEW(1:ni+1,1:nj)=0;
lamedaNS(1:ni,1:nj+1)=0;

for i=1:ni+1
    for j = 1:nj
        if j == 1
            lamedaEW(i,j) = Elameda(i,j+1);
        elseif j==nj
            lamedaEW(i,j) = Elameda(i,j);
        else
            lamedaEW(i,j) = (Elameda(i,j)+Elameda(i,j+1))*0.5;
        end
    end
end

for j=1:nj+1
    for i = 1:ni
        if i == 1
            lamedaNS(i,j) = Elameda(i+1,j);
        elseif i==ni
            lamedaNS(i,j) = Elameda(i,j);
        else
            lamedaNS(i,j) = (Elameda(i,j)+Elameda(i+1,j))*0.5;
        end
    end
end
%+++++++++++++++++++++++++++++++++++
for i=1:ni+1
    for j=1:nj
        kh(i,j) = lamedaEW(i,j);
        Ah(i,j)=kh(i,j)*Dy/Dx;
    end
end
for j=1:nj+1
    for i=1:ni
        kv(i,j) = lamedaNS(i,j);
        Av(i,j)=kv(i,j)*Dx/Dy;
    end
end


for i=1:ni
    for j=1:nj
        Ap(i,j) = Ah(i,j)+Ah(i+1,j)+Av(i,j)+Av(i,j+1);
    end
end
dk = p*(theta.^(p-1))*lameda0;
dc = zeros(ei,ej);
%% 目标函数倒数求解(边界是否要单独处理)
for i=2:ei-1
    for j=2:ej-1
        KK = [Ap(i-1,j-1)   -Ah(i,j-1)        0               -Av(i-1,j)
                -Ah(i,j-1)       Ap(i,j-1)        -Av(i,j)            0
                  0                 -Av(i,j)         Ap(i,j)            -Ah(i,j)
                -Av(i-1,j)              0           -Ah(i,j)         Ap(i-1,j)];
        dKK = dk(i,j)*[1  -1/2   0   -1/2
                     -1/2   1   -1/2   0
                       0   -1/2   1   -1/2
                     -1/2   0   -1/2  1];
        TT =[T(i-1,j-1),T(i,j-1),T(i,j),T(i-1,j)]';
        Lameda = -(1/4*ones(1,4)/KK)';
        %dc(i,j) = Lameda'*dKK*TT;
        dc(i,j) = TT'*dKK*TT;
    end
end

MMA程序更新设计变量

function [xnew]=MMA1(nelx,nely,x,volfrac,dc)
xlow=0.001; xhigh=1;
L=x-0.1*(xhigh-xlow)*ones(nelx,nely);
high=(x-L).^2.*-dc./(xlow-L).^2;
low=(x-L).^2.*-dc./(xhigh-L).^2; 
l2 = min(min(high));
l1 = max(max(low));
for i=1:50
lmid = 0.5*(l2+l1);
xnew=max(xlow,min(xhigh,L+abs(x-L).*sqrt(-dc./lmid)));
if sum(sum(xnew)) - volfrac*nelx*nely > 0;
l2 = lmid;
else
l1 = lmid;
end
end

过滤函数

function dfnew=SensitivityFilter(Ex,Ey,ei,ej,r,df,theta)
dfnew(1:ei,1:ej)=0;
for i=2:ei-1
    for j=2:ej-1
         sum1 = 0;
         sum2 = 0;
         temp1=0;
         temp2=0;
        for i2=2:ei-1
            for j2=2:ej-1
                if (((Ex(i2)-Ex(i))^2+(Ey(j2)-Ey(j))^2)<r^2)
                    temp1=theta(i2,j2)*df(i2,j2)*(r-sqrt((Ex(i2)-Ex(i))^2+(Ey(j2)-Ey(j))^2));
                    sum1 = sum1 + temp1;
                     temp2=theta(i,j)*(r-sqrt((Ex(i2)-Ex(i))^2+(Ey(j2)-Ey(j))^2));
                    sum2 = sum2 + temp2;
                else
                    ;
                end
            end
        end
        dfnew(i,j) = sum1/sum2;
    end
end     

参考文献

《Topology optimization using the finite volume method》
《Topology Optimization Using the SIMP Method for Multiobjective Conductive Problems》