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

       目录

💥1 概述

📚2 运行结果

🎉3 参考文献

👨‍💻4 Matlab代码

💥1 概述

膨胀波是流体力学的基本概念之一，指流体中扰动区与未扰动区的分界面，流体通过此界面会压力降低。膨胀波问题中一个典型例子是超音速气流绕凸角的定常流动，这种流动称为普朗特-迈耶尔流动。在这种流动中，气流逐渐从一个方向转到另一个方向，同时不断加速。

📚2 运行结果

🎉3 参考文献

[1]刘占芳,严波,唐录成.饱和多孔弹性材料中加速度波的传播[J].重庆大学学报(自然科学版),1998(02):9-17.

👨‍💻4 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;
  
  x=mesh.nodes{1}(1:nnd(1),1);
  addpath('ANALYTICAL')
  ana=analytical_biot(freq(f),x);
  x0=linspace(min(x),max(x),201);
  ana0=analytical_biot(freq(f),x0);
  Z=ana.p(1)/ana.ut(1);
  R=(Z-air.rho*air.c)/(Z+air.rho*air.c);
  alphaana(:,f)=1-abs(R)^2;  
  
  Err(1,f)=sqrt(sum(abs(us-ana.us).^2))./sqrt(sum(abs(ana.us.^2)));
  Err(2,f)=sqrt(sum(abs(ut-ana.ut).^2))./sqrt(sum(abs(ana.ut.^2)));
  Err(3,f)=sqrt(sum(abs(s-ana.s).^2))./sqrt(sum(abs(ana.s.^2)));
  Err(4,f)=sqrt(sum(abs(p-ana.p).^2))./sqrt(sum(abs(ana.p.^2)));
  
  PinvRes(:,f)=sqrt(sum(err.^2));
  Tol(:,f)=mean(tol);
  CondH(:,f)=mean(condH);