convcylqad.m

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

kappa=1.8*eye(2); % conductivity matrix
Q=0; % no heat source
h = 6;
radius = 1.5;
height =1.;
tol = (radius*height)/100000;
online_graphics = true;
mesh_sizes = 0.75.^(0:10);
% mesh_sizes = 2.2* 0.6.^(7:7);

Ts = [];
f=progressbar(' Calculating');
work0= sum (1./mesh_sizes.^2); work = work0;
for mesh_size = mesh_sizes
     [fens,gcells,groups,edge_gcells,edge_groups]=targe2_mesher({...
        ['curve 1 line 0 0 ' num2str(radius) ' 0'],...
        ['curve 2 line ' num2str(radius) ' 0 ' num2str(radius/0.75) ' ' num2str(height)],...
        ['curve 3 line ' num2str(radius/0.75) ' ' num2str(height) ' ' num2str(0) ' ' num2str(height) ],...
        ['curve 4 line ' num2str(0) ' ' num2str(height)  ' 0 0'],...
        'subregion 1  property 1 boundary 1 2 3 4',...
        ['m-ctl-point constant ' num2str(mesh_size)],...
         ['m-ctl-point 1 xy ' num2str(radius) ' 0 near ' num2str(mesh_size/100) ' influence ' num2str(mesh_size/30)],...
  
        }, 1.0, struct('axisymm', true,'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(3)));
    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));% everybody
    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', tol))];% hot face
    prescribed=ones(length(fenids),1);
    comp=[];
    val=zeros(length(fenids),1)+15;
    amb = set_ebc(amb, fenids, prescribed, comp, val);
   fenids=[fenode_select(fens,struct('box',[radius radius 0 0],...
        'inflate', tol))];% corner between the hot and cold face
    prescribed=ones(length(fenids),1);
    comp=[];
    val=zeros(length(fenids),1)+13.9;
    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])
        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
        view(2)
    end
    work= work-1/mesh_size^2; progressbar((work0-work)/work0,f)
end
progressbar(1.0,f);

c=['bo-';
    'c^-';
    'rx-';
    'gs-';
    'm*-';
    'kd-';];

ms = mesh_sizes; temps =Ts;
betas= []; xestims= [];
for i=1:min([length(c),length(Ts)-2])
%     ms
%     temps
    [xestim, beta] = richextrapol(temps(end-2:end),ms(end-2:end));
    xestims= [xestims,xestim]; betas= [betas,beta];
    loglog(ms,abs(temps-xestim),c(i,:)); grid on
    hold on;
    figure(gcf)
    ms=ms(1:end-1);
    temps=temps(1:end-1);
end
    xlabel('log(mesh size)')
    ylabel('log(|error in temperature|)')
Ts
xestims
betas
figure; plot(mesh_sizes,Ts,'bo-')