【物理应用】Matlab实现两端固支梁热力耦合的有限元分析

% ------------------------------------------------------------------------- % BASIC PARAMETERS (BRIDGE-LIKE STRUCTURE) % ------------------------------------------------------------------------- clear ; clc ; HoriTotal = 100 ; % horizontal VeriTotal = 70; % vertical IniVolfrac = 0.5 ; % initial volume fraction Volfrac = 0.25 ; % allowable volume fraction varimin = 1e-3 ; % smallest value of design variable ObjScale = 14e-3 ; rmin = 2 ; % filter radius F_uni = 35 ; % distributed load qnVari = 0.3 ; % STM applied to design variables RE = 28 ; % RAMP for stiffness Rbt = 16 ; % RAMP for thermal stress coefficient Talafa0 = 12.1e-6 ; % thermal expansion coefficient TwoDeVec = [ 1 1 0 ]; % ------------------------------------------------------------------------- % DISCRETIZATION FEATURES % ------------------------------------------------------------------------- %-ELEMENT STIFFNESS- a = 0.5 ; b = 0.5; h = 1 ; % element thickness v0 = (2*a) * (2*b) * h ; % volume of solid element E0 = 2.1e5 ; % elastic modulus Emin = 0 * E0 ; NU = 0.3 ; si = -1 ; ti = 1 ; sj = 1 ; tj = 1 ; sm = 1 ; tm = -1 ; sp = -1 ; tp = -1 ; ID = 1 ; % ID = 1: plane stress problem, ID = 2: plane strain problem dce = zeros( Syselem ,1 ) ; for loopi = 1 : Syselem c = c + 0.5 * (xPhys(loopi)/(1 + RE * (1-xPhys(loopi)))) * nodes_d(loopi,:) * E0 * KE * nodes_d(loopi,:)’ ; % structural compliance dce(loopi) = nodes_d(loopi,:) * (((1 + Rbt)/(1+Rbt*(1-xPhys(loopi)))^2 * E0 * Talafa0 ) * Feth0(:,loopi))... - 0.5 * (1+RE)/(1 + RE * (1-xPhys(loopi)))^2 * nodes_d(loopi,:) * E0 * KE * nodes_d(loopi,:)’ ; % sensitivity of structural compliance end Compli = c ; dce(:) = H * (dce(:)./Hs); dve = ones(Syselem,1).* v0/(Syselem * v0 * Volfrac); % volume constraint dve(:) = H * (dve(:)./Hs); %-METHOD OF MOVING ASYMPTOTES- m = 1; % number of constraint functions a1 = zeros(m,1); % Column vector with the constants a_i in the terms a_i*z. c_MMA = 10000*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i. d = zeros(m,1); % Columns vector with the constants d_i in the terms 0.5*d_i*(y_i)^2. xval = xDes; f0val = c; % compliance minimization df0dx = dce(:); fval = sum(xPhys.*v0)/(Syselem * v0 * Volfrac)-1 ; % volume constraint dfdx = dve’ ; [xmma, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ... mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ... f0val,df0dx,fval,dfdx,low,upp,a0,a1,c_MMA,d); % Update MMA Variables xnew = reshape(xmma , Syselem , 1 ); xold2 = xold1(:); xold1 = xDes(:); xDes = xnew ; xDes = xold1 + qnVari * (xDes-xold1) ; changeVari = max(abs(xDes-xold1)) ; %-UPDATE MODEL DATA- xPhys = (H * xDes(:))./Hs; % physical density volfrac = 100 * mean(xPhys(:)) ; % volume fraction v = sum(xPhys) * v0 * ObjScale ; %-PRINT RESULTS AND PLOT DENSITIES- disp([’ It.: ’ sprintf(’%-4i’,loop) ’ c.: ’ sprintf(’%-7.2f’,Compli)... ’ Volfrac.: ’ sprintf(’%-7.2f’,volfrac) ’Volume.: ’ sprintf(’%-7.2f’,v ) ... ’ chVari.: ’ sprintf(’%-7.3f’,changeVari ) ’MaxTdif.: ’ sprintf(’%-7.1f’,max(TEdiffe) ) ’Penal.: ’ sprintf(’%-7.0f’, RE )]) xPhysFig = reshape(xPhys ,nely,nelx) ; colormap(gray); imagesc(1-xPhysFig); axis equal; axis tight; axis off;pause(1e-6); end