alusteel.m

有限元热计算的算例

disp('Aluminum/steel assembly');

% Parameters:
Ea=70e9;
nua=0.33;
alphaa= 23e-6;
Es=200e9;
nus=0.3;
alphas= 12e-6;
thickness = 5.0;
integration_order =1;
scale = 100;
h=5;
% Mesh
[fens,gcells,groups,edge_gcells,edge_groups]=targe2_mesher({...
    'curve 1 line 0 0 75 0',...
    'curve 2 line 75 0 75 40',...
    'curve 3 line 75 40 50 40',...
    'curve 4 line 50 40 50 15',...
    'curve 5 line 50 15 0 15',...
    'curve 6 line 0 15 0 0',...
    'curve 7 line 50 40 0 40',...
    'curve 8 line 0 40 0 15',...
    ['subregion 1  property 1 boundary '...
    ' 1 2 3 4 5 6'],...
    ['subregion 2  property 2 boundary '...
    ' -5 -4 7 8'],...
    ['m-ctl-point constant ' num2str(h)],...
    ['m-ctl-point 1 xy 50 15 near ' num2str(h/10)...
    ' influence  ' num2str(h/4)]
    }, thickness);
% Material
propa = property_linel_iso ...
    (struct('E',Ea,'nu',nua,'alpha', alphaa));
matera = mater_defor_ss_linel_biax (struct('property',propa, ...
    'reduction','stress'));
props = property_linel_iso ...
    (struct('E',Es,'nu',nus,'alpha', alphas));
maters = mater_defor_ss_linel_biax (struct('property',props, ...
    'reduction','stress'));
% Finite element block
feba = feblock_defor_ss (struct ('mater',matera,...
    'gcells',gcells(groups{2}),...
    'integration_rule',tri_rule (integration_order)));
febs = feblock_defor_ss (struct ('mater',maters,...
    'gcells',gcells(groups{1}),...
    'integration_rule',tri_rule (integration_order)));
% Geometry
geom = field(struct ('name',['geom'], 'dim', 2, 'fens',fens));
% Define the displacement field
u   = clone(geom,'u');
u   = u*0; % zero out
% Apply EBC's
ebc_fenids=fenode_select (fens,struct('box',[0 0 0 40],'inflate', 0.01));
ebc_prescribed=ones(1,length (ebc_fenids));
ebc_comp=1*ones(1,length (ebc_fenids));
ebc_val=ebc_fenids*0;
u   = set_ebc(u, ebc_fenids, ebc_prescribed, ebc_comp, ebc_val);
ebc_fenids=fenode_select (fens,struct('box',[0 75 40 40],'inflate', 0.01));
ebc_prescribed=ones(1,length (ebc_fenids));
ebc_comp=2*ones(1,length (ebc_fenids));
ebc_val=ebc_fenids*0;
u   = set_ebc(u, ebc_fenids, ebc_prescribed, ebc_comp, ebc_val);
u   = apply_ebc (u);
% Number equations
u   = numbereqns (u);
% Temperature field
dT = field(struct ('name',['dT'], 'dim', 1, ...
    'data',zeros(length(fens),1)+70));
% Assemble the system matrix
K = start (sparse_sysmat, get(u, 'neqns'));
K = assemble (K, cat(2,stiffness(febs, geom, u),stiffness(feba, geom, u)));
% Load
F = start (sysvec, get(u, 'neqns'));
F = assemble (F, cat(2,thermal_strain_loads(feba, geom, u, dT),...
    thermal_strain_loads(febs, geom, u, dT)));
% Solve
u = scatter_sysvec(u, get(K,'mat')\get(F,'vec'));
% get(u,'values')

% Plot
gv=graphic_viewer;
gv=reset (gv,struct ('limits', [0, 100, -8, 40]));
set(gca,'FontSize', 12)
cmap=jet;
cmpn=3;
flda = field_from_integration_points(feba, geom, u, dT, 'Cauchy', cmpn);
flds = field_from_integration_points(febs, geom, u, dT, 'Cauchy', cmpn);
nvalsa=get(flda,'values');
nvalss=get(flds,'values');
nvalmin =min(min(nvalsa),min(nvalss));
nvalmax =max(max(nvalsa),max(nvalss));
dcm=data_colormap(struct ('range',[nvalmin,nvalmax], 'colormap',cmap));
colorfield=field(struct ('name', ['colorfield'], 'data',map_data(dcm, nvalsa)));
draw(feba, gv, struct ('x',geom,'u', scale*u, 'colorfield',colorfield, 'shrink',1));
% draw(feba, gv, struct ('x',geom,'u', 0*u, 'facecolor',' none'));
colorfield=field(struct ('name', ['colorfield'], 'data',map_data(dcm, nvalss)));
draw(febs, gv, struct ('x',geom,'u', scale*u, 'colorfield',colorfield, 'shrink',1));
% draw(febs, gv, struct ('x',geom,'u', 0*u, 'facecolor',' none'));
lighting  none;
title (['h=' num2str(h)]);
colormap(cmap);
cbh=colorbar;
set(cbh,...
    'Position',[0.72 0.33 0.05 0.5],...
    'YLim',[0,1],...
    'YTick',[0,1],...
    'YTickLabel',{[num2str(nvalmin)],[num2str(nvalmax)]},...
    'FontSize', 12);
set(get(cbh,'XLabel'),'String','\sigma_{xy}');
view (2)
 saveas(gcf, [mfilename '-' num2str(h) '.png'], 'png');