bcsquare.m

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

kappa=[5.5 0; 0 0.2]; % conductivity matrix
Q=0.; % no heat source
num_integ_pts=1; 
magn = -1;
[fens,gcells,groups,edge_gcells,edge_groups]=targe2_mesher({...
    ['curve 1 line 0 0 48 0'],...
    ['curve 2 line 48 0 48 48'],...
    ['curve 3 line 48 48 0 48'],...
    ['curve 4 line 0 48 0 26'],...
    ['curve 5 line 0 26 0 13'],...
    ['curve 6 line 0 13 0 0'],...
    'subregion 1  property 1 boundary 1 2 3 4 5 6',...
    ['m-ctl-point constant 1.5']
    }, 1.0);
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)));
efeb = feblock_diffusion (struct ('mater',mater,...
    'gcells',edge_gcells([edge_groups{5}]),...
    'integration_rule',simpson_1_3_rule));
geom = field(struct('name',['geom'], 'dim', 2, 'fens',fens));
theta=field(struct('name',['theta'], 'dim', 1, 'nfens',...
    get(geom,'nfens')));
for e= [1, 2, 3, 4, 6]
    eg=edge_gcells([edge_groups{e}]);
    for i= 1:length(eg)
        conn = get(eg(i),'conn');
        theta = set_ebc(theta, conn, conn*0+1, [], conn*0);
    end
end
% for e= [5]
%     eg=edge_gcells([edge_groups{e}]);
%     for i= 1:length(eg)
%         conn = get(eg(i),'conn');
%         theta = set_ebc(theta, conn, conn*0+1, [], conn*0-10);
%     end
% end
theta = apply_ebc (theta);
theta = numbereqns (theta);
K = start (sparse_sysmat, get(theta, 'neqns'));
K = assemble (K, conductivity(feb, geom, theta));
F = start (sysvec, get(theta, 'neqns'));
fi= force_intensity(struct('magn',magn));
F = assemble (F, flux_loads(efeb, geom, theta, fi));
theta = scatter_sysvec(theta, get(K,'mat')\get(F,'vec'));

% Plot
gv=graphic_viewer;
gv=reset (gv,[]);
camset(gv,[ 154.6458  -70.7169  173.1645   24.5670   18.6838    3.9035   -0.6027   0.4142 0.6820   10.3396])
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'), get(theta,'values')]));
for i=1:length (gcells)
    draw(gcells(i), gv, struct ('x',geomT, 'u',0*geomT,...
        'colorfield',colorfield, 'shrink',1));
end