boxgirder2.m

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

kappa_concrete=1.81*[44 0; 0 44]; % W/K/m
rho_concrete = 2350;% kg/m^3
cv_concrete =0.22*rho_concrete;
h=0;
times =(0:2:24);
Ti = 19;
Ta_inner=[19, 18, 17, 18, 20, 26, 30, 30, 29, 26, 22, 20, 19];
Ta_lower=[19, 18, 17, 25, 28, 35, 32, 31, 32, 27, 22, 20, 19];
Ta_upper=[19, 18, 17, 35, 65, 80, 80, 60, 45, 35, 22, 20, 19];
Ta_vert=[19, 18, 17, 52, 65, 56, 35, 31, 32, 27, 22, 20, 19];
hs=[6.2, 6.2, 6, 6.3, 7.0, 7.3, 9.5, 13.4, 17.1, 17.6, 15.4, 9.5, 6.2];
% dt=2; % time 0
% tend= 96; % length of the time interval
t=16;
% theta = 1.0; % backward Euler
online_graphics=true;% plot the solution as it is computed?
mesh_size=0.2;%0.0372; %
[fens,gcells,groups,edge_gcells,edge_groups]=targe2_mesher({...
    'curve 1 line 0 0 2.60 0',...
    'curve 2 line 2.60 0 2.6 2.00',...
    'curve 3 line 2.6 2.0 4.85 2.0',...
    'curve 4 line 4.85 2 4.85 2.2',...
    'curve 5 line 4.85 2.2 0.0 2.2',...
    'curve 6 line 0.0 2.2 0.0 2',...
    'curve 7 line 0 2 2.1 2',...
    'curve 8 line 2.1 2 2.1 0.2',...
    'curve 9 line 2.1 0.2 0.0 0.2',...
    'curve 10 line 0 0.2 0 0',...
    ['subregion 1  property 1 boundary '...
    ' 1 2 3 4 5 6 7 8 9 10 '],...
    ['m-ctl-point constant ' num2str(mesh_size)]
    }, 1.0);
prop_concrete=...
    property_diffusion (struct('conductivity',kappa_concrete,...
    'specific_heat',cv_concrete,'source',0.0));
mater_concrete=mater_diffusion (struct('property',prop_concrete));
feb_concrete = feblock_diffusion (struct ('mater',mater_concrete,...
    'gcells',gcells(groups{1}),...
    'integration_rule',tri_rule(1)));
egcells_inner=edge_gcells([edge_groups{[(7:9)]}]);
egcells_lower=edge_gcells([edge_groups{[1, 3]}]);
egcells_upper=edge_gcells([edge_groups{[5]}]);
egcells_vert=edge_gcells([edge_groups{[2, 4]}]);
geom = field(struct('name',['geom'], 'dim', 2, 'fens',fens));
tempn=field(struct('name',['tempn'], 'dim', 1, 'nfens',...
    get(geom,'nfens')));
amb = clone(tempn, ['amb']);

tempn = numbereqns (tempn);
tempn = scatter_sysvec(tempn,gather_sysvec(tempn)*0+Ti);
% C = start (sparse_sysmat, get(tempn, 'neqns'));
% C = assemble (C, capacity(feb_concrete, geom, tempn));
% Cm = get(C,'mat');
K = start (sparse_sysmat, get(tempn, 'neqns'));
K = assemble (K, conductivity(feb_concrete, geom, tempn));
Kmc = get(K,'mat');

if online_graphics
    gv=graphic_viewer;
    gv=reset (gv,struct('limits',[0, 4.85, 0, 2.2, 0, 2.25]));
    dcm=data_colormap(struct('range',[20, 42],'colormap',jet));
end
minmaxT = []
actualt= 0;
% while actualt<tend+0.1*dt % Time stepping
    tempn1 = tempn;
    h=interp1(times,hs,t);
    efeb = feblock_diffusion (struct ('mater',mater_concrete,...
        'gcells',edge_gcells([edge_groups{[(1:5),(7:9)]}]),...
        'integration_rule',gauss_rule(1,1),...
        'surface_transfer', h));
    K = start (sparse_sysmat, get(tempn, 'neqns'));
    K = assemble (K, surface_transfer(efeb, geom, tempn));
    Km = Kmc+get(K,'mat');
    F = start (sysvec, get(tempn, 'neqns'));
    for i= 1:length(egcells_inner)
        conn = get(egcells_inner(i),'conn');
        Ta = interp1(times,Ta_inner,t);
        amb = set_ebc(amb, conn, conn*0+1, [], conn*0+Ta);
    end
    for i= 1:length(egcells_lower)
        conn = get(egcells_lower(i),'conn');
        Ta = interp1(times,Ta_lower,t);
        amb = set_ebc(amb, conn, conn*0+1, [], conn*0+Ta);
    end
    for i= 1:length(egcells_upper)
        conn = get(egcells_upper(i),'conn');
        Ta = interp1(times,Ta_upper,t);
        amb = set_ebc(amb, conn, conn*0+1, [], conn*0+Ta);
    end
    for i= 1:length(egcells_vert)
        conn = get(egcells_vert(i),'conn');
        Ta = interp1(times,Ta_vert,t);
        amb = set_ebc(amb, conn, conn*0+1, [], conn*0+Ta);
    end
    amb = apply_ebc (amb);
    F=assemble(F, surface_transfer_loads(efeb,geom,tempn,amb));
    Tn1 = Km\ get(F,'vec');
    tempn = scatter_sysvec(tempn1,Tn1);
    Tn=gather_sysvec(tempn);
    minmaxT= [minmaxT; min(Tn),max(Tn)]
%     t=t+dt; actualt= actualt+dt;
%     if (t >24), t=t-24; end
% end
    if online_graphics
        gv=reset (gv,struct('limits',[0, 4.85, 0, 2.2, 0, 2.25]));
        set(gca,'FontSize', 14)
        T=get(tempn,'values');
        colorfield=field(struct('name',['colorfield'],'data',...
            map_data(dcm, T)));
        geomT=field(struct ('name', ['geomT'], ...
            'data',[get(geom,'values'), 0.1*(get(tempn,'values')-Ti)]));
        for i=1:length (gcells)
            draw(gcells(i),gv,struct('x',geomT, 'u',0*geomT,...
                'colorfield',colorfield, 'edgecolor','none'));
            draw(gcells(i),gv,struct('x',geom, 'u',0*geom, ...
                'facecolor','none'));
        end
        camset(gv,[25.5545   11.6402   34.6031    1.9206    1.1491    1.5038   -0.1301    0.9681   -0.2140    6.0276])
        xlabel('X [m]');        ylabel('Y [m]');
        zlabel('Temperature [degrees C]')
        figure(gcf);
        title (['Time =' num2str(t)]); hold off; pause(0.1);
        %         saveas(gcf, ['shrinkfit-' num2str(t) '.png'], 'png');
   rotate(gv)
    end
% plot((1:size(minmaxT, 1))*dt,minmaxT,'linewidth', 3)
% set(gca,'FontSize', 14); grid on
% xlabel('Time in seconds')
% ylabel('Temperature in degrees Celsius')