shockingss.m

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

kappa=[35.0]; % conductivity matrix
cm = 240.5;% specific heat per unit mass
rho=2500*1e-9;% mass density
cv =cm*1e3* rho;% specific heat per unit volume
h= 10e0;
Q=0; % uniform heat source
Tamb=20;
Tini=90;
L=1200;% thickness
dts=[10,5, 2.5, 1.25, 0.625]; % time step
theta = 1.0; % generalized trapezoidal method
online_graphics= true;% plot the solution as it is computed?
%igure;
n=4*5;% needs to be multiple of five
for dt=dts
tend=40; % length of the time interval
t=0;
    [fens,gcells] = block1d(L,n,1.0); % Mesh
    prop=property_diffusion(struct('conductivity',kappa,...
        'specific_heat',cv,'rho',rho,'source',Q));
    mater=mater_diffusion (struct('property',prop));
    feb = feblock_diffusion (struct (...
        'mater',mater,...
        'gcells',gcells,...
        'integration_rule',simpson_1_3_rule));
    point_gcells= [gcell_P1(struct('id',1,'conn',1));
        gcell_P1(struct('id',n+1,'conn',n+1))];
    pointfeb = feblock_diffusion (struct (...
        'mater',mater,...
        'gcells', point_gcells,...
        'integration_rule',point_rule,...
        'surface_transfer', h));
    geom = field(struct ('name',['geom'], 'dim', 1, 'fens',fens));
    tempn = field(struct ('name',['temp'], 'dim', 1,...
        'nfens',get(geom,'nfens')));
    tempn = numbereqns (tempn);
    tempn = scatter_sysvec(tempn,gather_sysvec(tempn)*0+Tini);
    amb = clone(tempn, ['amb']);
    for i= 1:length(point_gcells)
        conn = get(point_gcells(i),'conn');
        amb = set_ebc(amb, conn, conn*0+1, [], conn*0+Tamb);
    end
    amb = apply_ebc (amb);
    K = start (dense_sysmat, get(tempn, 'neqns'));
    K = assemble (K, conductivity(feb, geom, tempn));
    K = assemble (K, surface_transfer(pointfeb, geom, tempn));
    Km = get(K,'mat');
    C = start (dense_sysmat, get(tempn, 'neqns'));
    C = assemble (C, capacity(feb, geom, tempn));
    Cm = get(C,'mat');
    ts= []; sT= [];
    while t<tend+0.1*dt % Time stepping
        if online_graphics
            plot([0, 0],[0,Tini],'.'); hold on
            Ts=get(tempn,'values');
            ts= [ts,t]; sT= [sT,Ts(1)];
            plot(ts,sT,'r-', 'linewidth', 2);
            figure(gcf);
            title (['Time =' num2str(t)]); pause(0.2); %hold off;
        end
        tempn1 = tempn;
        F = start (sysvec, get(tempn, 'neqns'));
        F = assemble(F, surface_transfer_loads(pointfeb,geom,tempn,amb));
        Tn=gather_sysvec(tempn);
        Tn1 = (1/dt*Cm+theta*Km) \ ((1/dt*Cm-(1-theta)*Km)*Tn+...
            get(F,'vec'));
        tempn = scatter_sysvec(tempn1,Tn1);
        t=t+dt;
    end
end