squareinsquare.m

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

kappainner=[2.25 0; 0 0.06]; % orthotropic conductivity matrix
kappaouter=[0.25 0; 0 0.25]; % isotropic conductivity matrix
alpha =-45;% local material orientation angle
ca=cos(2*pi/360*alpha); sa=sin(2*pi/360*alpha);
Rm = [ca, -sa;sa, ca];% local material directions
num_integ_pts=1; % 1-point quadrature
[fens,gcells, groups] = targe2_mesher({...
    ['curve 1 line -48 -48 48 -48'],...
    ['curve 2 line 48 -48 48 48'],...
    ['curve 3 line 48 48 -48 48'],...
    ['curve 4 line -48 48 -48 -48'],...
    ['curve 5 line 0 -31 31 0'],...
    ['curve 6 line 31 0 0 31'],...
    ['curve 7 line 0 31 -31 0'],...
    ['curve 8 line -31 0 0 -31'],...
    ['subregion 1  property 1 ' ...
    '   boundary 1 2 3 4 -8 -7 -6 -5'],...
    ['subregion 2  property 2 '...
    '   boundary 5 6 7 8'],...
    ['m-ctl-point constant 4.75']
    }, 1.0);
propinner = property_diffusion(struct('conductivity',kappainner,...
    'source',0));
materinner=mater_diffusion(struct('property',propinner));
febinner = feblock_diffusion(struct ('mater',materinner,...
    'gcells',gcells(groups{2}),...
    'integration_rule',tri_rule(num_integ_pts),'Rm',Rm));
propouter = property_diffusion(struct('conductivity',kappaouter,...
    'source',0));
materouter=mater_diffusion(struct('property',propouter));
febouter = feblock_diffusion(struct ('mater',materouter,...
    'gcells',gcells(groups{1}),...
    'integration_rule',tri_rule(num_integ_pts)));
geom = field(struct('name',['geom'], 'dim', 2, 'fens',fens));
theta=field(struct('name',['theta'], 'dim', 1, 'nfens',get(geom,'nfens')));
fenids=[fenode_select(fens,struct('box',[-48 48 -48 -48],...
    'inflate', 0.01))];
prescribed=ones(length(fenids),1);
comp=[];
val=zeros(length(fenids),1)+20;% ambient temperature
theta = set_ebc(theta, fenids, prescribed, comp, val);
fenids=[fenode_select(fens,struct('box',[-48 48 48 48],...
    'inflate', 0.01))];
prescribed=ones(length(fenids),1);
comp=[];
val=zeros(length(fenids),1)+57;% hot inner surface
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(febinner, geom, theta));
K = assemble (K, conductivity(febouter, geom, theta));
F = start (sysvec, get(theta, 'neqns'));
F = assemble (F, nz_ebc_loads_conductivity(febouter, geom, theta));
theta = scatter_sysvec(theta, get(K,'mat')\get(F,'vec'));

% Plothot
gv=graphic_viewer;
gv=reset (gv,[]);
T=get(theta,'values');
dcm=data_colormap(struct ('range',[min(T),max(T)], 'colormap',hot(9)));
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));
    draw(gcells(i), gv, struct ('x',geom, 'u',0*geom, 'facecolor','none'));
end
view(2)