t4nafemsq.m

在有限元计算中作热交换分析的算例

kappa=[52 0; 0 52]; % conductivity matrix
Q=0.0; % uniform heat source
h=750;
num_integ_pts=3;
[fens,gcells,groups,edge_gcells,edge_groups]=targe2_mesher({... 
        ['curve 1 line 0 0 0.6 0'],...
        ['curve 2 line 0.6 0 0.6 0.2'],...
        ['curve 3 line 0.6 0.2 0.6 1.0'],...
        ['curve 4 line 0.6 1.0 0 1.0'],...
        ['curve 5 line 0 1.0 0 0'],...
        'subregion 1  property 1 boundary 1 2 3 4 5',...
        ['m-ctl-point constant 0.13']
        }, 1.0, struct('quadratic',true));
prop=property_diffusion(struct('conductivity',kappa,'source',Q));
mater=mater_diffusion (struct('property',prop));
feb = feblock_diffusion (struct ('mater',mater,...
    'gcells',gcells,...
    'integration_rule',tri_rule(num_integ_pts)));
edgefeb = feblock_diffusion (struct ('mater',mater,...
    'gcells',edge_gcells([edge_groups{[2, 3, 4]}]),...
    'integration_rule',gauss_rule(1,num_integ_pts),...
    'surface_transfer', h));
geom = field(struct('name',['geom'], 'dim', 2, 'fens',fens));
theta=field(struct('name',['theta'], 'dim', 1, 'nfens',...
    get(geom,'nfens')));
amb = clone(theta, ['amb']);
fenids=[
    fenode_select(fens,struct('box',[0.6 0.6 0 1],...
    'inflate', 0.01)),...
    fenode_select(fens,struct('box',[0 1 1 1],...
    'inflate', 0.01))]   ; 
prescribed=ones(length(fenids),1);
comp=[];
val=zeros(length(fenids),1)+0.0;
amb = set_ebc(amb, fenids, prescribed, comp, val);
amb = apply_ebc (amb);
fenids=[
    fenode_select(fens,struct('box',[0. 0.6 0 0],...
    'inflate', 0.01))]    ;
prescribed=ones(length(fenids),1);
comp=[];
val=zeros(length(fenids),1)+100.0;
theta = set_ebc(theta, fenids, prescribed, comp, val);
theta = apply_ebc (theta);
theta = numbereqns (theta);
K = start (sparse_sysmat, get(theta, 'neqns'));
K = assemble (K, conductivity(feb, geom, theta));
K = assemble (K, surface_transfer(edgefeb, geom, theta));
F = start (sysvec, get(theta, 'neqns'));
F = assemble(F, source_loads(feb, geom, theta));
F = assemble(F, nz_ebc_loads_conductivity(feb, geom, theta));
F = assemble(F, surface_transfer_loads(edgefeb,geom,theta,amb));
theta = scatter_sysvec(theta, get(K,'mat')\get(F,'vec'));
gather(theta,fenode_select(fens,...
    struct('box',[0.6 0.6 0.2 0.2],'inflate', 0.01)),'values')

% Plot
gv=graphic_viewer;
gv=reset (gv,[]);
set(gca,'FontSize', 14)
T=get(theta,'values');
dcm=data_colormap(struct('range',[min(T),max(T)],'colormap',jet));
colorfield=field(struct ('name', ['colorfield'], 'data',...
    map_data(dcm, T)));
geomT=field(struct ('name', ['geomT'], ...
    'data',[get(geom,'values'), 0.01*get(theta,'values')]));
for i=1:length (gcells)
    draw(gcells(i), gv, struct ('x',geomT, 'u',0*geomT,...
        'colorfield',colorfield, 'shrink',1));
    draw(gcells(i), gv, struct ('x',geom, 'u',0*geom, ...
        'facecolor','none'));
end
xlabel('X [m]')
ylabel('Y [m]')
zlabel('Temperature [degrees C]')
% for i=1:length (fens)
%     draw(fens(i), gv, struct ('x',geomT, 'u',0*geomT));
% end