多孔弹性材料中传播的膨胀波方法（Matlab代码实现）

%%
clear
%% inputs
freq=10^3;%*(2.^([-20:20]/3)); %frequency
%% setup domain and mesh
domain.dim=[0.05 0.025]; %dimensions of domain
domain.off=[0 0]; %offset of domain
domain.fce=[1 2 3 4]; %face ids (for boundary conditions) on [-x +x -y +y]
nnd=domain.dim./0.0025+1; %number of nodes in each direction
mesh=blockmesh(nnd,domain.dim,domain.off,domain.fce);
%% loop for frequency
for f=1:length(freq)
fprintf(['computing frequency ',num2str(f),' of ',num2str(length(freq)),' ... \n']);
%% setup boundary conditions
addpath('MESH')
bcs=initbcs(mesh.nodes,mesh.elems,mesh.faces,8);
bcs=addbcface(bcs,1,1,5);
bcs=addbcface(bcs,1,2,6);
bcs=addbcface(bcs,1,3,7);
bcs=addbcface(bcs,1,4,8);
bcs=addbcface(bcs,1,1,9);
% bcs=addbcface(bcs,1,2,10);
bcs=addbcface(bcs,1,3,11);
bcs=addbcface(bcs,1,4,12);
bcs.type{1}(:)=2;
bcs.coef{1}(bcs.face{1}==1,5)=1; %face 1
bcs.coef{1}(bcs.face{1}==2,1)=1; %face 2
bcs.coef{1}(bcs.face{1}==3,2)=1; %face 3
bcs.coef{1}(bcs.face{1}==4,2)=1; %face 4
bcs.coef{1}(bcs.face{1}==5,6)=1; %face 5
bcs.coef{1}(bcs.face{1}==6,3)=1; %face 6
bcs.coef{1}(bcs.face{1}==7,4)=1; %face 7
bcs.coef{1}(bcs.face{1}==8,4)=1; %face 8
bcs.coef{1}(bcs.face{1}==9,8)=1; bcs.rhds{1}(bcs.face{1}==9)=1; %face 9
% bcs.coef{1}(bcs.face{1}==10,1)=1; %face 10
bcs.coef{1}(bcs.face{1}==11,6)=1; %face 11
bcs.coef{1}(bcs.face{1}==12,6)=1; %face 12
%% setup physics
addpath('PLANES')
air=air_properties_generic;
medium=4003;
PEM.name_mat=['Mat_porous_' num2str(medium-1000*floor(medium/1000))];
PEM.typ_mat=floor(medium/1000);
eval(['PEM=Mat_porous_' num2str(medium-1000*floor(medium/1000)),'(PEM);'])
PEM=properties_JCA(PEM,air,freq(f));
PEM=properties_PEM(PEM,air,freq(f));
physics=PEM;
%% assemble
[stiff,force,err,tol,condH]=assemble(mesh,bcs,physics,freq(f));
%% solve
fprintf('solving ... \n');
operationtime=cputime;
fprintf(' solving directly ... ');
q=full(stiff\force);
fprintf('done\n');
%% postprocess
nnds=mesh.nnds;
soln.usx(:,f)=q((1:nnds)+nnds*0);
soln.usy(:,f)=q((1:nnds)+nnds*1);
soln.utx(:,f)=q((1:nnds)+nnds*2);
soln.uty(:,f)=q((1:nnds)+nnds*3);
soln.sxx(:,f)=q((1:nnds)+nnds*4);
soln.sxy(:,f)=q((1:nnds)+nnds*5);
soln.syy(:,f)=q((1:nnds)+nnds*6);
soln.p(:,f)=q((1:nnds)+nnds*7);
% soln.cond(:,f)=condest(stiff);
us=mean(reshape(soln.usx(:,f),nnd(1),nnd(2)),2);
ut=mean(reshape(soln.utx(:,f),nnd(1),nnd(2)),2);
s=mean(reshape(soln.sxx(:,f),nnd(1),nnd(2)),2);
p=mean(reshape(soln.p(:,f),nnd(1),nnd(2)),2);
Z=p(1)/ut(1);
R=(Z-air.rho*air.c)/(Z+air.rho*air.c);
alpha(:,f)=1-abs(R)^2;