two_circles_assembly.m

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

kappa=[0.2 0; 0 0.2]; % conductivity matrix
Q=0.1550; %  heat source
Ta= 1.3;
online_graphics= true;
h=0.061;
[fens,gcells,groups,edge_gcells,edge_groups]=targe2_mesher({...
    'curve 1 circle  Center 0 0 radius 5',...
    'curve 2 circle Center 0.5 -1.5 radius 2',...
    ['subregion 1  property 1 boundary 1 hole -2 '],...
    ['m-ctl-point constant 2.8']
    }, 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(1)));
efeb = feblock_diffusion (struct ('gcells',edge_gcells([edge_groups{1}]),...
    'integration_rule',gauss_rule(1,1),...
    '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']);
for i= 1:length(edge_gcells)
    conn = get(edge_gcells(i),'conn');
    amb = set_ebc(amb, conn, conn*0+1, [], conn*0+Ta);
end
amb = apply_ebc (amb);

theta = numbereqns (theta);

K = start (sparse_sysmat, get(theta, 'neqns'));
ems= conductivity(feb, geom, theta);
f=ars6(get(K,'mat'),'K',struct('xy',[1, 200],'maxcol',get(theta, 'neqns'),'maxrow',get(theta, 'neqns'),'butwidth',0.6));
pause (2);
for i=1:length(ems)
    get(ems(i),'eqnums')'
    pause (5);
    K = assemble (K,ems(i));
    close(f);
    f=ars6(get(K,'mat'),'K',struct('xy',[1, 200],'maxcol',get(theta, 'neqns'),'maxrow',get(theta, 'neqns'),'butwidth',0.6));
end
close(f);
ems = surface_transfer(efeb, geom, theta);
f=ars6(get(K,'mat'),'K',struct('xy',[1, 200],'maxcol',get(theta, 'neqns'),'maxrow',get(theta, 'neqns'),'butwidth',0.6));
pause (2);
for i=1:length(ems)
    get(ems(i),'eqnums')'
    pause (5);
    K = assemble (K,ems(i));
    close(f);
    f=ars6(get(K,'mat'),'K',struct('xy',[1, 200],'maxcol',get(theta, 'neqns'),'maxrow',get(theta, 'neqns'),'butwidth',0.6));
end
close(f);
F = start (sysvec, get(theta, 'neqns'));
F = assemble (F, source_loads(feb, geom, theta));
F = assemble(F, surface_transfer_loads(efeb,geom,theta,amb));

theta = scatter_sysvec(theta,get(K,'mat')\get(F,'vec'));

if online_graphics
    gv=graphic_viewer;
    gv=reset (gv,[]);
    set(gca,'FontSize', 14)
    dcm=data_colormap(struct('range',...
        [min(get(theta,'values')),max(get(theta,'values'))],...
        'colormap',jet));
    colorfield=field(struct('name',['colorfield'],'data',...
        map_data(dcm, get(theta,'values'))));
    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.0));
        draw(gcells(i),gv,struct('x',geom, 'u',0*geom, ...
            'facecolor','none'));
    end
    xlabel('X [m]');        ylabel('Y [m]');
    zlabel('Temperature [degrees C]')
    figure(gcf);
    title (['Temperature']);
    %         saveas(gcf, ['shrinkfit-' num2str(t) '.png'], 'png');
    rotate(gv)
end