[READ-ONLY] Mirror of https://github.com/jmrplens/FDTDexamples. Modelos FDTD de aeroacústica simples y complejos para MATLAB 2019b o superior
acoustic
fdtd
matlab
Upload
Jose M. Requena Plens
(Sep 29, 2019, 4:17 AM +0200)
7c6c8728
59b279f0
+178
-5
+4
-5
FDTD_2D_Basico.m
+4
-5
FDTD_2D_Basico.m
···
43
43
ux = zeros(nx+1,ny); % Matriz de velocidad de particula en x
44
44
uy = zeros(nx,ny+1); % Matriz de velocidad de particula en y
45
45
46
-
%% POSICIONES DE LA FUENTE
47
-
%
46
+
%% POSICION DE LA FUENTE
48
47
posNy = round((Origen(2))/dh) + nPML;
49
48
posNx = round(Origen(1)/dh) + nPML;
50
49
51
50
52
51
%% EXCITACI�N
53
52
lenT = ts*10^-3/dt; % Longitud del vector de tiempo
54
-
a = 2000/(sqrt(pi)/2)*4;
55
-
t=((1:lenT)/(1/dt)-4/a); % Vector de tiempos
56
-
w=-(exp(-a^2*(t.^2)/2).*(a^2*(t.^2)-1)); % Ricker
53
+
a = 2000/(sqrt(pi)/2)*4;
54
+
t = ((1:lenT)/(1/dt)-4/a); % Vector de tiempos
55
+
w = -(exp(-a^2*(t.^2)/2).*(a^2*(t.^2)-1)); % Ricker
57
56
58
57
%% PML
59
58
% Gradiente de impedancias desde el valor del medio hasta un porcentaje de
+174
FDTD_QRD_Unidimensional_2D.m
+174
FDTD_QRD_Unidimensional_2D.m
···
1
+
%% =========================================================================
2
+
%
3
+
% Modelo FDTD de difusor QRD en recinto 2D con excitaci�n Ricker o
4
+
% continua
5
+
%
6
+
%
7
+
% Jose Manuel Requena Plens (29/09/2019)
8
+
%
9
+
% =========================================================================
10
+
clc, clear, close all;
11
+
12
+
%% GENERAL
13
+
% Caracteristicas del medio
14
+
rho = 1.21; % Densidad del medio [kg/m3]
15
+
c = 341; % Velocidad de propagaci�n [m/s]
16
+
k = (c^2)*rho; % Modulo de compresibilidad
17
+
% Frecuencia de emision
18
+
freq = 3000; % (Hz)
19
+
% Tipo de excitaci�n ('pulse' o 'CW')
20
+
excit = 'pulse';
21
+
% Dise�o del difusor QRD
22
+
N = 7;
23
+
fd = 1000;
24
+
fw = 5*fd;
25
+
% �Difusor o Panel plano?
26
+
tipo = 'difusor'; % 'difusor' o 'plano'
27
+
% Dimensiones del mapa en metros
28
+
lx = 1; % [m]
29
+
ly = 1; % [m]
30
+
% Posicion de la fuente
31
+
Origen = [0.5,0.99]; % [m,m]
32
+
% Tama�o/Ancho de la PML en metros
33
+
lPML = 0.5; % m
34
+
% Par�metros de calculo
35
+
dh = .01; % Definici�n espacial [m]
36
+
dt = dh/341/2; % Definici�n temporal [m]
37
+
ts = 7; % Tiempo de simulaci�n [m]
38
+
39
+
%% MATRICES
40
+
% Dimensiones de la matriz del mapa
41
+
nPML = round(lPML/dh);
42
+
nx = round(lx/dh)+nPML*2; % Por 2 ya que la PML esta arriba y abajo
43
+
ny = round(ly/dh)+nPML*2; % Por 2 ya que la PML esta a izquierda y derecha
44
+
% Crear matriz del mapa
45
+
p = zeros(nx,ny); % Matrices de presion
46
+
px = zeros(nx,ny);
47
+
py = zeros(nx,ny);
48
+
ux = zeros(nx+1,ny); % Matriz de velocidad de particula en x
49
+
uy = zeros(nx,ny+1); % Matriz de velocidad de particula en y
50
+
51
+
%% POSICION DE LA FUENTE
52
+
posNy = round((Origen(2))/dh) + nPML;
53
+
posNx = round(Origen(1)/dh) + nPML;
54
+
55
+
%% EXCITACI�N
56
+
lenT = ts*10^-3/dt; % Longitud del vector de tiempo
57
+
a = freq/(sqrt(pi)/2)*4;
58
+
t = ((1:lenT)/(1/dt)-4/a); % Vector de tiempos
59
+
switch excit
60
+
case 'pulse'
61
+
w = -(exp(-a^2*(t.^2)/2).*(a^2*(t.^2)-1)); % Ricker
62
+
case 'CW'
63
+
w = cos(2*pi*freq*t);
64
+
end
65
+
66
+
%% PML
67
+
% Gradiente de impedancias desde el valor del medio hasta un porcentaje de
68
+
% este
69
+
gammamax = 0.5; % Maxima reducci�n de la impedancia / porcentaje
70
+
% PML Izquierda y derecha
71
+
gammaux = zeros(nx+1,ny);
72
+
gammaux(1:nPML,:) = repmat(gammamax*((nPML:-1:1)'/nPML).^2,1,ny); % Izquierda
73
+
gammaux(1+end-nPML:end,:) = repmat(gammamax*((1:1:nPML)'/nPML).^2,1,ny); % Derecha
74
+
gammax = (gammaux(1:end-1,:)+gammaux(2:end,:))/2; % Conjunto
75
+
% PML Superior e inferior
76
+
gammauy = zeros(nx,ny+1);
77
+
gammauy(:,1:nPML) = repmat(gammamax*((nPML:-1:1)/nPML).^2,nx,1); % Inferior
78
+
gammauy(:,1+end-nPML:end) = repmat(gammamax*((1:1:nPML)/nPML).^2,nx,1); % Superior
79
+
gammay = (gammauy(:,1:end-1)+gammauy(:,2:end))/2; % Conjunto
80
+
81
+
%% Condiciones de contorno
82
+
ux(1,:) = -p(1,:)/rho/c;
83
+
ux(end,:) = p(end,:)/rho/c;
84
+
uy(:,1) = -p(:,1)/rho/c;
85
+
uy(:,end) = p(:,end)/rho/c;
86
+
87
+
%% PAR�METROS DIFUSOR QRD
88
+
% Dimensiones
89
+
lambdad = c/fd; % Longitud de onda (profundidad) [m]
90
+
lambdaw = c/fw; % Longitud de onda (ancho) [m]
91
+
wide = lambdaw/2; % Ancho de cavidad [m]
92
+
widedh = round(wide/dh);
93
+
% Vector de profundidades
94
+
n0 = ceil(N/2); % Offset para centrar las secuencias
95
+
n = ceil((1:N))-(n0);
96
+
sn = mod(n.^2,N);
97
+
dn = (sn.*lambdad)/(N*2); % Profundidades [m]
98
+
% Otros
99
+
offsetdh = round(wide*(N/2)/dh); % Para centrar la posici�n del difusor
100
+
maxdh = round(max(dn)/dh); % Maxima profundidad
101
+
centdh = round(lx/2/dh); % Centro del mapa
102
+
103
+
%% MASCARA DEL DIFUSOR
104
+
mascdifusor = zeros(nx,ny); % Inicializaci�n de la mascara
105
+
for n = 1:N+1
106
+
% Lineas verticales
107
+
mascdifusor( (widedh*(n-1))+1+centdh-offsetdh+nPML,(1:maxdh)+nPML ) = 1;
108
+
% L�neas horizontales
109
+
if n<=N
110
+
switch tipo
111
+
case 'plano'
112
+
dnh = 0; % Profundidad de la cavidad n
113
+
case 'difusor'
114
+
dnh = round(dn(n)/dh); % Profundidad de la cavidad n
115
+
end
116
+
xhor = (((widedh*(n-1)):(widedh*(n)) )+1+centdh-offsetdh)+nPML;
117
+
yhor = maxdh-dnh+nPML;
118
+
mascdifusor( xhor,yhor ) = 1;
119
+
end
120
+
end
121
+
% L�nea horizontal en la parte inferior del difusor
122
+
mascdifusor( ((0:(widedh*N) )+1+centdh-offsetdh)+nPML,nPML ) = 1;
123
+
% Diferenciaci�n de la mascara
124
+
mascarax=1-(diff(mascdifusor).'~=0).';
125
+
mascaray=1-(diff(mascdifusor.')~=0).';
126
+
% Inversi�n de la mascara
127
+
mascdifusor = 1-mascdifusor;
128
+
129
+
%% REPRESENTACION MASCARAS Y PML
130
+
f = figure('Color',[1,1,1]);
131
+
subplot(1,3,1),pcolor(((1:nx)-nPML)*dh,((1:ny)-nPML)*dh,mascdifusor'), hold on
132
+
plot([0,0,nx-2*nPML,nx-2*nPML,0]*dh,[0,ny-2*nPML,ny-2*nPML,0,0]*dh,'r--','linewidth',2)
133
+
shading flat, axis equal,axis([-nPML,nx-nPML,-nPML,ny-nPML]*dh),title('M�scara')
134
+
subplot(1,3,2),pcolor(((1:nx-1)-nPML)*dh,((1:ny)-nPML)*dh,mascarax'), hold on
135
+
plot([0,0,nx-2*nPML,nx-2*nPML,0]*dh,[0,ny-2*nPML,ny-2*nPML,0,0]*dh,'r--','linewidth',2)
136
+
shading flat, axis equal,axis([-nPML,nx-nPML,-nPML,ny-nPML]*dh),title('M�scara en x')
137
+
subplot(1,3,3),pcolor(((1:nx)-nPML)*dh,((1:ny-1)-nPML)*dh,mascaray'), hold on
138
+
plot([0,0,nx-2*nPML,nx-2*nPML,0]*dh,[0,ny-2*nPML,ny-2*nPML,0,0]*dh,'r--','linewidth',2)
139
+
shading flat, axis equal,axis([-nPML,nx-nPML,-nPML,ny-nPML]*dh),title('M�scara en y')
140
+
fprintf('pulsa cualquier tecla para continuar');
141
+
pause
142
+
close(f)
143
+
144
+
%% Calculo
145
+
fig1 = figure('Color',[1,1,1]);
146
+
axis equal;axis([0,nx-nPML*2,0,ny-nPML*2]*dh)
147
+
xlabel('X [m]'),ylabel('Y [m]')
148
+
box on
149
+
colormap(flipud(bone(256)))
150
+
for tt=1:lenT
151
+
% Presion
152
+
px = px.*(1-gammax)-k*dt/dh*diff(ux).*mascdifusor;
153
+
py = py.*(1-gammay)-k*dt/dh*diff(uy')'.*mascdifusor;
154
+
155
+
% Excitaci�n
156
+
px(posNx,posNy) = w(tt);
157
+
py(posNx,posNy) = w(tt);
158
+
p = px+py;
159
+
160
+
% Velocidad
161
+
ux(2:nx,:) = ux(2:nx,:).*(1-gammaux(2:nx,:))-dt/rho/dh*diff(p).*mascarax;
162
+
uy(:,2:ny) = uy(:,2:ny).*(1-gammauy(:,2:ny))-dt/rho/dh*diff(p')'.*mascaray;
163
+
164
+
% Representaci�n gr�fica
165
+
if tt/10==round(tt/10)
166
+
map = abs(p);
167
+
surface(((1:nx)-nPML)*dh,((1:ny)-nPML)*dh,abs(mascdifusor-1)'*5)
168
+
hold on
169
+
pcolor(((1:nx)-nPML)*dh,((1:ny)-nPML)*dh,map'), hold off
170
+
set(gca,'Clim',[0 0.5]);shading interp
171
+
title(['Tiempo = ' num2str(round((tt)*1000*dt)) ' ms']);
172
+
drawnow
173
+
end
174
+
end