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

% -------------------------------------------------------------------------%                  BASIC PARAMETERS (BRIDGE-LIKE STRUCTURE)    % -------------------------------------------------------------------------clear ; clc ;HoriTotal  = 100 ;          % horizontalVeriTotal  = 70;            % verticalIniVolfrac = 0.5  ;         % initial volume fractionVolfrac    = 0.25 ;         % allowable volume fractionvarimin    = 1e-3 ;         % smallest value of design variableObjScale   = 14e-3 ;        rmin       = 2 ;            % filter radiusF_uni      = 35 ;           % distributed loadqnVari     = 0.3 ;          % STM applied to design variablesRE         = 28 ;           % RAMP for stiffnessRbt        = 16 ;           % RAMP for thermal stress coefficientTalafa0   = 12.1e-6 ;       % thermal expansion coefficientTwoDeVec  = [ 1 1 0 ];    % -------------------------------------------------------------------------%                         DISCRETIZATION FEATURES % -------------------------------------------------------------------------%-ELEMENT STIFFNESS-a = 0.5 ;                                              b = 0.5;                                              h = 1 ;                    % element thicknessv0 = (2*a) * (2*b) * h ;   % volume of solid elementE0 = 2.1e5 ;               % elastic modulusEmin = 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