t4nafems_convu.m

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

kappa=[52 0; 0 52]; % conductivity matrix
Q=0.0; % uniform heat source
h=750;
online_graphics = ~true;
num_integ_pts=[1,3];
options = [struct('quadratic',false),struct('quadratic',true)];
results = {[], []};
mesh_sizes = [0.2, 0.1, 0.05, 0.025];
for et= 1:2
    for mesh_size = mesh_sizes
        [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 ' num2str(mesh_size)]
            }, 1.0, options(et));
        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(et))));
        edgefeb = feblock_diffusion (struct ('mater',mater,...
            'gcells',edge_gcells([edge_groups{[2, 3, 4]}]),...
            'integration_rule',gauss_rule(1,4),...
            '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.0001)),...
            fenode_select(fens,struct('box',[0 1 1 1],...
            'inflate', 0.0001))]   ;
        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.0001))]    ;
        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'));
        results{et} =[results{et},gather(theta,fenode_select(fens,...
            struct('box',[0.6 0.6 0.2 0.2],'inflate', 0.0001)),'values')];

        % Plot
        if online_graphics
            gv=graphic_viewer;
            gv=reset (gv,[]);
            camset(gv, [-3.84897643467621  -4.70023392905539   4.85604797241542   0.20396556961523   0.39501248786488...
                0.29727878642496   0.35705972685897   0.44888559766667   0.81915204428899   6.42391124268001]);
            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',0.9));
                draw(gcells(i), gv, struct ('x',geom, 'u',0*geom, ...
                    'facecolor','none'));
            end
            xlabel('X [m]')
            ylabel('Y [m]')
            zlabel('Temperature [degrees C]')
        end
    end
    results{et}
end
xs =results{2};
xestim= 18.25396;
loglog(mesh_sizes,abs(results{2}-xestim)/xestim,'bo-','linewidth',3)
hold on
grid on
loglog(mesh_sizes,abs(results{1}-xestim)/xestim,'rs-','linewidth',3)