多孔弹性材料中传播的膨胀波方法（Matlab代码实现）
👨🎓个人主页：研学社的博客
💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁
目录
💥1 概述
📚2 运行结果
🎉3 参考文献
🌈4 Matlab代码实现
💥1 概述
膨胀波是流体力学的基本概念之一，指流体中扰动区与未扰动区的分界面，流体通过此界面会压力降低。膨胀波问题中一个典型例子是超音速气流绕凸角的定常流动，这种流动称为普朗特-迈耶尔流动。在这种流动中，气流逐渐从一个方向转到另一个方向，同时不断加速。
📚2 运行结果
部分代码：
%%
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;
🎉3 参考文献
部分理论来源于网络，如有侵权请联系删除。
[1]刘占芳,严波,唐录成.饱和多孔弹性材料中加速度波的传播[J].重庆大学学报(自然科学版),1998(02):9-17.