Onda longitudinal plana. Grupo 22

Trabajo realizado por estudiantes
Título	Onda longitudinal plana. Grupo 22
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Irene Delgado Felpeto Ana Sanz García Lucía Reneses Doncel Francisco Javier Vela Cobos Marta Escaso Camacho
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región [math] \left[ 0, 4 \right] \times \left[\frac{-1}{2}, \frac{1}{2}\right] [/math].

En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math] 𝑇(𝑥,𝑦,𝑡) [/math], que depende de las dos variables espaciales [math] (𝑥,𝑦) [/math]
y del tiempo [math] 𝑡 [/math], y los deplazamientos.
De esta forma, si definimos [math] 𝑟_0(𝑥,𝑦) [/math] el vector de posición de los puntos de la placa en reposo, la posición de cada punto [math] (𝑥,𝑦) [/math] de la placa en un
instante de tiempo [math] 𝑡 [/math] viene dada por

[math] \vec{𝑟}(𝑥, 𝑦, 𝑡) = \vec{𝑟_0}(𝑥,𝑦) + \vec{𝑢}(𝑥,𝑦,𝑡) [/math].

Vamos a suponer que sobre la placa se ha aplicado una fuerza que ha provocado una vibración, de manera que los desplazamientos vienen
dados por la onda

[math] \vec{𝑢}(𝑥, 𝑦, 𝑡) = \vec{𝑎}cos(\vec{𝑏} ⋅ \vec{𝑟_0} − 𝑐𝑡) [/math],

donde [math] \vec{𝑎} [/math] se conoce como amplitud, [math] \vec{𝑏} [/math] es la fase que indica la dirección de propagación y [math] \frac{𝑐}{|\vec{𝑏}|} [/math] es la velocidad de propagación.
Si [math] \vec{𝑎} [/math] es paralelo a [math] \vec{𝑏} [/math] diremos que la onda es longitudinal mientras que si es perpendicular hablaremos de onda transversal.
En este trabajo vamos a centrarnos en las ondas longitudinales. Supondremos lo siguiente:

[math] \vec{𝑎} = \frac{\vec{𝑖}}{10}, \qquad \vec{𝑏} = \pi\vec{i}, \qquad t = 0 [/math].

En este caso, [math] \vec{𝑢}(𝑥,𝑦,𝑡) = \frac{cos(\pi x)}{10} \vec{i} [/math].

1 Mallado del sólido

La placa que estudiamos es un rectángulo: [math] \left[ 0, 4 \right] \times \left[\frac{-1}{2}, \frac{1}{2}\right] [/math].

Este rectángulo representa la región donde viven todas nuestras funciones: desplazamientos, temperatura, tensiones, etc.

Hacer un mallado significa tomar muchos puntos dentro del dominio, separados por una distancia fija h para poder hacer cálculos numéricos sobre ellos.

Nosotros elegimos: [math] h=0.1 [/math]

Así generamos valores de (x,y) entre [math] \left[ 0, 4 \right] \times \left[\frac{-1}{2}, \frac{1}{2}\right] [/math].

Utilizamos: [math] \qquad x = -0.5:0.1:0.5; \qquad y = 0:0.1:4; \qquad [X,Y] = meshgrid(x,y); [/math]

Meshgrid crea dos matrices [math] X [/math] e [math] Y [/math]. Cada pareja [math] (X(i,j),Y(i,j)) [/math] es una coordenada del mallado.

Hemos representado el mallado de dos formas:

[math](1)[/math] Como nube de puntos, mostrando únicamente los puntos interiores donde se harán los cálculos.

[math](2)[/math] Como malla de líneas, dibujando las curvas [math] 𝑥=cte [/math] e [math] y=cte [/math], lo que permite visualizar la estructura completa del mallado.

1.1 Código y representación

h = 0.1;
y = -0.5:h:0.5;
x = 0:h:4;

[X, Y] = meshgrid(x, y);

%  Figura con puntos interiores
figure("Name","(1) Mallado - Puntos interiores");
plot(X(:), Y(:), ".k", "MarkerSize", 10);
axis equal tight;
grid on;
title("Mallado de puntos interiores (h = 0.1)");
xlabel("x"); ylabel("y");

%  Figura con líneas estilo 'curvas parametrizadas'
figure("Name","(1b) Mallado - Malla tipo curvas");
hold on;

% Líneas horizontales y verticales como parametrizaciones
for j = 1:size(X,1)
    plot(X(j,:), Y(j,:), "-k");
end
for i = 1:size(X,2)
    plot(X(:,i), Y(:,i), "-k");
end

axis equal tight;
grid on;
title("Mallado representado como familia de curvas");
xlabel("x"); ylabel("y");
hold off;

2 Temperatura

La temperatura no se da en coordenadas [math] X,Y [/math] sino en coordenadas polares:

[math]\rho[/math] = distancia al origen

[math]\theta[/math] = ángulo respecto al eje x

La fórmula [math] T(\rho ,\theta )=e^{-\theta } [/math] dice que la temperatura solo depende del ángulo, no de la distancia.

Por ello convertimos cada punto del mallado a coordenadas polares: [math] [TH, R] = cart2pol(X, Y) [/math]

[math]TH[/math] contiene todos los ángulos [math]\theta [/math]

[math]R[/math] contiene todos los radios [math]\rho [/math]

Evaluamos la temperatura en cada punto con [math] T = exp(-TH) [/math]

Hemos representado la temperatura de tres formas:

[math]1.[/math] Mapa de color (temperatura en 2D)

Representamos la temperatura con un código de colores. Cada punto del dominio tiene un color según su temperatura [math] T(\rho ,\theta )=e^{-\theta } [/math]. Permite ver rápidamente qué zonas son más calientes o más frías.

[math]2.[/math] Superficie 3D de la temperatura

Aquí representamos la temperatura como una “altura”. Esto muestra de forma intuitiva cómo cambia la temperatura según el ángulo.

[math]3.[/math] Curvas de nivel de la temperatura

Dibujamos líneas donde la temperatura es constante. Estas curvas permiten ver cómo varía 𝑇 en la placa y reflejan su dependencia angular

2.1 Código y representación

% Conversión a polares (como recomienda el PDF)
[TH, R] = cart2pol(X, Y);

T = exp(-TH);

%  Mapa 2D tipo 'heatmap' 
figure("Name","(2a) Temperatura - mapa 2D");
pcolor(X, Y, T);
shading interp; colorbar;
axis equal tight;
title("Temperatura T(\\rho,\\theta)=e^{-\\theta}");
xlabel("x"); ylabel("y");

%  Superficie 3D suave 
figure("Name","(2b) Temperatura - superficie 3D");
surf(X, Y, T, "EdgeColor","none");
shading interp; colorbar;
view(45,30);
axis tight;
title("Superficie de T(\\rho,\\theta)");
xlabel("x"); ylabel("y"); zlabel("T");

% Curvas de nivel (estilo 'curvas planas' del PDF)
figure("Name","(2c) Temperatura - curvas de nivel");
contourf(X, Y, T, 30);
colorbar;
axis equal tight;
title("Curvas de nivel de T");
xlabel("x"); ylabel("y");


.

3 Gradiente

El campo gradiente de la función escalar [math]T(x,y)[/math] es el campo vectorial formado por las derivadas parciales de [math]T[/math]:

[math]\Delta T=\left ( \frac{\partial T}{\partial x},\frac{\partial T}{\partial y}\right )[/math]

En este caso, las curvas de nivel [math]T[/math] son líneas de [math]\theta [/math] constante, con origen en el origen de coordenadas (dirección radial) y el gradiente apunta en diercción angular.

3.1 Cálculo del gradiente en coordenadas polares

[math]\Delta T=\frac{\partial T}{\partial \rho }\hat{e_{\rho }}+\frac{1}{\rho }\frac{\partial T}{\partial \theta }\hat{e_{\rho }}[/math]

-Derivada respecto de [math]\theta[/math]: Ya que [math]T[/math] no depende de [math]\theta[/math], la derivada parcial respecto [math]\theta[/math] es nula [math]\frac{\partial T}{\partial \rho }=0 [/math]

-Derivada respecto de [math]\rho[/math]: [math]\frac{\partial T}{\partial \theta}=-e^{-\theta }[/math]

El gradiente queda de la forma: [math] \Delta T=-\frac{e^{-\theta }}{\rho }\hat{e_{\theta }} [/math]

Se observa que únicamente cuenta con componente en dirección de [math]\hat{e_{\theta }} [/math]

3.2 Conclusión

Se deduce de la expresión que define la temperatura que, para una temperatura concreta, las curvas de nivel resultan líneas a [math]\rho[/math] constante. Es decir, semirrectas que parten del origen (dirección radial)

[math]e^{-\theta }=T\rightarrow \theta =-lnT[/math]

Por otro lado, el gradiente apunta en dirección angular [math] \Delta T=-\frac{e^{-\theta }}{\rho }\hat{e_{\theta }} [/math]

Por lo que ambos campos son perpendiculares entre sí en cada punto, como se puede observar en la gráfica obtenida con octave.

3.3 Código y representación

% Mallado polar para T(rho,theta)
h = 0.1
rho = 0.1:h:2        % evitando radio nulo
theta = 0:h:2*pi;    
[R,Th] = meshgrid(rho, theta);   % mallado polar

% Temperatura T = exp(-theta)
T=exp(-Th);

% Gradiente en coordenadas polares
dpr = 0.*R;                       % dT/dr=0
dpth = -exp(-Th);                   % dT/dth =-e^{-th}

% Convertir gradiente a coordenadas cartesianas
Ex = dpr.*cos(Th)-(1./R).*dpth.*sin(Th); 
Ey = dpr.*sin(Th)+(1./R).*dpth.*cos(Th);

% Convertir puntos a coordenadas cartesianas
X = R.*cos(Th);
Y = R.*sin(Th);

% Graficar campo vectorial
figure;
quiver(X, Y, Ex, Ey);
hold on;
title('Campo vectorial ∇T');
xlabel('x');
ylabel('y');
axis equal;

% Dibujar curvas de nivel T
contour(X, Y, T, 10,'LineWidth',1.2);
title('Gradiente y curvas de nivel de T');

4 Campo vectorial [math]u[/math]

El campo vectorial [math] \bar{u}(x,y)=\frac{cos(\pi x)}{10}\hat{i} [/math] es un campo horizontal (sólo tiene componente en el eje x) ya que se trata de una onda longitudinal.

La función [math]cos(\pi x)[/math] es periódica con periodo [math]2[/math], por lo que la dirección del flujo cambia (derecha, cera, izquierda, cero, derecha...) a lo largo de [math]x[/math] repitiéndose la magnitud cada [math]2[/math] unidades.

[math]x=0 \rightarrow cos(\pi x)=1 \qquad [/math] (Flecha máxima a la derecha)

[math]x=n*\frac{1}{2} \qquad n=1,2,3... \rightarrow cos(\pi x)=0 \qquad [/math] (Flecha nula)

[math]x=n*1 \qquad n=1,2,3... \rightarrow cos(\pi x)=-1 \qquad [/math] (Flecha máxima a la izquierda)

4.1 Código y representación

% Definición campo vectorial
Ux=cos(pi*X)/10;    % componente x
Uy=zeros(size(Ux));    % componente y

% Dibujo campo vectorial
figure; 
hold on;
plot(X(:), Y(:), ".b", "MarkerSize", 5);
quiver(X, Y, Ux, Uy,"k","Markersize",8);
title('Campo vectorial u');
xlabel('x'); ylabel('y'); axis equal;
grid on;
hold off;

5 Deformación de la placa inducida por una onda longitudinal

Procedemos al estudio de la deformación de la placa rectangular, que ocupa la región [math] \left[ 0, 4 \right] \times \left[\frac{-1}{2}, \frac{1}{2}\right] [/math], cuando sobre ella actúa una onda longitudinal.

En este caso el desplazamiento de la onda viene dado por: [math] \vec{u} (x,y)=\frac{cos(πx)}{10}\vec{i} [/math].

El desplazamiento es paralelo a la dirección de propagación (eje x).

Al no haber desplazamiento vertical, la placa solo se estira o comprime horizontalmente. Los puntos de la placa no se desplazan todos igual, podemos estudiarlo de la siguiente manera:

· Observamos la deformación máxima en los puntos donde cos(πx) es máximo, como x=0, x=2 y x=4

· La deformación es nula en los puntos donde cos(πx)=0, como x=1 y x=2

· La placa se comprime cuando la onda es negativa y se dilata cuando es positiva

5.1 Código y representación

Valor del desplazamiento Ux en cada punto (x,y) y malla física deformada

%% Dominio 
Nx = 160; Ny = 40;
x = linspace(0,4,Nx);
y = linspace(-0.5,0.5,Ny);
[X, Y] = meshgrid(x,y);

%% Desplazamiento longitudinal (onda)
ux = (1/10) * cos(pi * X);   % ux(x,y) (no depende de y)
uy = zeros(size(X));

%% Malla deformada (solo en x)
Xdef = X + ux;
Ydef = Y;

%% Plot: antes (campo sobre malla original) y despues (malla deformada coloreada) 
figure('Units','normalized','Position',[0.05 0.15 0.9 0.6]);

%% Izquierda: campo ux sobre la malla ORIGINAL 
subplot(1,2,1);
% coloreado de ux sobre la malla original
surf(X, Y, zeros(size(X)), ux, 'EdgeColor', 'none');
view(2);
axis equal tight;
colormap(subplot(1,2,1), parula);
colorbar;
xlabel('x'); ylabel('y');
title('ANTES: desplazamiento u_x sobre malla original');
hold on;
% añadir flechas muestreadas para ver dirección/amplitud (no saturan)
step = max(1, round(Nx/20));
quiver(X(1:step:end,1:step:end), Y(1:step:end,1:step:end), ...
       ux(1:step:end,1:step:end), uy(1:step:end,1:step:end), ...
       0.6, 'k', 'LineWidth', 0.8);
hold off;

%% Derecha: malla DEFORMADA coloreada por el mismo ux
subplot(1,2,2);
% superficie coloreada sobre la malla deformada
surf(Xdef, Ydef, zeros(size(Xdef)), ux, 'EdgeColor', 'none');
view(2);
axis equal tight;
colormap(subplot(1,2,2), parula);
colorbar;
xlabel('x'); ylabel('y');
title('DESPUÉS: malla deformada coloreada por u_x');

hold on;
% superponer contorno de la malla original para comparar (líneas finas negras)
for i=1:Ny
    plot(X(i,:), Y(i,:), 'k-', 'LineWidth', 0.4);   % horizontales originales
end
for j=1:Nx
    plot(X(:,j), Y(:,j), 'k-', 'LineWidth', 0.4);   % verticales originales
end
% y superponer contorno de la malla deformada (líneas rojas más gruesas)
for i=1:Ny
    plot(Xdef(i,:), Ydef(i,:), 'r-', 'LineWidth', 1.0);
end
for j=1:Nx
    plot(Xdef(:,j), Ydef(:,j), 'r-', 'LineWidth', 1.0);
end
hold off;

6 Divergencia del campo de desplazamiento

Partimos del campo de desplazamiento [math] \vec{u}(x, y)= \frac{cos(πx)}{10} \vec{i} [/math]

La divergencia en cartesianas la calculamos de la siguiente manera: [math] \triangledown \circ \vec{u} = \frac{\partial u_{x}}{\partial x} + \frac{\partial u_{y}}{\partial y} [/math] . En este caso en específico Uy=0 por lo que [math] \triangledown \circ \vec{u} = \frac{\partial u_{x}}{\partial x} [/math].

Finalmente nos queda que la divergencia del campo de desplazamiento en todos los puntos del sólido es:

[math] \triangledown \circ \vec{u} = \frac{\partial }{\partial x} (\frac{cos(πx)}{10}) = \frac{-π }{10}sin(πx) [/math]

6.1 Código y representación

En las gráficas proporcionadas a continuación observamos la superficie de la divergencia tanto en 3D como en 2D. Sin nos fijamos en la primera imagen, podemos identificar los máximos y mínimos gráficamente, ubicando los máximos en x=1.5 y x=3.5, y los mínimos en x=0.5 y x=2.5.

Superficies de divergencia 2D y 3D

%% Malla
Nx = 200;
Ny = 60;

x = linspace(0, 4, Nx);
y = linspace(-0.5, 0.5, Ny);
[X, Y] = meshgrid(x, y);

%%  Divergencia 
div_u = - (pi/10) * sin(pi * X);

%%  Superficie limpia
figure;

surf(X, Y, div_u, 'EdgeColor', 'none');
colormap turbo;      % (puedes cambiarlo si quieres)
colorbar;

xlabel('x');
ylabel('y');
zlabel('\nabla \cdot u');
title('Superficie de la divergencia \nabla \cdot u(x,y)');

view(45, 30);  % ángulo 3D

En la siguiente gráfica observamos de otra manera donde se encuentran los máximos, marcados con una linea vertical blanca. Para encontrar estos puntos analíticamente, partimos de la base de que sabemos que los puntos de mayor divergencia se encontraran cuando -sin(πx)=1 con [math] x \epsilon [0, 4] [/math], esto se da en x=1.5 y x=3.5.

Máximos de la divergencia

%% Malla
Nx = 200;
Ny = 60;

x = linspace(0, 4, Nx);
y = linspace(-0.5, 0.5, Ny);
[X, Y] = meshgrid(x, y);

%%  Divergencia
div_u = - (pi/10) * sin(pi * X);

%%  Posiciones teóricas de los máximos
x_max = [1.5, 3.5];  % donde sin(pi x) = -1

%% Encontrar los puntos de malla más cercanos a los máximos
[~, idx1] = min(abs(x - x_max(1)));
[~, idx2] = min(abs(x - x_max(2)));

Xmax_points = [X(:, idx1); X(:, idx2)];
Ymax_points = [Y(:, idx1); Y(:, idx2)];

%% Dibujo 
figure;

surf(X, Y, zeros(size(X)), div_u, 'EdgeColor', 'none');
view(2);
axis equal tight;
colorbar;

title('Divergencia \nabla \cdot u con máximos marcados');
xlabel('x');
ylabel('y');
hold on;

% Marcar líneas verticales de máximos
plot3(X(:, idx1), Y(:, idx1), ones(Ny,1), ...
      'wo', 'MarkerSize', 8, 'MarkerFaceColor', 'w');
plot3(X(:, idx2), Y(:, idx2), ones(Ny,1), ...
      'wo', 'MarkerSize', 8, 'MarkerFaceColor', 'w');

7 Rotacional del campo de desplazamiento

Siendo [math] \vec{u}(x, y)= \frac{cos(πx)}{10} \vec{i} [/math] el campo de desplazamiento con el que estamos trabajando y partiendo de la fórmula del rotaciones en cartesianas [math] \triangledown [/math] x [math] \vec{u} = (\frac{\partial u_{3}}{\partial x_{2}}-\frac{\partial u_{2}}{\partial x_{3}})\vec{i} - (\frac{\partial u_{3}}{\partial x_{1}}-\frac{\partial u_{1}}{\partial x_{3}})\vec{j} [/math] , procedemos a analizar y sustituir para realizar el estudio del rotación del campo de desplazamiento.

En este caso no contamos con componente [math] \vec{j} [/math], por los que el rotacional se nos queda como: [math] \triangledown [/math] x [math] \vec{u} = (\frac{\partial u_{3}}{\partial x_{2}}-\frac{\partial u_{2}}{\partial x_{3}})\vec{i} [/math]

Sustituimos en la ecuación y analizamos la solución:

[math] \triangledown [/math] x [math] \vec{u} = \frac{\partial 0}{\partial x}-0 [/math]

[math] \triangledown [/math] x [math] \vec{u} = \vec{0} [/math]

Al ser el resultado nulo sabemos que el campo en todo el sólido es nulo, por no que ho hay giro local y denominamos el campo como irrotacional.

7.1 Código y representación

Rotacional

%% Apartado 7: rotacional de u en Omega = [0,4] x [-0.5,0.5]
Nx = 300; Ny = 120;
x = linspace(0, 4, Nx);
y = linspace(-0.5, 0.5, Ny);
[X, Y] = meshgrid(x, y);

% Campo u
ux = (1/10) * cos(pi * X);   % ux(x,y)
uy = zeros(size(ux));        % uy = 0

% Cálculo numérico del rotacional (componente z en 2D)
% gradient devuelve derivadas: [dF/dx, dF/dy] si se pasan vectores x,y
[dux_dx, dux_dy] = gradient(ux, x, y);
[duy_dx, duy_dy] = gradient(uy, x, y);

curl_z = duy_dx - dux_dy;    % should be all zeros
curl_abs = abs(curl_z);

% Estadísticas numéricas
max_abs = max(curl_abs(:));
fprintf('Máx abs(curl_z) en la malla = %.3e\n', max_abs);


%% Dibujo 2: superficie 3D del curl (debe ser plano en 0) 
figure('Name','curl_z (superficie)','NumberTitle','off');
surf(X, Y, curl_z, 'EdgeColor', 'none');
colormap parula; colorbar;
xlabel('x'); ylabel('y'); zlabel('(\nabla\times u)_z');
title('Componente z del rotacional (numérico)');
view([45 30]); axis tight; grid on;

8 Modelado de desplazamientos y tensiones en la placa vibrante.

Primero, se define el campo de desplazamiento sabiendo que este se modela mediante: [math] \vec{𝑢}(𝑥, 𝑦, 𝑡) = \vec{𝑎}cos(\vec{𝑏} ⋅ \vec{𝑟_0} − 𝑐𝑡) [/math].
Teniendo en cuenta las condiciones iniciales [math] (t = 0) [/math] dadas: [math] \vec{𝑎} = \frac{\vec{𝑖}}{10}, \vec{𝑏} = \pi\vec{i} [/math], este resulta [math] \vec{𝑢}(𝑥,𝑦,𝑡) = \frac{cos(\pi x)}{10} \vec{i} [/math].

Para caracterizar localmente la deformación de la placa, se define el tensor de deformaciones:

[math] ε(\vec{𝑢})= \frac{1}{2}(∇\vec{𝑢}+(∇\vec{𝑢})^T) [/math].

Este tensor, es la parte simétrica de [math] \vec{𝑢} [/math], y mide elongaciones, compresiones, etc. en cada punto de la placa.

Por otro lado, se cuenta con el tensor de tensiones, el cual sirve para materiales elásticos, homogéneos e isotrópicos. Este se relaciona con el anterior mediante la siguiente expresión:

[math] σ = λ(∇·\vec{𝑢})I + 2με [/math],

donde [math] λ [/math] y [math] μ [/math] son los conocidos como coeficientes de Lamé, los cuales dependen de las propiedades elásticas de cada material; y donde [math] I [/math], es el tensor identidad de [math] R^3 [/math].
Tomando [math] λ = μ = 1 [/math], este tensor se simplifica y queda:

[math] σ = (∇·\vec{𝑢})I + ∇\vec{𝑢}+(∇\vec{𝑢})^T [/math]

8.1 Cálculo de tensiones normales

Para estudiar el efecto de la deformación en cada dirección, se buscan las tensiones normales a los ejes:

Para las del eje [math]x[/math] se elige la posición en la que coinciden la primera fila y la primera columna de la matriz del tensor de tensiones:

[math] σ_{xx} = \vec{i}·σ·\vec{i} \rightarrow (∇·\vec{𝑢})\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} + \begin{pmatrix} \frac{\partial u_{x}}{\partial x} & \frac{\partial u_{x}}{\partial y} \\ \frac{\partial u_{y}}{\partial x} & \frac{\partial u_{y}}{\partial y} \end{pmatrix} + \begin{pmatrix} \frac{\partial u_{x}}{\partial x} & \frac{\partial u_{y}}{\partial x} \\ \frac{\partial u_{x}}{\partial y} & \frac{\partial u_{y}}{\partial y} \end{pmatrix} = \begin{pmatrix} (∇·\vec{𝑢})+\frac{\partial u_{x}}{\partial x}+\frac{\partial u_{x}}{\partial x} & \frac{\partial u_{x}}{\partial y} + \frac{\partial u_{y}}{\partial x} \\ \frac{\partial u_{y}}{\partial x} + \frac{\partial u_{x}}{\partial y} & (∇·\vec{𝑢})+ \frac{\partial u_{y}}{\partial y} + \frac{\partial u_{y}}{\partial y} \end{pmatrix} \rightarrow (∇·\vec{𝑢}) + 2\frac{\partial u_{x}}{\partial x} [/math]

Y para el caso del eje [math]y[/math] se elige la posición en la que coinciden la segunda fila y la segunda columna: [math] σ_{yy} = \vec{j}·σ·\vec{j} = (∇·\vec{𝑢}) + 2\frac{\partial u_{y}}{\partial y} [/math]

En este caso, se tiene: [math] \qquad u_{x} = \frac{1}{10}cos(\pi x), \qquad u_{y} = 0, \qquad [/math] por lo que:

[math] \frac{\partial u_{x}}{\partial x} = -\frac{\pi}{10}sin(\pi x), \qquad \frac{\partial u_{y}}{\partial y} = 0 [/math]

[math] ∇·\vec{𝑢} = \frac{\partial u_{x}}{\partial x} + \frac{\partial u_{y}}{\partial y} = -\frac{\pi}{10}sin(\pi x) [/math]

Así, las tensiones normales quedan:

[math] σ_{xx} = (∇·\vec{𝑢}) + 2\frac{\partial u_{x}}{\partial x} = -\frac{\pi}{10}sin(\pi x) + 2(-\frac{\pi}{10}sin(\pi x)) = -\frac{3\pi}{10}sin(\pi x) [/math]

[math] σ_{yy} = (∇·\vec{𝑢}) + 2\frac{\partial u_{y}}{\partial y} = -\frac{\pi}{10}sin(\pi x) + 0 = -\frac{\pi}{10}sin(\pi x) [/math]

8.2 Código y representación

La gráfica de la izquierda muestra cómo la onda longitudinal genera tensiones internas en la dirección de propagación. Los valores positivos se representan con colores cálidos que indican tracción (es decir, los lugares en los que el material se estira) y los valores negativos, con colores fríos que indican compresión. En el caso de la gráfica de la derecha, esta aparece casi sin variaciones, ya que el desplazamiento no depende de la variable [math]y[/math], solo de la [math]x[/math]. Esto indica que en la dirección vertical no hay apenas tensiones normales inducidas por la onda longitudinal.

Representación de las tensiones normales en las direcciones de i y j en la placa

% Se define la malla de la placa
Nx = 400; Ny = 100;   % Ajustamos la resolución para mantener una buena relación x/y
x = linspace(0, 4, Nx);
y = linspace(-0.5, 0.5, Ny);
[X,Y] = meshgrid(x,y);

% Desplazamiento longitudinal:
% u(x,y) = ( (1/10) * cos(pi*x), 0 )
Ux = (1/10) * cos(pi*X);
Uy = zeros(size(X));

% Gradientes de u
[dUx_dx, dUx_dy] = gradient(Ux, x, y);
[dUy_dx, dUy_dy] = gradient(Uy, x, y);

% Tensor de deformaciones epsilon = 1/2 (grad u + grad u^T)
eps_xx = dUx_dx;
eps_xy = 0.5*(dUx_dy + dUy_dx);
eps_yy = dUy_dy;

% Divergencia de u
div_u = dUx_dx + dUy_dy;

% Coeficientes de Lamé
lambda = 1; mu = 1;

% Tensor de tensiones sigma = lambda*(div u) I + 2 mu eps
sigma_xx = lambda*div_u + 2*mu*eps_xx;
sigma_yy = lambda*div_u + 2*mu*eps_yy;
sigma_xy = 2*mu*eps_xy;

% Tensiones normales en dirección i y j
sigma_i = sigma_xx;   % Dirección x
sigma_j = sigma_yy;   % Dirección y

% Gráficas
figure;
set(gcf, 'Position', [100, 100, 1000, 300]); 
subplot(1,2,1);
imagesc(x,y,sigma_i);
axis xy;
axis equal tight; 
colorbar;
title('Tensión normal en dirección x (σ_{xx})');
xlabel('x'); ylabel('y');

subplot(1,2,2);
imagesc(x,y,sigma_j);
axis xy;
axis equal tight; 
colorbar;
title('Tensión normal en dirección y (σ_{yy})');
xlabel('x'); ylabel('y');

9 Tensiones tangenciales respecto al plano ortogonal a [math] \vec{i} [/math]

Para estudiar las tensiones tangenciales respecto al plano ortogonal a [math] \vec{i} [/math] , es decir: [math] \left | \sigma·\vec{i} - (\vec{i}·\sigma·\vec{i}) \vec{i} \right | [/math]

Necesitamos los valores de la divergencia (apartado 6) y tensor de tensiones [math] \sigma [/math] con [math] \lambda =\mu = 1 [/math] (apartado 8):

Divergencia: [math] \bigtriangledown · u = - \frac{\pi }{10} sin (\pi x) [/math]

Tensor de tensiones con [math] \lambda =\mu = 1 [/math] [math] \qquad \qquad \sigma = \lambda (\bigtriangledown · u)\mathbb{I} + 2\mu \varepsilon \qquad \Rightarrow \qquad \sigma (x,y)=\begin{pmatrix} - 3\frac{\pi }{10} sin (\pi x) & 0 \\ 0 & - \frac{\pi }{10} sin (\pi x) \\ \end{pmatrix} [/math]

[math] \sigma [/math] es diagonal, ya que los valores de [math] σ_{12},σ_{21} [/math] son nulos.

Ahora sí podemos calcular componente a componente, al ser el plano ortogonal a [math] \vec{i} [/math], tenemos que evaluar [math] σ_{11},σ_{21} [/math]:

[math] \sigma·\vec{i} = \begin{pmatrix} σ_{11}\\σ_{21} \end{pmatrix} = \begin{pmatrix} - 3\frac{\pi }{10} sin (\pi x) \\ 0 \end{pmatrix} [/math]

[math] (\vec{i}·\sigma·\vec{i})= σ_{11} = - 3\frac{\pi }{10} sin (\pi x) [/math]

Por lo tanto, las tensiones tangenciales respecto al plano ortogonal a [math] \vec{i} [/math] , [math] \left | \sigma·\vec{i} - (\vec{i}·\sigma·\vec{i}) \vec{i} \right | = \begin{pmatrix} 0 \\ 0 \end{pmatrix} [/math]

Es decir, la componente tangencial, sobre las caras normales al eje [math] x [/math] es 0 en todo el dominio. Por lo tanto, las tracciones sobre esas caras son puramente normales, no hay esfuerzo cortante.

No hay puntos donde sean mayores las tensiones tangenciales, porque son nulas.

9.1 Comparación con puntos de mayor deformación de la malla

La deformación relevante para este apartado es [math] \varepsilon _{11} = - \frac{\pi }{10} sin (\pi x) [/math]

[math] \left | \varepsilon _{11} \right | [/math] es máxima cuando [math] \left | sin (\pi x) \right | = 1 [/math] esto ocurre para [math] \pi x = \frac{\pi }{2} + k \pi \Rightarrow x = \frac{1}{2} + k, (k\epsilon \mathbb{Z}) [/math]

En el intervalo [math] x\epsilon [0,4] [/math] los valores son: [math] x = 0.5, 1.5, 2.5, 3.5 [/math]

Estos vectores verticales (en la dirección [math] y [/math]) tendrá la deformación máxima en el valor absoluto. Las máximas deformaciones ([math] x = 0.5, 1.5, 2.5, 3.5 [/math]) no implican que exista tensiones tangenciales. Como hemos demostrado en este apartado, las tensiones tangenciales son [math] 0 [/math].

Lo que sí varía y tiene máximos en esos puntos [math] x [/math] son las tensiones normales [math] σ_{11},σ_{21} [/math] y sus valores absolutos alcanzan su máximo donde [math] \left | sin (\pi x) \right | = 1 [/math]. Concretamente:

[math] σ_{11} = - 3\frac{\pi }{10} sin (\pi x) [/math] , [math] σ_{22} = - \frac{\pi }{10} sin (\pi x) [/math]

Conclusión, la onda longitudinal produce cambios de volumen en la dirección de propagación [math] x [/math], no genera esfuerzos cortantes. Por eso las tensiones tangenciales son nulas incluso donde la deformación absoluta es máxima.

9.2 Código y representación

Simplificamos el código mostrando el apartado 9. Para realizar este código matlab consideramos los apartados anteriores en los que ya están detallados los: parámetros del problema, mallado, campo de desplazamientos, derivadas del desplazamiento, tensor deformación, divergencia y tensiones.

% Tensiones tangenciales respecto al plano normal a i (apartado 9)

%% Tensión tangencial respecto a i (normal = (1,0))
% t = sigma*i = (sigma11, sigma12)
% t_tan = t - (i*sigma*i)*i = (0,0)
t_tan_i_x = zeros(size(XX));
t_tan_i_y = zeros(size(YY));
mag_tan_i = zeros(size(XX));

%% Figuras del apartado 9
figure('Name','Apartado 9: Tensiones tangenciales respecto a i');

% σ11
subplot(2,2,1);
imagesc(x,y,sigma11); axis xy; colorbar;
title('\sigma_{11}'); xlabel('x'); ylabel('y'); axis equal;

% eps11
subplot(2,2,2);
imagesc(x,y,abs(eps11)); axis xy; colorbar;
title(' \epsilon_{11} '); xlabel('x'); ylabel('y'); axis equal;

% Magnitud tensión tangencial
subplot(2,2,3);
imagesc(x,y,mag_tan_i); axis xy; colorbar;
title('Magnitud tensión tangencial wrt i'); xlabel('x'); ylabel('y'); axis equal;

% Quiver tensión tangencial
subplot(2,2,4);
plot(XX(:),YY(:),'.','Color',[0.7 0.7 0.7]); hold on;
quiver(XX,YY,t_tan_i_x,t_tan_i_y,0,'b');
title('Quiver tensión tangencial wrt i'); xlabel('x'); ylabel('y'); axis equal;

fprintf("Máxima tensión tangencial (i): %.6e\n", max(mag_tan_i(:)));

10 Tensiones tangenciales respecto al plano ortogonal a [math] \vec{j} [/math]

Para estudiar las tensiones tangenciales respecto al plano ortogonal a [math] \vec{j} [/math] , es decir: [math] \left | \sigma·\vec{j} - (\vec{j}·\sigma·\vec{j}) \vec{j} \right | [/math]

Utilizamos el mismo procedimiento y calculamos componente a componente, considerando de la misma forma [math] \sigma (x,y) [/math] del apartado (8) y evaluando en [math] σ_{12},σ_{22} [/math]:

[math] \sigma·\vec{j} = \begin{pmatrix} σ_{12}\\σ_{22} \end{pmatrix} = \begin{pmatrix} 0 \\ - \frac{\pi }{10} sin (\pi x) \end{pmatrix} [/math]

[math] (\vec{j}·\sigma·\vec{j}) = σ_{22} = - \frac{\pi }{10} sin (\pi x) [/math]

Por lo tanto, las tensiones tangenciales respecto al plano ortogonal a [math] \vec{j} [/math], [math] \left | \sigma·\vec{j} - (\vec{j}·\sigma·\vec{j}) \vec{j} \right | = \begin{pmatrix} 0 \\ 0 \end{pmatrix} [/math]

En este caso, la traccion tangencial sobre planos normales al eje [math] y [/math] también es nula en todo el dominio. Por lo que podemos decir que [math] \sigma [/math] es diagonal en este problema.

De forma similar al apartado anterior, no hay puntos donde sean mayores las tensiones tangenciales, porque son nulas.

Esto es completamente lógico para una onda longitudinal cuyo desplazamiento va exclusivamente en dirección x y que produce sólo componentes normales de tensión en esta formulación lineal simple.

10.1 Comparación con puntos de mayor deformación de la malla

De la misma forma que en el apartado 9.1, al ser matriz diagonal, no hay puntos con tensiones tangenciales mayores

10.2 Código y representación

Simplificamos el código mostrando el apartado 10, para realizar este código matlab consideramos los apartados anteriores en los que ya están detallados los: parámetros del problema, mallado, campo de desplazamientos, derivadas del desplazamiento, tensor deformación, divergencia y tensiones.

%% Tensión tangencial respecto a j
% t = sigma*j = (sigma12, sigma22)
% t_tan = (sigma12,0)  → pero sigma12 = 0 siempre
t_tan_j_x = sigma12;
t_tan_j_y = zeros(size(YY));
mag_tan_j = abs(t_tan_j_x);

%% Gráficas del apartado 10
figure('Name','Apartado 10: Tensiones tangenciales respecto a j');

% σ22
subplot(2,2,1);
imagesc(x,y,sigma22); axis xy; colorbar;
title('\sigma_{22}'); xlabel('x'); ylabel('y'); axis equal;

% eps11
subplot(2,2,2);
imagesc(x,y,abs(eps11)); axis xy; colorbar;
title('\epsilon_{11}'); xlabel('x'); ylabel('y'); axis equal;

% Magnitud tensión tangencial
subplot(2,2,3);
imagesc(x,y,mag_tan_j); axis xy; colorbar;
title('Magnitud tensión tangencial wrt j'); xlabel('x'); ylabel('y'); axis equal;

% Quiver tangencial
subplot(2,2,4);
plot(XX(:),YY(:),'.','Color',[0.7 0.7 0.7]); hold on;
quiver(XX,YY,t_tan_j_x,t_tan_j_y,0,'r');
title('Quiver tensión tangencial wrt j'); xlabel('x'); ylabel('y'); axis equal;

fprintf("Máxima tensión tangencial (j): %.6e\n", max(mag_tan_j(:)));

11 Masa de la placa

Si la densidad de la placa viene dada en coordenadas polares por [math] \delta (\rho ,\theta )= 1 + e^{\rho ^2 cos \theta} [/math]

Queremos calcular la masa total de la placa dada por el dominio [math] D = \left\{ (x,y); x\epsilon [0,4], y\epsilon [-\frac{1}{2},\frac{1}{2}] \right\} [/math]

La masa de la integral dobles es: [math] M = \int\int_{D} d(x,y)dA [/math]

Donde en coordenadas cartesianas, para poder integrar en el dominio del rectangulo, [math] \rho = {\sqrt{x^2+y^2}} [/math] y [math] \theta = atan2(x,y) \rightarrow cos \theta =\frac{x}{\rho} [/math]

[math] \delta (\rho ,\theta )= 1 + e^{\rho ^2 cos \theta} [/math]

[math] \delta (x,y)= 1 + e^{(x^2 + y^2)·\frac{x }{\sqrt{x^2+y^2}}} [/math] [math](para \rho \gt 0 ) [/math]

En la práctica vamos a usar la definición continua en [math] \rho = 0 [/math]; [math] \rho = 0 \Rightarrow e = 1 ; d(0,0)= 1 + 1 = 2 [/math]

De esta ecuación de la densidad podemos sacar las primeras conclusiones, como el factor [math] cos \theta =\frac{x}{\rho} [/math] hace que la parte exponencial sea mayor donde [math] x \gt 0 [/math], ya que el dominio siempre es positivo en [math] x \epsilon [0,4] [/math] y va creciendo rápidamente para [math] \rho [/math] grandes. Por este motivo tenemos la conclusión de que la masa será grande.

Aproximamos la integral por suma de Riemann con [math] h= \frac{1}{10} [/math] obteniendo [math] M\approx 1.912271 [/math] x [math] 10^6 [/math]

Dibujamos a continuación la densidad de la placa, considerando dominio, parámetro y mallado

%% Cálculo de rho y theta
rho = sqrt(XX.^2 + YY.^2);
theta = atan2(YY, XX);

%% Densidad
dens = 1 + exp( rho.^2 .* cos(theta) );

%% Masa mediante sumatorio numérico
area_celda = h^2;
masa = sum(dens(:)) * area_celda;

fprintf("Masa aproximada de la placa = %.6f\n", masa);

%% FIGURA 2D
figure('Name','Densidad en 2D');
imagesc(x, y, dens); 
axis xy; axis equal;
colorbar;
xlabel('x'); ylabel('y');
title('Densidad d(\rho,\theta) = 1 + e^{\rho^2 \cos\theta} (2D)');

%% FIGURA 3D
figure('Name','Densidad en 3D');
surf(XX, YY, dens, 'EdgeColor','none'); 
colormap('jet');
colorbar;
xlabel('x'); ylabel('y'); zlabel('d(x,y)');
title('Densidad d(\rho,\theta) = 1 + e^{\rho^2 \cos\theta} (3D)');
view(45,35);  % Ángulo 3D agradable

La función puede tomar valores muy grandes en la esquina derecha del rectángulo [math] x=4 [/math] porque [math] \rho ^2 cos \theta [/math] se aproxima a [math] x^2 [/math] cuando [math] y [/math] es pequeño y [math] x [/math] grande, por eso [math] e^{\rho ^2 cos \theta} [/math] puede alcanzar valores de [math] e^{16} [/math]

Físicamente este resultado puede ser muy elevado, pero matemáticamente es razonable por las razones que hemos dado.

12 Aplicaciones de las ondas longitudinales planas

12.1 Ingeniería estructural y vibraciones

Para diseñar, por ejemplo, un panel de avión, se necesita conocer cómo responden a las vibraciones y cómo estas afectan a la distribución de tensiones. Las piezas aeronáuticas como las de fuselaje, las alas o las cubiertas interiores suelen estar fabricadas generalmente de placas delgadas de materiales metálicos o compuestos. Estas deben soportar numerosos esfuerzos, como las vibraciones creadas por fuerzas variables aerodinámicas, los cambios de temperatura del aire o la fatiga mecánica que puede producir microfisuras en el material.
El estudio realizado en este trabajo es crucial para identificar las zonas de mayor tensión y para conocer como se comporta el material frente a los anteriores efectos. De esta forma, se pueden prevenir fallos y mejorar el diseño y selección de materiales previa a su fabricación.

12.2 Control térmico

Detectar los cambios de temperatura y su distribución en un material es imprescindible en el caso de los microchips o placas electrónicas, ya que afectan a la integridad física y la función electrónica de estos. El calor y las vibraciones juntas aumentan su degradación, sin embargo, realizando el estudio previo, se puede optimizar el diseño térmico y mecánico de las placas para evitar los microdesplazamientos de capas internas y así controlar la disipación de las ondas y sus tensiones internas.

12.3 Sensores de deformación térmica

Este tipo de ondas suelen utilizarse en sensores ultrasonido para determinar defectos como grietas o variaciones en materiales. El modelo estudiado funciona muy bien para la evaluación de la integridad de estructuras en ingeniería civil de forma rápida, precisa y no destructiva. Suelen utilizarse sobre vigas de acero, losas de hormigón o columnas y pilares; si existe una grieta, la onda se refleja sin causar daños en el elemento estructural. También pueden determinar el espesor de materiales, o incluso detectar vibraciones longitudinales inducidas por cargas reales, muy útiles en el caso de las presas, puentes o edificios altos. De esta forma, se crea una herramienta de inspección al controlar la onda y su dispersión o desplazamiento a lo largo del material.

13 Póster

A continuación compartimos el póster, basado en los contenidos publicados en este trabajo, en el que resumimos las labores desarrolladas por el grupo: