Ver el código de «Deformaciones de una columna.Grupo 28»

{{ TrabajoED | Deformaciones de una columna . Grupo 28 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2024-25]] |
*Mateo Navarro Díaz
*María Victoria González Junco
*Fernando Benítez Pérez
*Rodrigo Prado Fornos
*Claudia Elimar Manrique }}


[[Categoría:Teoría de Campos]]
[[Categoría:TC24/25]]


Para el siguiente articulo, consideraremos una placa rectangular plana ocupando la región en el espacio plano<br> <math>(x, y) ∈ [-1, 1]×[0, 10]</math>. A continuación definimos dos cantidades físicas; por un lado la temperatura dada <br>como 
<center><math> T(\rho, \theta)=\sin(2π\rho)</math></center>
en cartesianas como
<center><math> T(x, y)=\sin(2π\sqrt{x^2 + y^2})</math></center> 

También, definimos la densidad de la placa como <center><math> d(\rho, \theta)=(2-\rho)(4-\cos(4(\rho+π/2)))</math></center>

Por otro lado, hemos de tener en cuenta los desplazamientos <math>\vec{u}(\rho, \theta)</math> producidos por la acción de una fuerza <br>determinada. 
Estos desplazamientos se corresponden con el campo <math>\vec{u}(\rho, \theta)= \rho sin(\theta)sin(2π\rho/50)e_\vec{\theta}</math> .







==  Representación del mallado del sólido ==
[[Archivo:Mallado Placa.png|944 × 766px|miniaturadeimagen|derecha|Mallado de la placa]]
Definimos el mallado que representa el interior del sólido mediante dos vectores x e y los cuales representan el intervalo en el que está definida nuestra placa, tomando un muestreo h=1/10. Al crear la malla, al representarla con el comando ''surf'', tenemos que tener en cuenta que es plana dando valor 0 a la altura de los puntos de la misma. 

{{matlab|codigo=
h=0.1;                       %muestreo
x=-1:h:1;                    %eje x de la placa
y=0:h:10;                     %eje y de la placa
[xx,yy]=meshgrid(x,y);       %mallado de la placa
figure(1)
surf(xx,yy,0*xx);            %representación de la placa
axis([-2,2,0,10])             %ejes del rectángulo
view(2)

}}

==Gradiente de la temperatura==
El cálculo del gradiente de un campo escalar se expresa como la derivada parcial respecto de cada coordenada de dicho campo en función de la base ortonormal orientada positiva  <math>\vec{i},\vec{j},\vec{k}</math> (COORDENADAS CARTESIANAS). Es decir:
<center><math> \nabla T=\frac{\partial{T}}{\partial{x}}\vec{i}+\frac{\partial{T}}{\partial{y}}\vec{j}</math></center>


<center><math> T=\sin(2π\sqrt{x^2 + y^2})</math></center>


<center><math> \nabla T=\frac{2πx cos(2π\sqrt{x^2 + y^2})}{\sqrt{x^2 + y^2}}\vec{i}+\frac{2πy cos(2π\sqrt{x^2 + y^2})}{\sqrt{x^2 + y^2}}\vec{j}</math></center>

A continuación se muestra el código completo de matlab, del cual resultan las siguientes imágenes en las que podemos ver las curvas de nivel y donde se encuentra su máximo y mínimo, a partir de los colores de las gráficas y a partir de el propio matlab. Este nos índica al ejecutar el programa que el máximo es: 0.99996
[[Archivo:Ctempcnivel.jpg|1120 × 840px|miniaturadeimagen|C.Temperatura C.Nivel]]

{{matlab|codigo=
clc
clear
h = 0.2; %MUESTREO
x =-1.0:h:1.0; %DOMINIO DE X
y =0.0:h:10.0; %DOMINIO DE Y
[X,Y]= meshgrid(x,y); % CREACIÓN DEL MALLADO
T=sin(2.*pi.*(sqrt(X.^2+Y.^2))); %TEMPERATURA EN CADA PUNTO DEL MALLADO
mesh(X,Y,0*X) %CREACIÓN DE LA MALLA
subplot(1,3,1) %GRÁFICO SUPERIOR
surf(X,Y,T); %DEGRADACIÓN DE COLORES
axis ([-1,1,-0.5,10.5]) 
view(2)
colorbar %BARRA DE COLORES
%TÍTULO PRIMERA GRÁFICA Y EJES
title('CAMPO DE TEMPERATURAS')
xlabel('eje X')
ylabel('eje Y')
subplot(1,3,2) %GRÁFICO INFERIOR
contour(X,Y,T,11);
axis ([-1.0,1.0,-0.5,10.5]);
%TÍTULO SEGUNDA GRÁFICA Y EJES 
title('CURVAS DE NIVEL')
xlabel('eje X')
ylabel('eje Y')
colorbar %BARRA DE COLORES
Tmax=max(max(T)) %PUNTO MÁXIMO
subplot(1,3,3)
axis ([-1.25,1.25,-0.5,10.5]) 
view(2)
%CURVAS DE NIVEL ENUMERADAS
[M,c]=contour(X,Y,T,[0:0.75:10],'ShowText','on') 
%TÍTULO TERCERA GRÁFICA Y EJES
title('CURVAS DE NIVEL NUMERADAS')
xlabel('eje X')
ylabel('eje Y')
colorbar %BARRA DE COLORES
}}

A continuación, se va a demostrar como el gradiente, previamente calculado, es perpendicular a las líneas de nivel anteriores.
Al no verse claro en la imagen original, con la cantidad de líneas de flujo, se opta por aplicar un zoom para una correcta visualización.


[[Archivo:Gradtemp3.jpeg|400px|centro|thumb|Mallado de la placa]]

{{matlab|codigo=
clear;clc;
h = 2/10; %MUESTREO
x = [-1:h:1]; %DOMINIO DE X
y = [0:h:12]; %DOMINIO DE Y
[X,Y]= meshgrid(x,y); %CREACIÓN DEL MALLADO
T=sin(2.*pi.*(sqrt(X.^2+Y.^2))); %TEMPERATURA EN CADA PUNTO DEL MALLADO
contour(X,Y,T,11); %CURVAS DE NIVEL
dx=2*pi*X.*cos(2.*pi.*sqrt(X.^2+Y.*2))./sqrt(X.^2+Y.*2); %PARCIAL DE X
dy=2*pi*Y.*cos(2.*pi.*sqrt(X.^2+Y.*2))./sqrt(X.^2+Y.*2); %PARCIAL DE Y
%TÍTULO Y EJES
title('Gradiente de temperatura');
xlabel('Eje X');
ylabel('Eje Y');
% Representación de la temperatura y las curvas de nivel
hold on
quiver(x,y,dx,dy);
axis equal
colorbar
}}

==Ley de Fourier==
Gracias al segundo principio de la termodinámica conocemos el calor emitido por dos focos, uno cálido y otro frío, solo puede ir en un sentido, de mayor a menor calor.
Y la ley expresada por el matemático y físico Jean-Baptiste Joseph Fourier trata de establecer una relación entre la distancia, tiempo y el flujo de calor entre focos. Dicha ley viene expresada <math>\vec{Q}=−K∇T</math>, donde K es la constante de conductividad térmica que dependerá del material a estudio. En nuestro caso, la sección del solido tiene por  K el valor asociado 1. 

Se observa que las flechas van en sentido contrario al gradiente. 

[[Archivo:enercalorif2.png|1254 × 1247px|miniaturadeimagen|derecha|Figura 3]]

{{matlab|codigo=

clear;clc;
%Definimos las variables x e y usando como paso de muestreo h
 h=2/10;
 x=[-1:h:1];
 y=[0:h:10];
%Creación del mallado
 [X,Y]=meshgrid(x,y);
 mesh(X,Y,0.*X);
%Definimos la función T
 T=sin(2.*pi.*(sqrt(X.^2+Y.^2)));
%Curvas de nivel
 contour(X,Y,T,11);
hold on
%Calculo del gradiente de T
 dx=2*pi*X.*cos(2.*pi.*sqrt(X.^2+Y.*2))./sqrt(X.^2+Y.*2);
 dy=2*pi*X.*cos(2.*pi.*sqrt(X.^2+Y.*2))./sqrt(X.^2+Y.*2);
%Constante de conductividad térmica y ley de Fourier 
 k=1;
  Q1=-k.*(dx);
  Q2=-k.*(dy);
%Representación del campo vectorial
 quiver(x,y,Q1,Q2,'r');
 axis equal;
hold off
%Título y nombre de los ejes
 title('Energía calorífica');
  xlabel('Eje X');
  ylabel('Eje Y');
}}

== Variación de temperatura ==
La temperatura <math>T(\rho) = \sin(2\pi \rho)</math> varía cuando cambia <math>T(\rho)</math>, y su mayor cambio ocurre cuando la derivada <math>\frac{\partial T}{\partial \rho} = 2\pi \cos(2\pi \rho)</math> alcanza valores máximos o mínimos. Esto sucede cuando <math>\rho</math> es un número entero, ya que en esos puntos el coseno toma sus valores máximos.
La variación es radial porque <math>T</math> solo depende de <math>\rho</math>. El gradiente de temperatura apunta hacia afuera o hacia adentro, dependiendo del signo de <math>\cos(2\pi \rho)</math>, lo que indica la dirección de mayor cambio. Por lo tanto, la dirección de máximo crecimiento es perpendicular a las circunferencias de <math>\rho=n</math>, donde 𝑛 es un número entero. Esto es lo que se observa en las flechas ubicadas en la circunferencia de radio <math>\rho = 1</math>.

[[Archivo:Vartemp4bien.jpeg|500 px|miniaturadeimagen|derecha| Variación de temperatura. ]]

{{matlab|codigo=
% Parámetros en coordenadas cartesianas
x = linspace(-1, 1, 200);
y = linspace(0, 10, 200);
% Crear una malla para la sección
[X, Y] = meshgrid(x, y);
% Conversión a coordenadas cilindricas
R = sqrt(X.^2 + Y.^2); % rho
Theta = atan2(Y, X); % theta (arctan(y/x))
% Temperatura en coordenadas polares
T = sin(2 * pi * R);
% Gradiente de T con respecto a rho
dT_drho = 2 * pi * cos(2 * pi * R);
% Dibujar el mapa de temperatura en coordenadas cartesianas
figure;
contourf(X, Y, T, 20, 'LineColor', 'none'); % Mapa de temperatura
colorbar;
hold on;
% Definir las circunferencias de máxima variación
max_rhos = 1:10;
for rho_val = max_rhos
    num_puntos = 300; 
    theta_circle = linspace(0, 2*pi, num_puntos);
    x_circle = rho_val * cos(theta_circle); 
    y_circle = rho_val * sin(theta_circle);
    valid_idx = (x_circle >= -1 & x_circle <= 1 & y_circle >= 0 & y_circle <= 10);
    x_circle = x_circle(valid_idx);
    y_circle = y_circle(valid_idx);
    
    % Marcar los puntos de máxima variación con rojo
    plot(x_circle, y_circle, 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 2);
    
    
    if rho_val == 1
        num_flechas = 20; 
        theta_flechas = linspace(0, 2*pi, num_flechas + 1); 
        
        for t = theta_flechas(1:end-1)
            x_start = rho_val * cos(t); 
            y_start = rho_val * sin(t); 
            dx = 0.4 * cos(t);
            dy = 0.4 * sin(t); 
            quiver(x_start, y_start, dx, dy, 'r', 'LineWidth', 1, 'MaxHeadSize', 1.5);
        end
    end
end
title('Temperatura y Dirección de Máxima Variación (Coordenadas Cartesianas)');
xlabel('x');
ylabel('y');
axis equal;
xlim([-1.5, 1.5]);
ylim([0, 10]);
hold off;
}}

== Campo de desplazamientos ==

Los puntos de desplazamiento nulo son aquellos en los que el campo de desplazamiento se hace cero, es decir, donde el valor de <math>u(\rho, \theta)</math> es igual a 0.
Dado que el desplazamiento en coordenadas cilíndricas está dado por:
<center><math> u(\rho, \theta) = \rho \sin(\theta) \sin\left( \frac{50}{2\pi \rho} \right)</math></center>
Las condiciones completas para los puntos de desplazamiento nulo son:
* Cuando <math>\rho =0 </math>es decir, el origen de coordenadas (<math>x =0</math> <math>y =0</math>.

* Cuando <math>\sin(\theta) = 0</math>, es decir, cuando <math>\theta =0 </math> o <math>\theta = \pi</math>, esto sucede en cuando <math>y = 0</math>).

* Cuando <math>\sin\left( \frac{50}{2 \pi \rho} \right) = 0</math>, esto ocurre cuando <math>\rho</math> es multiplo de 50.

 



[[Archivo:Campodesplazcolumnag28.jpeg|500px||miniaturadeimagen|derecha| Campo de desplazamientos y puntos fijos. ]]

{{matlab|codigo=
% Definir la malla con menos flechas pero mayor densidad de puntos rojos
h_flechas = 0.3;
h_puntos_rojos = 0.05; 

% Crear la malla en coordenadas cartesianas para las flechas
x_flechas = -1:h_flechas:1;  
y_flechas = 0:h_flechas:10;  
x_rojos = -1:h_puntos_rojos:1; 
y_rojos = zeros(size(x_rojos));
[Mx_flechas, My_flechas] = meshgrid(x_flechas, y_flechas);

% Convertir a coordenadas cilíndricas
rho = sqrt(Mx_flechas.^2 + My_flechas.^2);  
theta = atan2(My_flechas, Mx_flechas);  

% Desplazamiento en coordenadas cilíndricas
u_theta = rho .* sin(theta) .* sin((2 * pi * rho) / 50);
% Convertir a coordenadas cartesianas
up = u_theta .* (-sin(theta)); 
uth = u_theta .* cos(theta); 
% Graficar el campo de desplazamientos
figure;
h_quiver = quiver(Mx_flechas, My_flechas, up, uth, 1.5, 'Color', [0, 0, 1], 'AutoScale', 'on', 'LineWidth', 1.0);  
h_quiver.AutoScaleFactor =2.8;  
hold on;
% Graficar la región [−1, 1] × [0, 10]
plot([-1 1], [0 0], 'k', 'LineWidth', 2); 
plot([-1 1], [10 10], 'k', 'LineWidth', 2);  
plot([-1 -1], [0 10], 'k', 'LineWidth', 2);  
plot([1 1], [0 10], 'k', 'LineWidth', 2);  
% Etiquetas y título
axis equal;
axis([-2, 2, -0.5, 12.5]);
xlabel('Eje de las X');
ylabel('Eje de las Y');
title('Campo de Desplazamientos de la Columna');
% Detectar puntos fijos (donde los desplazamientos son cero)
% 1. Cuando y = 0, sin(theta) = 0 (θ = 0 y θ = π), es decir, para todos los x en [-1, 1]
fixed_points = (My_flechas == 0); 
% 2. También cuando sin(2πρ/50) = 0, es decir, cuando ρ = 50n, donde n es un número entero
fixed_rho_points = mod(rho, 50) == 0;  
% Combinamos ambas condiciones para encontrar los puntos nulos de desplazamiento
combined_fixed_points = fixed_points | fixed_rho_points;
% Graficar los puntos fijos (solo aquellos donde u_theta es cero)
h_points = plot(x_rojos, y_rojos, 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 4); % Mostrar los puntos rojos con alta densidad
legend([h_quiver, h_points], {'Campo de desplazamientos', 'Puntos fijos'}, 'Location', 'best');
grid on;
hold off;
}}

== Sólido antes y después del desplazamiento ==

En esta ocasión representaremos el desplazamiento de la columna definido por el campo vectorial <math>u(\rho, \theta) = \rho \sin(\theta) \sin\left( \frac{2\pi\rho}{50} \right) e^\theta</math>.

Hemos ejecutado el mallado para generar los puntos de dominio, aplicamos la fórmula del campo vectorial a dicho mallado para obtener la nueva posición del sólido,y por último, nos hemos apoyado en subgráficos para comparar el sólido antes y después del desplazamiento.

[[Archivo:Antesydespuescolumna28.jpeg|500px||miniaturadeimagen|derecha| Solido antes y después del desplazamiento. ]]

{{matlab|codigo=
% Parametrización 
h = 0.1; 
rho = -1:h:1; 
tt = 0:h:10;

% Mallado
[Mrho, Mtt] = meshgrid(rho, tt);

% Campo de desplazamientos 
u_theta = (Mrho .* sin(Mtt) .* sin(Mrho .* 2 * pi ./ 50)); 
rdrho = Mrho; 
rdtt = Mtt + u_theta;

% Grafica sólido antes del desplazamiento
subplot(1, 2, 1);
mesh(Mrho, Mtt, zeros(size(Mrho))); 
axis equal;
axis([-1.5, 1.5, 0, 10]); 
title('Sólido antes del desplazamiento');
xlabel('(Eje X)');
ylabel('(Eje Y)');
view(2); 
grid on;

% Grafica sólido después del desplazamiento
subplot(1, 2, 2);
mesh(rdrho, rdtt, zeros(size(rdrho))); 
axis equal;
axis([-1.5, 1.5, 0, 10]); 
title('Sólido después del desplazamiento');
xlabel('(Eje X)');
ylabel('(Eje Y)');
view(2); 
grid on;

}}

== Divergencia del campo de desplazamientos ==

La divergencia de un campo vectorial es una medida de la cantidad de flujo que sale o entra en un punto de la región definida.


En este caso la divergencia resultante es <math>
\nabla\cdot\vec u = cosθsin( \frac{2πρ}{50})\
</math>

Puntos mínimo y máximo de la divergencia : <math> \frac{\partial \nabla\cdot\vec u}{\partial θ}=0 </math> ; <math> \frac{\partial \nabla\cdot\vec u}{\partial ρ }=0 </math>

Mínimo(θ,ρ)=(0,-1) ; Máximo(θ,ρ)=(0,1)

Punto nulo de la divergencia : <math>  \nabla\cdot\vec u = 0 </math> ; si ρ=0 el <math> sin( \frac{2πρ}{50})\ = 0 </math>  lo que anularía el valor de la divergencia.


[[Archivo:Divergenciacolumna.g28.jpg|600 px|miniaturadeimagen|derecha|Divergencia]]

[[Archivo:Divergenciacolumna2.g28.jpeg|600 px|miniaturadeimagen|derecha|Divergencia]]

{{matlab|codigo=

% Definición de la malla en coordenadas cartesianas
x = -1:0.2:1;         % Coordenadas x
y = 0:0.1:10;         % Coordenadas y
[xx, yy] = meshgrid(x, y);

% Conversión de la malla cartesiana (x, y) a coordenadas cilíndricas (r, theta)
r = sqrt(xx.^2 + yy.^2);         % Radio en coordenadas cilíndricas
theta = atan2(yy, xx);           % Ángulo en coordenadas cilíndricas

% Definir la divergencia en coordenadas cilíndricas
diver = cos(theta) .* sin(2 * pi * r / 50);

% Encontrar los puntos de divergencia máxima, mínima y nula
[minVal, minIdx] = min(diver(:));       % Valor mínimo y su índice lineal
[maxVal, maxIdx] = max(diver(:));       % Valor máximo y su índice lineal
nullIdx = find(abs(diver) < 1e-6);      % Índices donde la divergencia es nula (aproximación)

% Convertir los índices lineales a índices de matriz (para la malla)
[rowMin, colMin] = ind2sub(size(diver), minIdx); % Índices del mínimo
[rowMax, colMax] = ind2sub(size(diver), maxIdx); % Índices del máximo
[rowNull, colNull] = ind2sub(size(diver), nullIdx); % Índices donde la divergencia es nula

% Coordenadas correspondientes en (x, y)
minPoint = [xx(rowMin, colMin), yy(rowMin, colMin)];
maxPoint = [xx(rowMax, colMax), yy(rowMax, colMax)];
nullPoints = [xx(rowNull, colNull), yy(rowNull, colNull)];

% Mostrar los resultados
fprintf('Divergencia mínima:\n Valor: %.3f, Punto: (%.3f, %.3f)\n', minVal, minPoint(1), minPoint(2));
fprintf('Divergencia máxima:\n Valor: %.3f, Punto: (%.3f, %.3f)\n', maxVal, maxPoint(1), maxPoint(2));
fprintf('Divergencia nula en %d puntos:\n', size(nullPoints, 1));
disp(nullPoints);

% Representación de la divergencia
figure;
surf(xx, yy, diver);             % Malla original, pero con valores de divergencia corregidos
shading interp;                  % Suavizado
colorbar;                        % Barra de colores
hold on;

% Marcar los puntos de interés
plot3(minPoint(1), minPoint(2), minVal, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % Mínimo
plot3(maxPoint(1), maxPoint(2), maxVal, 'go', 'MarkerSize', 8, 'MarkerFaceColor', 'g'); % Máximo
for i = 1:size(nullPoints, 1)  % Marcar puntos nulos
    plot3(nullPoints(i, 1), nullPoints(i, 2), 0, 'bo', 'MarkerSize', 6, 'MarkerFaceColor', 'b');
end

% Etiquetas del gráfico
title('Divergencia en coordenadas cilíndricas con puntos de interés');
xlabel('x');
ylabel('y');
zlabel('Divergencia');
legend('Divergencia', 'Mínimo', 'Máximo', 'Nulo');
}}

== Rotacional del campo ==

*<math>\triangledown \times \vec{u}(\rho,\theta,z)=\frac{1}{\rho}\begin{vmatrix}
\vec{e}_{\rho} & \rho \vec{e}_{\theta } & \vec{e}_{z}\\ 
\frac{\partial }{\partial \rho } & \frac{\partial }{\partial \theta } & \frac{\partial }{\partial z}\\ 
u_{\rho } & \rho u_{\theta } & u_{z}
\end{vmatrix} </math> 

Aplicando nuestro campo escalar, obtenemos:

*<math>\triangledown \times \vec{u}(\rho,\theta,z)=\frac{1}{\rho}\begin{vmatrix}
\vec{e}_{\rho} & \rho \vec{e}_{\theta } & \vec{e}_{z}\\ 
\frac{\partial }{\partial \rho } & \frac{\partial }{\partial \theta } & \frac{\partial }{\partial z}\\ 
0 & p^2 sin(\theta)sin( \frac{2πρ}{50}) & 0
\end{vmatrix} = (2sin(\theta)sin( \frac{2πρ}{50}) + psin(\theta)cos( \frac{2πρ}{50})*\frac{2πρ}{50})\vec{e}_{z}\\ </math>


*<math>\left \|\triangledown \times \vec{u}(\rho,\theta,z)\right \| = sin(\theta)\sqrt{4sin^2(\frac{2πρ}{50})+p^2cos^2(\frac{2πρ}{50})*\frac{4π^2}{50^2}+4psin(\frac{2πρ}{50})cos(\frac{2πρ}{50})*\frac{2π}{50}} </math>


[[Archivo:Rotacionalcolumnag28.jpeg|600px||miniaturadeimagen|derecha| Rotacional ]]
{{matlab|codigo=

% Dimensiones de la sección transversal
x = -1:0.1:1;  % Coordenadas x
y = 0:0.1:10;  % Coordenadas y
[xx, yy] = meshgrid(x, y);

% Convertir coordenadas cartesianas a cilíndricas
r = sqrt(xx.^2 + yy.^2);         % Radio en coordenadas cilíndricas
theta = atan2(yy, xx);           % Ángulo en radianes (atan2 para cuadrante correcto)

% Definir los límites de r y theta
r(r == 0) = eps;  % Evitar división por cero en r

% Cálculo del rotacional
% (nabla x u)_z = 2 * sin(theta) * sin(2*pi*r/50) + (2*pi*r/50) * sin(theta) * cos(2*pi*r/50)
rot_z = 2 .* sin(theta) .* sin(2 * pi * r / 50) + (2 * pi * r / 50) .* sin(theta) .* cos(2 * pi * r / 50);

% Calcular el valor absoluto del rotacional
abs_rot_z = abs(rot_z);

% Representación gráfica del valor absoluto del rotacional
figure;
contourf(xx, yy, abs_rot_z, 20, 'LineColor', 'none');  % Gráfico de contornos rellenos
colorbar;                                             % Barra de colores
title('Valor absoluto del rotacional (∇ × ⃗u_z) en la sección transversal');
xlabel('x (m)');
ylabel('y (m)');

% Añadir información sobre los ejes
hold on;
quiver(xx, yy, -yy ./ r, xx ./ r, 0.5, 'k');  % Direcciones del campo vectorial (rotacional)
hold off;
}}

== Tensiones Normales ==
 El tensor de tensiones <math>σ_{ij}</math> se puede escribir a través de la fórmula: <center> <math>σ = λ∇ · \vec{u} 1 + 2µϵ</math></center>
Donde <math>ϵ(\vec{u})= (∇\vec{u} + ∇\vec{u}^t)/2  </math>,      es la parte simétrica del tensor gradiente de <math> \vec{u}</math>.



Empezaremos calculando  <center><math>ε(\vec u)=
\begin{pmatrix}
0 & 0 & 0\\ 0 & cos(\theta)sin(\frac{2\pi*p}{50}) & 0 \\ 0 & 0 & 0
\end{pmatrix}</math> </center>
<br />
<center><math>σ=λ\nabla·\vec u 1+2με = λ*cos(\theta)sin(\frac{2\pi*p}{50})·\begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} + 2μ\begin{pmatrix}0 & 0 & 0 \\ 0 & cos(\theta)sin(\frac{2\pi*p}{50}) & 0 \\ 0 & 0 & 0 \end{pmatrix} </math>



<math>= cos(\theta)sin(\frac{2\pi*p}{50}) (\begin{pmatrix} λ & 0 & 0\\ 0 & λ & 0 \\ 0 & 0 & λ \end{pmatrix} + \begin{pmatrix} 0 & 0 & 0\\ 0 & 2μ & 0 \\ 0 & 0 & 0 \end{pmatrix}) = cos(\theta)sin(\frac{2\pi*p}{50})\begin{pmatrix} λ & 0 & 0\\ 0 & 2μ+λ & 0 \\ 0 & 0 & λ \end{pmatrix}→{λ=μ=1}→cos(\theta)sin(\frac{2\pi*p}{50})\begin{pmatrix} 1 & 0 & 0\\ 0 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix}</math></center>


Las tensiones normales en la dirección de los ejes <math>(\vec e_ρ,\vec e_\theta,\vec e_z)</math> son:

::* <math>\vec e_ρ · σ ·\vec e_\rho = \begin{pmatrix} 1 & 0 & 0 \end{pmatrix}·\begin{pmatrix} cos(\theta)sin(\frac{2\pi*p}{50}) & 0 & 0\\ 0 & 3*cos(\theta)sin(\frac{2\pi*p}{50}) & 0 \\ 0 & 0 & cos(\theta)sin(\frac{2\pi*p}{50}) \end{pmatrix}·\begin{pmatrix} 1 \\0 \\ 0 \end{pmatrix} = cos(\theta)sin(\frac{2\pi*\rho}{50}) </math>


::*<math>\vec e_\theta · σ ·\vec e_\theta = \begin{pmatrix} 0 & 1 & 0 \end{pmatrix}·\begin{pmatrix} cos(\theta)sin(\frac{2\pi*p}{50}) & 0 & 0\\ 0 & 3*cos(\theta)sin(\frac{2\pi*p}{50}) & 0 \\ 0 & 0 & cos(\theta)sin(\frac{2\pi*p}{50}) \end{pmatrix}·\begin{pmatrix} 0 \\1 \\ 0 \end{pmatrix} = 3*cos(\theta)sin(\frac{2\pi*\rho}{50}) </math>


::* <math>\vec e_z · σ ·\vec e_z = \begin{pmatrix} 0 & 0 & 1 \end{pmatrix}·\begin{pmatrix} cos(\theta)sin(\frac{2\pi*p}{50}) & 0 & 0\\ 0 & 3*cos(\theta)sin(\frac{2\pi*p}{50}) & 0 \\ 0 & 0 & cos(\theta)sin(\frac{2\pi*p}{50}) \end{pmatrix}·\begin{pmatrix} 0 \\0 \\ 1 \end{pmatrix} = cos(\theta)sin(\frac{2\pi*\rho}{50}) </math>

== Tensiones tangenciales al plano ortogonal al \( \vec{e_\rho} \)==
Para calcular las tensiones tangenciales respecto al plano ortogonal a \( \vec{e_\rho} \), utilizaremos la fórmula:<math>\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| </math>. <br />  en la que sustituiremos \(\vec{i}\) por \( \vec{e_\rho} \)
 
Entonces para comenzar a operar sabemos que el tensor tensiones es el:
σ=\begin{pmatrix}\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right) & 0 & 0 \\ 0 & 3\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right) & 0 \\ 
0 & 0 & \cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right)
\end{pmatrix}.

Primero, calculamos el producto \( \sigma \cdot \vec{e_\rho} \):=\begin{pmatrix}\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right) \\ 0 \\ 0\end{pmatrix}
y luego el segundo sumando \((\vec{e_\rho} \cdot \sigma \cdot \vec{e_\rho}) \vec{e_\rho}\)=
\begin{pmatrix}
\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right) \\ 
0 \\ 
0
\end{pmatrix}


Finalmente, restamos para obtener las tensiones por medio de la fórmula \( \vec{\tau} \)=
\begin{pmatrix}
\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right) \\ 
0 \\ 
0
\end{pmatrix}
-
\begin{pmatrix}
\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right) \\ 
0 \\ 
0
\end{pmatrix}
=
\begin{pmatrix}
0 \\ 
0 \\ 
0
\end{pmatrix}.


Por lo tanto, no existen tensiones tangenciales no nulas en este caso.



== Tensión de Von Mises ==

La tensión de Von Mises se define mediante la fórmula: <br /> <math> \sigma_{VM} = \sqrt{ \frac { ( \sigma_1 - \sigma_2 )^2 + ( \sigma_2 - \sigma_3 )^2 + ( \sigma_3 - \sigma_1 )^2  }{2} }</math>;<br />
donde \( \sigma_1 \), \( \sigma_2 \), y \( \sigma_3 \) son los autovalores del tensor de tensiones \( \sigma \).

Este resultado permite calcular la tensión de Von Mises indicando los puntos donde la tensión alcanza valores máximos y así podemos determinar si un material soportará las cargas aplicadas o si llegará a fallar.


[[Archivo:VonMisesapartado11campos.jpeg|450px||miniaturadeimagen|derecha| Tensión Von Mises ]]

{{matlab|codigo=
% Definir la sección transversal
x = linspace(-1, 1, 200);
y = linspace(0, 10, 500);
[X, Y] = meshgrid(x, y);
rho = sqrt(X.^2 + Y.^2);
theta = atan2(Y, X);

% Componentes del tensor de tensiones en coordenadas polares
T_rr = rho .* cos(theta) .* sin(2 * pi * rho / 50);
T_rtheta = 3 * cos(theta) .* sin(2 * pi * rho / 50);
T_zz = T_rr;

% Inicializar los autovalores
sigma_1 = zeros(size(X));
sigma_2 = zeros(size(X));
sigma_3 = zeros(size(X));

% Calcular los autovalores en cada punto
for i = 1:size(X, 1)
    for j = 1:size(X, 2)
        % Construir el tensor de tensiones
        stress_tensor = [T_rr(i, j), T_rtheta(i, j), 0;
                         T_rtheta(i, j), T_rr(i, j), 0;
                         0, 0, T_zz(i, j)];
                     
        % Calcular los autovalores
        eigenvalues = eig(stress_tensor);
        sigma_1(i, j) = eigenvalues(1);
        sigma_2(i, j) = eigenvalues(2);
        sigma_3(i, j) = eigenvalues(3);
    end
end

% Calcular la tensión de Von Mises
VM = sqrt(0.5 * ((sigma_1 - sigma_2).^2 + (sigma_2 - sigma_3).^2 + (sigma_3 - sigma_1).^2));

% Encontrar el valor máximo y su posición
[VM_max, idx] = max(VM(:));
[max_row, max_col] = ind2sub(size(VM), idx);
max_x = X(max_row, max_col);
max_y = Y(max_row, max_col);

% Graficar en 3D
figure;
surf(X, Y, VM, 'EdgeColor', 'none'); % Gráfico de superficie
colormap('parula');
colorbar;
hold on;

% Asegurar que el punto rojo aparece en la gráfica
z_offset = VM_max + 0.01 * (max(VM(:)) - min(VM(:))); % Evitar problemas de superposición
plot3(max_x, max_y, z_offset, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');

% Añadir texto para mayor claridad
text(max_x, max_y, z_offset, sprintf('VM Máx: %.2f', VM_max), 'Color', 'red', 'FontSize', 10);

% Personalizar el gráfico
xlabel('x');
ylabel('y');
zlabel('Tensión de Von Mises (VM)');
title('Distribución de Tensión de Von Mises en 3D');
legend(sprintf('Máximo VM: %.2f en (%.2f, %.2f)', VM_max, max_x, max_y), 'Location', 'best');
grid on;
view(45, 30); % Ajustar ángulo de la vista

}}

== Campo de fuerzas ==
El desplazamiento de la placa viene dado por un vector fuerza <math>\ \vec F </math>  que actua en el interior de la placa, El cálculo de dicho vector de fuerzas se realiza haciendo una aproximación de la ecuación de la elasticidad lineal, que es : <math>\ \vec F = \nabla ·\sigma_{ij} </math>, es decir que habrá que aplicar la divergencia a la matriz del antes mencionado tensor tensiones  
<math> \ \vec{F}_\rho = -\frac{\partial \sigma_{\rho \rho}}{\partial \rho} \cdot \vec{e}_\rho - \frac{1}{\rho}\frac{\partial \sigma_{\rho \theta}}{\partial \theta} \cdot \vec{e}_\rho - \frac{\sigma_{\theta \theta}}{\rho} \cdot \vec{e}_\rho = -\left( \frac{\partial \cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right)}{\partial \rho} + \frac{1}{\rho} \frac{\partial 0}{\partial \theta} + \frac{3\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right)}{\rho} \right) \cdot \vec{e}_\rho = -\left(\frac{2\pi}{50}\cos(\theta) \cos\left(\frac{2\pi \rho}{50}\right) + \frac{3\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right)}{\rho}\right) \cdot \vec{e}_\rho </math><br />


<math> \ \vec{F}_\theta = -\frac{\partial \sigma_{\theta \rho}}{\partial \rho} \cdot \vec{e}_\theta - \frac{1}{\rho}\frac{\partial \sigma_{\theta \theta}}{\partial \theta} \cdot \vec{e}_\theta - \frac{\sigma_{\rho \theta}}{\rho} \cdot \vec{e}_\theta = -\left( \frac{\partial 0}{\partial \rho} + \frac{1}{\rho} \frac{\partial 3\cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right)}{\partial \theta} + \frac{0}{\rho} \right) \cdot \vec{e}_\theta = -\left(\frac{-3\sin(\theta) \sin\left(\frac{2\pi \rho}{50}\right)}{\rho}\right) \cdot \vec{e}_\theta </math><br />

<math> \ \vec{F}_z = -\frac{\partial \sigma_{z \rho}}{\partial \rho} \cdot \vec{e}_z - \frac{1}{\rho}\frac{\partial \sigma_{z \theta}}{\partial \theta} \cdot \vec{e}_z - \frac{\partial \sigma_{z z}}{\partial z} \cdot \vec{e}_z = -\left( \frac{\partial 0}{\partial \rho} + \frac{1}{\rho} \frac{\partial 0}{\partial \theta} + \frac{\partial \cos(\theta) \sin\left(\frac{2\pi \rho}{50}\right)}{\partial z} \right) \cdot \vec{e}_z = 0 \cdot \vec{e}_z </math><br />


Después de hacer la divergencia vemos que la fuerza no tiene componente \( \vec{F}_z \) ya que es nula es por eso que solo actuará en  \( \vec{F}_\rho \) y \( \vec{F}_\theta \)

== Masa Total de la Placa == 
Queremos calcular la masa de la placa rectangular [−1,1]×[0,10]. La masa se obtiene integrando el campo escalar: <math> d(\rho, \theta)=(2-\rho)(4-\cos(4(\theta+π/2)))</math>  (densidad de la placa) a lo largo de la superficie.
Como la placa es rectangular, la podemos parametrizar utilizando coordenadas cartesianas.<math> d(x, y)=(2-\sqrt{x^2+y^2})(4-\cos(4(arctan(\frac{x}{y})+π/2)))</math>


<math> Masa=\int_{-1}^{1} \int_{0}^{10}  (2-\sqrt{x^2+y^2})(4-\cos(4(arctan(\frac{x}{y})+π/2)))dxdy = 184.771 </math>