convcylflux.m

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

kappa=1.8*eye(2); % conductivity matrix
Q=0; % no heat source
h = 6;
radius = 1.5;
height =1.;
online_graphics = true;
% mesh_sizes = 2.2* 0.6.^(2:9);
mesh_sizes = 2.2* 0.6.^(5:5);

Ts = [];
f=progressbar(' Calculating');
work0= sum (1./mesh_sizes.^2); work = work0;
for mesh_size = mesh_sizes
    [fens,gcells] = mesh_quadrilateral([[0, 0];...
        [radius, 0]; [radius/0.75, height];...
        [0, height]],...
        radius/mesh_size,height/mesh_size,...
        struct('axisymm', 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',gauss_rule(2, 2)));
    edge_gcells=mesh_bdry(gcells, struct('axisymm', true));
    el=[gcell_select(fens, edge_gcells, struct('box', [0, radius, 0, 0],'inflate', 0.0001)), ...
        gcell_select(fens, edge_gcells, struct('box', [0.001, 100*radius, 0, height],'inflate', 0.0001)),...
        gcell_select(fens, edge_gcells, struct('box', [0, 100*radius, height, height],'inflate', 0.0001))];
    efeb = feblock_diffusion (struct ('mater',mater,...
        'gcells',edge_gcells(unique(el)),...
        'integration_rule',gauss_rule(1,2),...
        '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']);
    theta = numbereqns (theta);
    K = start (sparse_sysmat, get(theta, 'neqns'));
    K = assemble (K, conductivity(feb, geom, theta));
    K = assemble (K, surface_transfer(efeb, geom, theta));
    F = start (sysvec, get(theta, 'neqns'));
    fenids=(1:length(fens));
    prescribed=ones(length(fenids),1);
    comp=[];
    val=zeros(length(fenids),1)+13;
    amb = set_ebc(amb, fenids, prescribed, comp, val);
    fenids=[fenode_select(fens,struct('box',[0 radius 0 0],...
        'inflate', 0.01))];
    prescribed=ones(length(fenids),1);
    comp=[];
    val=zeros(length(fenids),1)+15;
    amb = set_ebc(amb, fenids, prescribed, comp, val);
    amb = apply_ebc (amb);
    F = assemble (F, surface_transfer_loads(efeb, geom, theta, amb));
    F = assemble(F, nz_ebc_loads_conductivity(feb, geom, theta));
    F = assemble(F, nz_ebc_loads_surface_transfer(efeb, geom, theta));
    sol=get(K,'mat')\get(F,'vec');
    theta = scatter_sysvec(theta, sol);
    % energies=[energies,sol'*get(K,'mat')*sol];
    T_E=gather(theta,fenode_select(fens,struct('box',[radius, radius, 0, 0],...
        'inflate', 0.0001)),'values');
    Ts=[Ts,T_E]

    if online_graphics %  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])
        for i=1:length (gcells)
            draw(gcells(i), gv, struct ('x',geom, 'u',0*geom,...
                'facecolor','none', 'shrink',1));
        end
        draw_integration_points(feb, gv, struct ('x',geom,'u',0*geom, ...
            'theta', theta, 'scale', 0.05,'color','red'));
        view(2)
    end
    work= work-1/mesh_size^2; progressbar((work0-work)/work0,f)
end
delete(f);