[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
Create FDTD_2D_Basico.m
Jose M. Requena Plens
(Sep 9, 2019, 4:12 AM +0200)
5bf69596
a06cc6ce
+153
+153
FDTD_2D_Basico.m
+153
FDTD_2D_Basico.m
···
1
+
%% =========================================================================
2
+
%
3
+
% Modelo simple de FDTD de un recinto en 2D con excitaci�n Ricker y PML
4
+
%
5
+
% Se incluyen 3 puntos de recepci�n y un ejemplo b�sico para introducir
6
+
% objetos rigidos que interactuan con el campo ac�stico.
7
+
%
8
+
% Jose Manuel Requena Plens (09/09/2019)
9
+
%
10
+
% =========================================================================
11
+
clc, clear, close all;
12
+
13
+
%% GENERAL
14
+
% Caracteristicas del medio
15
+
rho = 1.21; % Densidad del medio [kg/m3]
16
+
c = 341; % Velocidad de propagaci�n [m/s]
17
+
K = (c^2)*rho; % Modulo de compresibilidad
18
+
% Dimensiones del mapa en metros
19
+
lx = 5; % Ancho [m]
20
+
ly = 3; % Alto [m]
21
+
%
22
+
Origen = [1,2.5]; % Origen de la fuente
23
+
% Posicion de los receptores
24
+
rec1 = [4,0.25]; % Posicion del receptor 1 [m,m]
25
+
rec2 = [1,0.25];% Posicion del receptor 2 [m,m]
26
+
rec3 = [2,2]; % Posicion del receptor 3 [m,m]
27
+
% Tama�o/Ancho de la PML en metros
28
+
lPML = 0.25; % m
29
+
% Par�metros de calculo
30
+
dh = .01; % Definici�n espacial [m]
31
+
dt = dh/341/2; % Definici�n temporal [s]
32
+
ts = 15; % Tiempo de simulaci�n [ms]
33
+
34
+
%% MATRICES
35
+
% Dimensiones de la matriz del mapa
36
+
nPML = round(lPML/dh);
37
+
nx = round(lx/dh)+nPML*2; % Por 2 ya que la PML esta arriba y abajo
38
+
ny = round(ly/dh)+nPML*2; % Por 2 ya que la PML esta a izquierda y derecha
39
+
% Crear matriz del mapa
40
+
p = zeros(nx,ny); % Matrices de presion
41
+
px = zeros(nx,ny);
42
+
py = zeros(nx,ny);
43
+
ux = zeros(nx+1,ny); % Matriz de velocidad de particula en x
44
+
uy = zeros(nx,ny+1); % Matriz de velocidad de particula en y
45
+
46
+
%% POSICIONES DE LA FUENTE
47
+
%
48
+
posNy = round((Origen(2))/dh) + nPML;
49
+
posNx = round(Origen(1)/dh) + nPML;
50
+
51
+
52
+
%% EXCITACI�N
53
+
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
57
+
58
+
%% PML
59
+
% Gradiente de impedancias desde el valor del medio hasta un porcentaje de
60
+
% este
61
+
gammamax = 0.5; % Maxima reducci�n de la impedancia / porcentaje
62
+
% PML Izquierda y derecha
63
+
gammaux = zeros(nx+1,ny);
64
+
gammaux(1:nPML,:) = repmat(gammamax*((nPML:-1:1)'/nPML).^2,1,ny); % Izquierda
65
+
gammaux(1+end-nPML:end,:) = repmat(gammamax*((1:1:nPML)'/nPML).^2,1,ny); % Derecha
66
+
gammax = (gammaux(1:end-1,:)+gammaux(2:end,:))/2; % Conjunto
67
+
% PML Superior e inferior
68
+
gammauy = zeros(nx,ny+1);
69
+
gammauy(:,1:nPML) = repmat(gammamax*((nPML:-1:1)/nPML).^2,nx,1); % Inferior
70
+
gammauy(:,1+end-nPML:end) = repmat(gammamax*((1:1:nPML)/nPML).^2,nx,1); % Superior
71
+
gammay = (gammauy(:,1:end-1)+gammauy(:,2:end))/2; % Conjunto
72
+
73
+
%% Inicializacion de los vectores de los puntos de escucha
74
+
h = zeros(1,lenT);
75
+
h2 = zeros(1,lenT);
76
+
h3 = zeros(1,lenT);
77
+
78
+
%% Condiciones de contorno
79
+
ux(1,:) = -p(1,:)/rho/c;
80
+
ux(end,:) = p(end,:)/rho/c;
81
+
uy(:,1) = -p(:,1)/rho/c;
82
+
uy(:,end) = p(:,end)/rho/c;
83
+
84
+
%% Calculo
85
+
fig1 = figure('Color',[1,1,1]);
86
+
87
+
colormap(jet(256))
88
+
for tt=1:lenT
89
+
% Presion
90
+
px = px.*(1-gammax)-K*dt/dh*diff(ux);
91
+
py = py.*(1-gammay)-K*dt/dh*diff(uy')';
92
+
93
+
% Excitaci�n
94
+
px(posNx,posNy) = w(tt);
95
+
py(posNx,posNy) = w(tt);
96
+
p = px+py;
97
+
98
+
% Velocidad
99
+
ux(2:nx,:) = ux(2:nx,:).*(1-gammaux(2:nx,:))-dt/rho/dh*diff(p);
100
+
uy(:,2:ny) = uy(:,2:ny).*(1-gammauy(:,2:ny))-dt/rho/dh*diff(p')';
101
+
102
+
% Objetos (Se puede eliminar o modificar sin introducir valores fuera
103
+
% de las dimensiones del mapa)
104
+
% Barrera de 30cm de ancho
105
+
orX = 1; % Origen en X
106
+
orY = 0; % Origen en Y
107
+
wid = 0.3; % Ancho de la barrera
108
+
hei = 1.5; % Altura de la barrera
109
+
vecX = (round(orX/dh):round((orX+wid)/dh))+nPML;
110
+
vecY = (round(orY/dh)+1:round(hei/dh))+nPML;
111
+
ux(vecX,vecY) = 0;
112
+
uy(vecX,vecY) = 0;
113
+
114
+
% Respuesta al impulso
115
+
h(tt) = p(round(rec1(1)/dh)+nPML,round(rec1(2)/dh)+nPML);
116
+
h2(tt) = p(round(rec2(1)/dh)+nPML,round(rec2(2)/dh)+nPML);
117
+
h3(tt) = p(round(rec3(1)/dh)+nPML,round(rec3(2)/dh)+nPML);
118
+
119
+
% Representaci�n gr�fica
120
+
if tt/10==round(tt/10)
121
+
splmap = 10*log10(abs(p).^2/(2e-5)^2);
122
+
%splmap = splmap - max(splmap(:));
123
+
pcolor(((1:nx)-nPML)*dh,((1:ny)-nPML)*dh,splmap');
124
+
hold on;
125
+
plot([0,0,nx-2*nPML,nx-2*nPML,0]*dh,[0,ny-2*nPML,ny-2*nPML,0,0]*dh,'r--','linewidth',2)
126
+
plot([rec1(1),rec2(1),rec3(1)],[rec1(2),rec2(2),rec3(2)], 'ro', 'MarkerSize', 2,'MarkerFaceColor','r');
127
+
text([rec1(1),rec2(1),rec3(1)]+0.05,[rec1(2),rec2(2),rec3(2)],{'Rec1','Rec2','Rec3'},'Color','white')
128
+
hold off
129
+
axis equal;axis([-nPML,nx-nPML,-nPML,ny-nPML]*dh)
130
+
set(gca,'Clim',[50 110]);shading flat,cc = colorbar;
131
+
title(['Tiempo = ' num2str(round((tt)*1000*dt)) ' ms']);
132
+
xlabel('X [m]'),ylabel('Y [m]'),cc.Label.String = 'SPL (dB)';
133
+
drawnow
134
+
end
135
+
end
136
+
137
+
% Respuesta al impulso / Punto de escucha
138
+
time = seconds(t);
139
+
maxx = max([abs(h),abs(h2),abs(h3)]);
140
+
h = h./maxx; h2 = h2./maxx; h3 = h3./maxx;
141
+
T = array2timetable([h',h2',h3'],...% Datos iniciales
142
+
'RowTimes',time); % Columna de tiempo
143
+
fig= figure;
144
+
s = stackedplot(T,...
145
+
'Title','Se�al recibida',...
146
+
'DisplayLabels',{'Receptor 1', 'Receptor 2', 'Receptor 3'},...
147
+
'XLabel','Se�al recibida');
148
+
s.GridVisible = 1;
149
+
s.AxesProperties(1).YLimits = [-1.2,1.2];
150
+
s.AxesProperties(2).YLimits = [-1.2,1.2];
151
+
s.AxesProperties(3).YLimits = [-1.2,1.2];
152
+
fig.Children.FontName = 'Arial';
153
+
fig.Children.FontSize = 12;