1 简介
Matlab实现两端固支梁热力耦合的有限元分析
2 部分代码
% -------------------------------------------------------------------------
% 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
3 仿真结果
4 参考文献
[1]黄雅洁. 圆筒热力耦合冲击下热应力场的有限元分析及MATLAB编程[D]. 昆明理工大学, 2015.
博主简介:擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真,相关matlab代码问题可私信交流。
部分理论引用网络文献,若有侵权联系博主删除。