Time-varying boundary condition and impedance boundary condition for the wave equation

4 次查看（过去 30 天）
tri stan
tri stan 2018-2-20
I would like to model the propagation of an acoustic impulse in the duct which is opened on the right side. So my governing equation is wave equation:
m d2u/dt2 - div c grad u = f
where
m = 1/(rho*ss^2),
c = 1/rho,
f= 0,
where rho is the air density, and ss -- speed of sound.
The hard part for me is with proper imposing boundary conditions (BC):
  1. How to model a time-changing BC? In the script below I want to set u=1 for the first half of a second, then I would like to have zero Neumann condition. I'm trying to achieve this by using generalized Neumann condition
n*(c*grad(u)) + q*u = g
and setting q and g to very large value at first and then to zero (c = const = ~0.81).
fg = @(region,state) 100 * (state.time< 0.5) + 0 * (state.time >= 0.5) ;
fq = @(region,state) 100 * (state.time< 0.5) + 0 * (state.time >= 0.5) ;
Unfortunately, I read here that pde toolbox can't handle discontinuous functions with respect to time. In fact, the script is really slow and what is more it does not give expected results (experiments with "a bit more continuous" atan didn't succeed either, yet for the problem with q and g defined as linear functions of time the solution fulfills BC). Can it be fixed?
  1. I have found somewhere that impedance boundary condition should be imposed as follows:
$ \frac{1}{\rho} \frac{\partial u}{\partial n} = -\frac{1}{Z_0} \frac{\partial u}{\partial t}$
where
Z_0 = rho * ss
is the air characteristic impedance. Is it possible to implement this kind of BC in Matlab?
clear all
close all
% physics:
rho = 1.23 ; % [kg/m] air density
ss = 340 ; % [m/s] speed of sound
%%geometry & mesh
model = createpde();
R1 = [3,4,-1,1,1,-1,.4,.4,-.4,-.4]';
g = decsg(R1);
geometryFromEdges(model,g);
figure, pdegplot(model,'EdgeLabels','on'); axis equal
generateMesh(model);
figure, pdemesh(model); axis equal
%%boundary conditions
% _E1________du/dn=0_____________
% | E2
% u=1 |
% for a moment then air_impedance ??
% du/dn = 0 |
% | |
% E4_________du/dn=0___________E3_|
% top,bottom -- rigid wall
wallBC = applyBoundaryCondition(model,'Edge',[1 3], 'g', 0, 'q', 0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TODO ->
left -- impulse bc ?
fg = @(region,state) 100 * (state.time< 0.5) + 0 * (state.time >= 0.5) ;
fq = @(region,state) 100 * (state.time< 0.5) + 0 * (state.time >= 0.5) ;
% % % fg = @(region,state) 100*(-atan(100*(state.time-.5))/(pi/2)+1)/2 ;
% % % fq = @(region,state) 100*(-atan(100*(state.time-.5))/(pi/2)+1)/2 ;
% fg = @(region,state) 100*(state.time) ;
% fq = @(region,state) 100*(state.time) ;
impulseBC = applyBoundaryCondition(model,'Edge',4, 'g', fg, 'q', fq);
% right -- open space -- (should be) air characteristic acoustic impedance (rigid-wall BC at the moment)
air_impBC = applyBoundaryCondition(model,'Edge',[2], 'g', 0, 'q', 0);
%%%<- TODO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%initial conditions
u0 = @(locations) zeros(size(locations.x)) ;
ut0 = @(locations) zeros(size(locations.x)) ;
setInitialConditions(model,u0,ut0);
%%model coefficients
% $$ m \frac{\partial^2 u}{\partial t^2} -\nabla\cdot( c \nabla u) = f,$$
m = 1/(rho*ss^2) ;
c = 1/rho ;
f = 0;
specifyCoefficients(model,'m',m,'d',0,'c',c,'a',0,'f',f);
%%solution
n = 4; % number of frames in eventual animation
tlist = linspace(0,1,n) % list of times
tic
results = solvepde(model,tlist);
toc
%%Animate the solution.
u = results.NodalSolution;
umax = max(max(u));
umin = min(min(u));
figure;
M = moviein(n);
for i=1:n,
pdeplot(model,'XYData',u(:,i),'ZData',u(:,i),'ColorBar','off');
% colormap(jet)
% axis([-1 1 -1 1 umin umax]);
% caxis([umin umax]);
view([0 90])
M(:,i) = getframe;
end
%%Plot solution tlist = [0 0.33 0.66 1]
figure ;
for tt = 1:4
subplot(2,2,tt)
pdeplot(model,'XYData',u(:,tt),'ZData',u(:,tt),'ColorBar','off');
colormap(jet); view([0 90])
end

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 General PDEs 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by