actuator1.m

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

x0= 0;
x1=x0+5e-6;
x2=x1+10e-6;
x3=x2+10e-6;
x4=x3+10e-6;
y0=0;
y4=250e-6;
y3=y4-10e-6;
y2=y3-10e-6;
y1=y2-10e-6;
t=2e-6;
h=0.1e-6;
z0=0;
z3=2*t+h;
z2=z3-t;
z1=z2-h;
m1=2;
m2=2;
m3=2;
m4=3;
n1=20;
n2=4;
n3=2;
n4=2;
n5=7;
p1=1;
p2=1;
p3=1;
kappa=175*eye(3, 3); % Wconductivity matrix
DV=5;% voltage drop in volt
l =2*(y1+y2)/2;% length of the conductor
resistivity = 1.1e-5;% Ohm m
Q=DV^2/resistivity/l^2;% rate of Joule heating
T_substrate=293;% substrate temperature in degrees Kelvin

[fens,gcells] = mesh_hexahedron([x1,y0,z0;x2,y1,z1],m2,n1,p1);
[fens1,gcells1] = mesh_hexahedron([x1,y1,z0;x2,y2,z1],m2,n2,p1);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
[fens1,gcells1] = mesh_hexahedron([x0,y1,z0;x1,y2,z1],m1,n2,p1);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
[fens1,gcells1] = mesh_hexahedron([x0,y1,z1;x1,y2,z2],m1,n2,p2);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
[fens1,gcells1] = mesh_hexahedron([x0,y1,z2;x1,y2,z3],m1,n2,p3);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
[fens1,gcells1] = mesh_hexahedron([x0,y2,z2;x1,y3,z3],m1,n3,p3);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
[fens1,gcells1] = mesh_hexahedron([x0,y3,z2;x1,y4,z3],m1,n4,p3);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
[fens1,gcells1] = mesh_hexahedron([x1,y3,z2;x3,y4,z3],m4,n4,p3);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
[fens1,gcells1] = mesh_hexahedron([x3,y3,z2;x4,y4,z3],m3,n4,p3);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
[fens1,gcells1] = mesh_hexahedron([x3,y0,z2;x4,y3,z3],m3,n5,p3);
[fens,gcells1,gcells2] = merge_meshes(fens1, gcells1, fens, gcells, eps);
gcells= cat(2,gcells1,gcells2);
% gv=drawmesh({fens,gcells},'gcells','facecolor','red');
[fens,gcells] = H8_to_H20(fens,gcells);

hotprop=property_diffusion(struct('conductivity',kappa,'source',Q));
hotmater=mater_diffusion (struct('property',hotprop));
coldprop=property_diffusion(struct('conductivity',kappa,'source',0.0));
coldmater=mater_diffusion (struct('property',coldprop));
cl= gcell_select(fens, gcells, struct('box',[x0,x2,y0,y2,z0,z1],'inflate',t/100));
hotfeb = feblock_diffusion (struct ('mater',hotmater,...
    'gcells',gcells(cl),...
    'integration_rule',gauss_rule(3,2)));
coldfeb = feblock_diffusion (struct ('mater',coldmater,...
    'gcells',gcells(setdiff((1:length(gcells)),cl)),...
    'integration_rule',gauss_rule(3,2)));
geom = field(struct('name',['geom'], 'dim', 3, 'fens',fens));
theta=field(struct('name',['theta'], 'dim', 1, 'nfens',...
    get(geom,'nfens')));
fenids=fenode_select(fens,struct('box',[x0,x4,y0,y0,z0,z3],...
    'inflate', t/1000)) ; % prescribed temperature on substrate
prescribed=fenids*0+1;
comp=[];
val=zeros(length(fenids),1)+T_substrate;
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, cat (2,conductivity(hotfeb, geom, theta),...
    conductivity(coldfeb, geom, theta)));
F = start (sysvec, get(theta, 'neqns'));
F = assemble (F, source_loads(hotfeb, geom, theta));
F = assemble (F, cat (2,nz_ebc_loads_conductivity(hotfeb, geom, theta),...
    nz_ebc_loads_conductivity(coldfeb, geom, theta)));
theta = scatter_sysvec(theta, get(K,'mat')\get(F,'vec'));

y_i=gather(geom,fenode_select(fens,struct('box',[x1,x1,y0,y1,z1,z1],...
    'inflate', t/1000)),'values','noreshape');
T_i=gather(theta,fenode_select(fens,struct('box',[x1,x1,y0,y1,z1,z1],...
    'inflate', t/1000)),'values');
[ignore,ix]=sort(y_i(:,2));
figure;
plot(y_i(ix,2),T_i(ix),'bd','linewidth',3); hold on
y_o=gather(geom,fenode_select(fens,struct('box',[x3,x3,y0,y3,z2,z2],...
    'inflate', t/1000)),'values','noreshape');
T_o=gather(theta,fenode_select(fens,struct('box',[x3,x3,y0,y3,z2,z2],...
    'inflate', t/1000)),'values');
[ignore,ix]=sort(y_o(:,2));
plot(y_o(ix,2),T_o(ix),'bo','linewidth',3); hold off

% Plot
gv=graphic_viewer;
gv=reset (gv,[]);
T=get(theta,'values');
cmap= hot;
dcm=data_colormap(struct('range',[min(T),max(T)],'colormap',cmap));
colorfield=field(struct ('name', ['colorfield'], 'data',...
    map_data(dcm, T)));
draw(hotfeb, gv, struct ('x',geom, 'u',0*geom,...
    'colorfield',colorfield, 'shrink',1));

draw(coldfeb, gv, struct ('x',geom, 'u',0*geom,...
    'colorfield',colorfield, 'shrink',1));
draw_colorbar(gv, struct('colormap',cmap,'label','Temperature',...
    'minmax',[min(T),max(T)]))
rotate(gv);