Diferencia entre revisiones de «Estudio de la temperatura y deformación sobre una placa plana en forma de sección de anillo. (Grupo 31A)»

Revisión actual del 21:27 17 dic 2023

En este articulo se analizaran los efectos de distintos campos sobre una placa plana, entra otras cuestiones. Dicha placa tendrá la forma de una sección de anillo circular, la cual será correctamente definida en el primer apartado. Los campos que se estudiaran serán dos: un campo escalar, la temperatura; y un campo vectorial, la deformación de la placa.

Para una mejor representación visual de los conceptos expuestos en el proyecto se usara el apoyo del software Matlab. Todo el código utilizado durante el trabajo será debidamente expuesto al final del mismo. La organización se realizara en sucesivos apartados respondiendo a la preguntas propuestas por el profesorado, analizando las soluciones y sacando conclusiones sobre estas.

1 Sección de anillo

Empezaremos delimitando la sección sobre la que se va a trabajar, una placa plana que está definida por dos anillos de radios 1 y 2 que se corta sobre un plano [math]y≥|x|/2 [/math], esto se podrá observar en la Figura 1. En este caso se usaran coordenadas cilíndricas, donde nuestras variables serán: [math]\rho[/math] y [math]\theta [/math] que se encontraran entre los valores: [math]1 \leq \rho \leq 2[/math] y [math]0,464 \leq \theta \leq 2,678[/math]. La sección del anillo se puede observar con una z genérica, la placa plana se puede en dos dimensiones, en los ejes X e Y que tienen esta delimitada en los intervalos [-3,3]x[-1,3] y en tres dimensiones donde el eje Z no esta delimitada.

Figura 1. Mallado de la sección de anillo

2 Curvas de Nivel de la Temperatura

La temperatura en el semi-anillo está definida por el campo escalar: [math]T(x, y) = \cos ((y − 3)^2 + x) [/math], para mantener la coherencia en el trabajo, pasaremos el campo de temperaturas a coordenadas cilíndricas, de tal manera que quedaría así:

\begin{eqnarray*} %x = ρcosθ %y = ρsenθ %z = z x &=& \rho \cos \theta \\ y &=& \rho \sin \theta \\ z &=& z \end{eqnarray*}

De tal forma que la temperatura se definiese [math] T(\rho ,\theta )=\cos (\left( \rho \sin \theta -3\right) ^{2}+\rho \cos \theta )[/math]. Ahora se podrá visualizar en la figura 2 las curvas de nivel del campo de temperaturas sobre la placa, de la cual se descubre que la temperatura máxima que se alcanza es de: 0.99904 ºC

Figura 2. La temperatura y sus curvas de nivel

Figura3. La temperatura en tres dimensiones

2.1 Gradiente

Después de se procederá a calcular el gradiente del campo de temperaturas ([math]\nabla T [/math]) ,es decir, un campo vectorial que se podrá ver en la Figura 4

\begin{eqnarray*} \nabla T (\rho ,\theta ), \ \ &=&\left \{ \left(- 2\left( \rho \sin \theta -3\right) \sin \theta +\cos \theta \right) \sin (\left( \rho \sin \theta -3\right) ^{2}+\rho \cos \theta )\right \} \cdot \overrightarrow{e}_{\rho } \\ &&+\left \{ \left( -2\left( \rho \sin \theta -3\right) \cos \theta -\sin \theta \right) \sin (\left( \rho \sin \theta -3\right) ^{2}+\rho \cos \theta )\right \} \cdot \overrightarrow{e}_{\theta } \end{eqnarray*}

Figura 4. Gradiente de la temperatura

3 Ley de Fourier

Para finalizar este tema se calcular la energía calorífica mediante la ley de Fourier, donde k=1 y que se podrá apreciar en la Figura 5. En este caso, Q es la transferencia de energía en forma de calor y k la constante de conductividad térmica del material del cual está hecha la placa.

\begin{eqnarray*} \bar Q &=& -k \cdot \nabla T (\rho ,\theta )\\ &=&\left \{ \left( 2\left( \rho \sin \theta -3\right) \sin \theta +\cos \theta \right) \sin (\left( \rho \sin \theta -3\right) ^{2}+\rho \cos \theta )\right \} \cdot \overrightarrow{e}_{\rho } \\ &&+\left \{ \left( 2\left( \rho \sin \theta -3\right) \cos \theta -\sin \theta \right) \sin (\left( \rho \sin \theta -3\right) ^{2}+\rho \cos \theta )\right \} \cdot \overrightarrow{e}_{\theta } \end{eqnarray*}

Figura 5. Ley de Fourier sobre la placa

4 Campo vectorial \begin{eqnarray*} \bar u \end{eqnarray*}

El campo vectorial \begin{eqnarray*} \bar u \end{eqnarray*} representa los desplazamientos producidos por la acción de una fuerza determinada. A continuación, en la Figura 6 se representará el solido antes y después del desplazamiento dado por el campo de vectores [math]\bar u\ [/math] en t=0:

\begin{equation*} \overrightarrow{u}\left( \rho ,\theta \right) =\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \overrightarrow{e}_{\theta }. \end{equation*} \begin{eqnarray*} \end{eqnarray*}

Figura 6. Comparación de la deformación de la placa

4.1 Divergencia

Después, se calculará la divergencia del campo vectorial \begin{eqnarray*} \bar u \end{eqnarray*} (la cual es una medida del cambio de volumen local debido al desplazamiento) y se representa en las Figuras 7 y 8 para t = 0. Los puntos en los que la divergencia de \begin{eqnarray*} \bar u \end{eqnarray*} es máxima, mínima y nula son: el P1 (0,0) es el punto donde la divergencia es nulo, el P2 (1.3738, 1.4535) es el punto de la divergencia es máxima y el P3(-1.534, 1.2833) es el punto de la divergencia es mínima.

\begin{eqnarray*} \nabla \cdot \overrightarrow{u} &=&\frac{1}{\rho }\left( \frac{\partial }{% \partial \rho }\left( \rho u_{\rho }\right) +\frac{\partial }{\partial \theta }\left( u_{\theta }\right) +\frac{\partial }{\partial z}\left( \rho u_{z}\right) \right) \\ &=&\frac{1}{\rho }\left( 0+\frac{1}{2}e^{\rho -1}2\cos \left( 2\theta -\pi /2\right) +0\right) \\ &=&\frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \end{eqnarray*}

Figura8. Divergencia en tres dimensiones

Figura 7. Divergencia de la placa

4.2 Rotacional

Por último, calcularemos el valor absoluto del rotacional de \begin{eqnarray*} \bar u \end{eqnarray*} que se podra observar el resultado en las Figuras 9 y 10. Podemos concluir diciendo que los siguientes puntos P1(-0.3482, 1.9694), P2(-0.1168, 1.9965) y P3(0.1160, 1.9966) son los que sufren un mayor rotacional.

\begin{eqnarray*} |\nabla \times \overrightarrow{u}\left( \rho ,\theta ,z\right) | &=&\frac{1}{% \rho }\left \vert \begin{array}{ccc} \overrightarrow{e}_{\rho } & \rho \overrightarrow{e}_{\theta } & \overrightarrow{e}_{z} \\ \frac{\partial }{\partial \rho } & \frac{\partial }{\partial \theta } & \frac{\partial }{\partial z} \\ u_{\rho } & \rho u_{\theta } & u_{z}% \end{array}% \right \vert =\frac{1}{\rho }\left \vert \begin{array}{ccc} \overrightarrow{e}_{\rho } & \rho \overrightarrow{e}_{\theta } & \overrightarrow{e}_{z} \\ \frac{\partial }{\partial \rho } & \frac{\partial }{\partial \theta } & \frac{\partial }{\partial z} \\ 0 & \rho \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 0% \end{array}% \right \vert \\ &=&\frac{1}{\rho }\left( 0\cdot \overrightarrow{e}_{\rho }+0\cdot \rho \overrightarrow{e}_{\theta }+\frac{\partial }{\partial \rho }\left \{ \rho \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \right \} \cdot \overrightarrow{e}_{z}\right) \\ &=&\frac{1}{\rho }\left \{ \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) +\rho \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \right \} \cdot \overrightarrow{e}_{z} \\ &=&\frac{1+\rho }{2\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \cdot \overrightarrow{e}_{z} \end{eqnarray*}

El resultado del rotacional es un vector, pero tras calcular su valor absoluto, tenemos un escalar:

\begin{eqnarray*} &=&\frac{1+\rho }{2\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \cdot \end{eqnarray*}

Figura10. Rotacional en tres dimensiones

Figura 9. Rotacional sobre la placa

5 Tensor de tensiones

En este apartado, calcularemos las tensiones normales en la dirección de ρ, θ y z las cuales se pueden ver gráficamente en las Figuras 11 y 12. Para ello, calcularemos el tensor de tensiones [math]\sigma (\vec{u}) [/math], el cual en un medio elástico lineal, isótropo y homogéneo se puede escribir como se especifica más abajo:

• [math]\overrightarrow{u}\left( \rho ,\theta \right) =\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \overrightarrow{e}% _{\theta }[/math]
• [math] \sigma =\lambda \triangledown \cdot \vec{u}\cdot 1 + 2\mu \epsilon [/math]

Para calcular el tensor de tensiones, [math]\sigma (\vec{u}) [/math], primero necesitamos calcular el tensor de deformaciones, [math]\epsilon (\vec{u}) [/math]. Este es igual a la parte simétrica del tensor gradiente del campo vectorial \begin{eqnarray*} \bar u \end{eqnarray*}. A continuación se muestra el proceso para calcularlo:

• [math]\epsilon (\vec{u})=\frac{\triangledown \vec{u}+\triangledown \vec{u}^t}{2}[/math]

Empezamos por el gradiente del campo vectorial:

\begin{eqnarray*} \nabla \overrightarrow{u}\left( \rho ,\theta \right) &=&\left( \ \frac{% \partial \overrightarrow{u}}{\partial \rho }\ \ \vdots \ \ \frac{1}{\rho }% \frac{\partial \overrightarrow{u}}{\partial \theta }\ \ \vdots \ \ \frac{% \partial \overrightarrow{u}}{\partial z}\ \right) \\ &=&\left( \begin{array}{ccc} \frac{\partial \overrightarrow{u}}{\partial \rho } \cdot \overrightarrow{e}_{\rho } & \frac{1}{\rho }\frac{\partial \overrightarrow{u}}{\partial \theta } \cdot \overrightarrow{e}_{\rho } & \frac{\partial \overrightarrow{u}}{\partial z} \cdot \overrightarrow{e}_{\rho } \\ \frac{\partial \overrightarrow{u}}{\partial \rho } \cdot \overrightarrow{e}_{\theta } & \frac{1}{\rho }\frac{\partial \overrightarrow{u}}{\partial \theta }\cdot \overrightarrow{e}_{\theta } & \frac{\partial \overrightarrow{u}}{\partial z} \cdot \overrightarrow{e}_{\theta }\\ \frac{\partial \overrightarrow{u}}{\partial \rho }\cdot \overrightarrow{e}_{z }& \frac{1}{\rho }\frac{\partial \overrightarrow{u}}{\partial \theta }\cdot \overrightarrow{e}_{z } & \frac{\partial \overrightarrow{u}}{\partial z} \cdot \overrightarrow{e}_{z } \end{array} \right) \end{eqnarray*}

\begin{eqnarray*} \frac{\partial \overrightarrow{e}_{i}}{\partial x_{j}} &=&\Gamma _{ij}^{k}\cdot \overrightarrow{e}_{k}=\Gamma _{ij}^{1}\cdot \overrightarrow{e}_{1}+\Gamma _{ij}^{2}\cdot \overrightarrow{e}_{2}+\Gamma _{ij}^{3}\cdot \overrightarrow{e% }_{3} \end{eqnarray*}

\begin{eqnarray*} \frac{\partial \overrightarrow{u}}{\partial \rho } &=&\left \{ \frac{\partial }{\partial \rho }\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \right \} \overrightarrow{e}_{\theta }+\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \left \{ \frac{\partial }{\partial \rho }% \overrightarrow{e}_{\theta }\right \} \rightarrow \Gamma _{21}^{k}\cdot \overrightarrow{e}_{k} \\ &=&\left \{ \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \right \} \overrightarrow{e}_{\theta }+\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \left \{ 0\right \} \\ &=&\left \{ \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \right \} \overrightarrow{e}_{\theta } \end{eqnarray*}

\begin{eqnarray*} \frac{\partial \overrightarrow{u}}{\partial \theta } &=&\left \{ \frac{% \partial }{\partial \theta }\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \right \} \overrightarrow{e}_{\theta }+\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \left \{ \frac{\partial }{\partial \theta }% \overrightarrow{e}_{\theta }\right \} \rightarrow \Gamma _{22}^{k}\cdot \overrightarrow{e}_{k} \\ &=&\left \{ e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \right \} \overrightarrow{e}_{\theta }+\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \left \{ -\overrightarrow{e}_{\rho }\right \} \\ &=&\left \{ e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \right \} \overrightarrow{e}_{\theta }-\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \overrightarrow{e}_{\rho } \end{eqnarray*}

\begin{eqnarray*} \frac{\partial \overrightarrow{u}}{\partial z} &=&0\overrightarrow{e}% _{\theta }+\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \left \{ \frac{\partial }{\partial z}\overrightarrow{e}_{\theta }\right \} \rightarrow \Gamma _{23}^{k}\cdot \overrightarrow{e}_{k} \\ &=&\frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \left \{ 0\right \} \\ &=&0 \end{eqnarray*}

\begin{eqnarray*} \nabla \overrightarrow{u}\left( \rho ,\theta \right) &=& \left( \begin{array}{ccc} 0 & -\frac{1}{2\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 0 \\ \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & \frac{1}{\rho }% e^{\rho -1}\cos \left( 2\theta -\pi /2\right) & 0 \\ 0 & 0 & 0% \end{array}% \right) \end{eqnarray*}

Una vez calculado el gradiente, ya podemos obtener el valor de [math]\epsilon (\vec{u}) [/math]:

\begin{eqnarray*} \epsilon \left( \overrightarrow{u}\right) &=&\left( \nabla \overrightarrow{u% }+\nabla \overrightarrow{u}^{t}\right) /2 \\ &=&\frac{1}{2}\left( \begin{array}{ccc} 0 & -\frac{1}{2\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 0 \\ \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & \frac{1}{\rho }% e^{\rho -1}\cos \left( 2\theta -\pi /2\right) & 0 \\ 0 & 0 & 0% \end{array}% \right) +\frac{1}{2}\left( \begin{array}{ccc} 0 & \frac{1}{2}e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 0 \\ -\frac{1}{2\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & \frac{1}{% \rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) & 0 \\ 0 & 0 & 0% \end{array}% \right) \\ &=&\left( \begin{array}{ccc} 0 & \frac{1}{4}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 0 \\ \frac{1}{4}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & \frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) & 0 \\ 0 & 0 & 0% \end{array}% \right) . \end{eqnarray*}

Ahora que ya tenemos [math]\epsilon (\vec{u}) [/math] calculado, podemos por fin calcular [math]\sigma (\vec{u}) [/math]:

\begin{eqnarray*} \sigma &=&\lambda \nabla \cdot \overrightarrow{u}\mathbf{1}+2\mu \epsilon \\ &=&\lambda \left( \frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \right) \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1% \end{array}% \right) +2\mu \left( \begin{array}{ccc} 0 & \frac{1}{4}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 0 \\ \frac{1}{4}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & \frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) & 0 \\ 0 & 0 & 0% \end{array}% \right) . \end{eqnarray*}

Para [math] \lambda =\mu =1 [/math] (coeficientes de Lamé, dependen de las propiedades elásticas de cada material):

\begin{eqnarray*} \sigma &=&\frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1% \end{array}% \right) +\left( \begin{array}{ccc} 0 & \frac{1}{2}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 0 \\ \frac{1}{2}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 2\frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) & 0 \\ 0 & 0 & 0% \end{array}% \right) \\ &=&\left( \begin{array}{ccc} \frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) & \frac{1}{2}% \frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 0 \\ \frac{1}{2}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) & 3\frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) & 0 \\ 0 & 0 & \frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \end{array}% \right) \end{eqnarray*}

Para acabar, las tensiones normales a cada dirección son:

\begin{eqnarray*} \overrightarrow{e_\rho}\cdot \sigma \cdot \overrightarrow{e_\rho} &=&\frac{1}{\rho }% e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \\ \overrightarrow{e_\theta}\cdot \sigma \cdot \overrightarrow{e_\theta} &=&3\frac{1}{\rho }% e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \\ \overrightarrow{e_z}\cdot \sigma \cdot \overrightarrow{e_z} &=&\frac{1}{\rho }% e^{\rho -1}\cos \left( 2\theta -\pi /2\right) . \end{eqnarray*}

A continuación, se adjunta la representación gráfica de dichas tensiones:

Figura11. Tensiones en cada uno de los ejes en tres dimensiones

Figura 12. Visualización de las tensiones

6 Tensiones tangenciales

En este apartado, se calcularán las tensiones tangenciales respecto al plano ortogonal a [math]\vec{u} [/math], en este caso el plano ortogonal a ρ. En la Figura 13 se dibujaran solo las que no sean nulas.

\begin{eqnarray*} \left\|\sigma \cdot \overrightarrow{e_\rho}-\left( \overrightarrow{e_\rho} \cdot \sigma \cdot \overrightarrow{e_\rho} \right) \overrightarrow{e_\rho} \right\| &=& \left| \left( \begin{array}{c} \frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \\ \frac{1}{2}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) \\ 0 \end{array} \right) - \left( \begin{array}{c} \frac{1}{\rho }e^{\rho -1}\cos \left( 2\theta -\pi /2\right) \\ 0 \\ 0 \end{array} \right) \right| \\ &=&|\frac{1}{2}\frac{\rho -1}{\rho }e^{\rho -1}\sin \left( 2\theta -\pi /2\right) | \\ &=&\frac{1}{2}\frac{\rho -1}{\rho }e^{\rho -1} | \sin \left( 2\theta -\pi /2\right) |. \end{eqnarray*}

Figura14. Tensión tangencial en tres dimensiones

Figura 13. Tensión tangencial

7 Tensión de Von Mises

Para la tensión de Von Mises seguiremos la siguiente fórmula:

• [math] \sigma_V= \sqrt \frac{(\sigma _1-\sigma _2)^2+(\sigma _2-\sigma _3)^2+(\sigma _3-\sigma _1)^2}{2} [/math]

Donde σ1, σ2 y σ3 (también conocidos como tensiones principales) son los autovalores de σ (tensor de tensiones). Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro).

En la Figura 15 se puede ver que σ1, σ2 y σ3 representan las tensiones principales como autovalores de σ, siendo estos

Dirección ρ: [math] σ_{1}=▽⋅u + 2ϵ_{11}= \frac{ 3e^{ρ-1} \cos(2θ-\frac{π}{2})}{ρ} [/math]

Dirección θ: [math] σ_{2}=▽⋅u +2ϵ_{22}= \frac{ 7e^{ρ-1} \cos(2θ-\frac{π}{2})}{ρ} [/math]

Dirección Z: [math] σ_{3}=▽⋅u +2ϵ_{33}= \frac{ 3e^{ρ-1} \cos(2θ-\frac{π}{2})}{ρ} [/math]

Figura 16.Visualización de la Tensión de Von Mises en el eje Z

Figura 15.Tensión de Von Mises

8 Códigos Matlab

En este apartado se podrá ver los códigos en Matlab que se han utilizado a lo largo del estudio para representar gráficamente los apartados propuestos.

8.1 Representación Gráfica de la Sección de anillo

%Apartado 1
%Con estos dos comandos, nos aseguramos de no tener ningún valor guardado en la memoria, ni representado en la interfaz gráfica
clear;
clc;

%Parametrizamos la superficie
%h es el paso de nestra malla
h=2/10;
%ro y tt, son los dos vectores, que usaremos para definir una malla en la que representar nuestra figura
ro=1:h:2;
tt=linspace(0.464,2.678,20);

%Con Meshgrid, creamos una malla, que ocupa un "área" definida por los vectores ro y tt
[U,V]=meshgrid(ro,tt);

%Es importante saber en que tipo de coordenadas trabajamos, por ello y para más comodidad, vamos a pasar de cilíndricas a cartesianas
X=U.*cos(V);
Y=U.*sin(V);



%Representamos el mallado de la placa en dos y tres dimensiones

%Nos abre una ventana, nueva, en este caso la número 1
figure(1);

%Nos permite abrir subventanas dentro de una misma ventana
subplot(1,2,1);

%Representa gráficamente la malla de la superficie definida por las matrices X, Y y Z, en este caso, Z es una matriz nula,
mesh(X,Y,0*X);
%Nota; Una solución más elegante a 0*X, es usar "zeros(size(X))", que da el mismo resultado

%Nos da una perspectiva predefinida, en este caso la vista superior de la pieza
view(2);

%Nos limita la longitud de los ejes, a X=[-3,3] e Y=[-1,3]
axis([-3,3,-1,3]);

%Le da titulo al gráfico
title('Mallado');

%Le da titulo a los ejes del gráfico
xlabel('Eje X');
ylabel('Eje Y');

%Aqui simplemente se repiten los comandos que ya he explicado con anterioridad, con el extra de zlabel que tambien da titulo al eje Z y view(3) que nos da una vista "ortogonal"
subplot(1,2,2);
mesh(X,Y,0*X);
view(3);
axis([-3,3,-1,3]);
title('Mallado');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');

8.2 Representación Gráfica de las Curvas de Nivel de la Temperatura

%Apartado 2
%AL igual que antes, con estos dos comandos, nos aseguramos de no tener ningún valor guardado en la memoria, ni representado en la interfaz gráfica, no volveré a comentarlo
clear;
clc;

%Volvemos a parametrizar la superficie, creando un mallado y pasando las coordenadas de cilindricas a cartesianas
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);

[U,V]=meshgrid(ro,tt);

X=U.*cos(V);
Y=U.*sin(V);
%-----------------------------------
%Hasta aquí nada nuevo


%Nuestra función de Temperatura viene dada por T(x,y)=cos((y-3)^2+x)
T=cos((Y-3).^2+X);
%Estamos guardando los puntos de la temperatura en una matriz T, que luego usaremos para graficarla
%Nota; Se puede apreciar el uso de un punto antes del exponencial, esto se debe a que se está elevando una matriz y no un número normal, hay que indicarlo a Matlab

%Representamos la temperatura y las curvas de nivel

figure(1)
subplot(1,2,1);
%Este nuevo comando, es muy similar a "mesh", pero nos representa una superficie, y no una malla de la misma
surf(X,Y,T);

%Esta linea, nos genera una leyenda del mapa de colores usado para colorear la superficie
colorbar

title('Temperatura en placa');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
axis ([-3,3,-1,3]);

view(2);
subplot(1,2,2);

%Contour sirve para representar las curvas de nivel, en este caso, el "30" determina la densidad de dichas lineas
contour(X,Y,T,30);

colorbar
title('Curvas de nivel de la Temperatura');
xlabel('Eje X');
ylabel('Eje Y');
axis ([-3,3,-1,3]);

%El comando max busca el valor más alto de un vector, en este caso, al ser T una matriz, hay que ejecutarlo 2 veces seguidas
Tmax=max(max(T));

%Con este comando, mostramos en la consola el valor hallado
fprintf('La temperatura máxima que se alcanza es de: %1.5f ºC\n',Tmax);

8.2.1 Representación Gráfica del Gradiente

%Apartado 2, 2ªparte

clear;clc;
%Parametrizamos la superficie
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);
%Creamos el mallado y pasamos a coordenas cilindricas
[U,V]=meshgrid(ro,tt);
X=U.*cos(V);
Y=U.*sin(V);
%-----------------------------------
%Hasta aquí nada nuevo

%Volvemos a definir T
T=cos((Y-3).^2+X);

%Definimos las dos componentes del gradiente de T(X,Y), para luego poder representar el campo escalar
Tx=-sin((Y-3).^2+X)*1;
Ty=-sin((Y-3).^2+X).*(2*Y-6);

%A partir de aquí solo definimos como queremos que quede el gráfico
figure(1)
mesh (X,Y,0*X);
contour(X,Y,T);
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar

%Hold on, sirve para que Matlab no borre los datos de la ventana, si no que vuelva escribir sobre ella
hold on

%Quiver representa el campo escalar dado por el gradiente de T
quiver(X,Y,Tx,Ty);

axis equal
axis([-3,3,-1,3]);
view(2);
title('Gradiente de la temperatura');
xlabel('Eje X');
ylabel('Eje Y');

%Hold off, sirve para cerrar el hold on, no es técnicamente necesario en este caso, pero denota una falta de elegancia
hold off

8.3 Representación Gráfica de la Ley de Fourier

%Apartado 3
clear;clc;
%Parametrización de la superficie
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);
%Creamos el mallado y pasamos a coordenas cilindricas
[U,V]=meshgrid(ro,tt);
X=U.*cos(V);
Y=U.*sin(V);

%-----------------------------------
%Hasta aquí nada nuevo

%Definimos la Ley de Fourier donde la energia calorifica es Q=-k∇T

%Al igual que antes, definifmos la Temperatura y su gradiente, solo que esta vez, hemos multiplicado Tx y Ty por -k, siendo k=1
T=cos((Y-3).^2+X);
Tx=sin((Y-3).^2+X)*1;
Ty=sin((Y-3).^2+X).*(2*Y-6);


%Representamos el gradiente del campo vactorial
%A partir de aquí solo definimos como queremos que quede el gráfico
figure(1)
mesh (X,Y,0*X);
contour(X,Y,T);
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar
hold on
quiver(X,Y,Tx,Ty);
axis equal
axis([-3,3,-1,3]);
view(2);
title('Ley de Fourier');
xlabel('Eje X');
ylabel('Eje Y');
hold off

8.4 Representación Gráfica del Campo vectorial \begin{eqnarray*} \bar u \end{eqnarray*}

%Apartado 5
clear;
clc;
%Parametrización de la superficie
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);
%Creamos el mallado y pasamos a coordenas cilindricas
[U,V]=meshgrid(ro,tt);
X=U.*cos(V);
Y=U.*sin(V);
%Función
FX=-1/2.*exp(U-1).*sin(2.*V-pi/2).*sin(V);
FY=1/2.*exp(U-1).*sin(2.*V-pi/2).*cos(V);
%Representación gráfica
%Antes de la deformación
subplot(1,3,1)
surf(X,Y,0*X);
view(2)
axis equal
axis([-3,3,-1,3])
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
title('Antes de la deformación');
%Después de la deformación
subplot(1,3,2)
UX=X+FX;
UY=Y+FY;
surf(UX,UY,0*UX);
view(2)
axis equal
axis([-3,3,-1,3])
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
title('Después de la deformación');
 %Comparación
subplot(1,3,3)
hold on
plot(X,Y,'r');
plot(UX,UY,'b');
view(2)
axis equal
axis([-3,3,-1,3])
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
title('Comparación');
hold off

8.4.1 Representación Gráfica de la Divergencia

%Apartado 6
clear;
clc;
%Parametrización de la superficie
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);
%Creamos el mallado y pasamos a coordenas cilindricas
[U,V]=meshgrid(ro,tt);
X=U.*cos(V);
Y=U.*sin(V);
%-----------------------------------
%Hasta aquí nada nuevo


%Calculamos los valores de la divergencia, en t =0
Div=(exp(U-1)).*cos(2.*V-(pi/2)))./U;


%Creamos la gráfica
surf(X,Y,Div);
view(2);
axis equal
axis([-3,3,-1,3]);
title('Divergencia del campo vectorial')
colorbar
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');

%Encontramos la divergencia máxima y mínima
MaxDiv=max(Div);
MinDiv=min(Div);
%El comando find, nos permite hallar los valores que busquemos en una matriz
NullDiv= find(Div == 0);


fprintf('Los puntos máximos de la divergencia son: %1.5f \n',MaxDiv);
fprintf('Los puntos mínimos de la divergencia son: %1.5f \n',MinDiv);
fprintf('Los puntos nulos de la divergencia son: %1.5f \n',NullDiv);

8.4.2 Representación Gráfica del Rotacional

%Apartado 7
clear;
clc;

%Parametrizamos la superficie
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);
%Creamos el mallado y pasamos a coordenas cilindricas
[U,V]=meshgrid(ro,tt);
X=U.*cos(V);
Y=U.*sin(V);
%-----------------------------------
%Hasta aquí nada nuevo

%Calculamos el valor de el rotacional
Rot=abs(((1+U)./2.*U).*(exp(U-1).*cos(2.*V-(pi/2)))./(2.*U));

%Creamos la gráfica
surf(X,Y,Rot);
view(2);
axis equal
axis([-3,3,-1,3]);
title('Rotacional del campo vectorial')
colorbar
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');

%Encontramos los puntos que sufren mayor rotación
MaxRot=max(Rot);
fprintf('Los puntos que sufren mayor rotación son: %1.5f \n',MaxRot);

8.5 Representación Gráfica del Tensor de tensiones

%Aparatdo 8
clear;
clc;
%Parametrizamos la superficie
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);
%Creamos el mallado y pasamos a coordenas cilindricas
[U,V]=meshgrid(ro,tt);
X=U.*cos(V);
Y=U.*sin(V);
%-----------------------------------
%Hasta aquí nada nuevo


%Tensiones normales en la dirección del eje Ro, para ello calculamos los valores de Sigma1
Sigma1=(1./(U)).*(exp(U-1)).*cos((2.*V)-(pi/2)); 
subplot(1,3,1);
surf(X,Y,Sigma1);
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');
ylabel('Eje Y');
colorbar;
view(2);
title('Tensiones normales en la dirección del eje Ro')


%Tensiones normales en la dirección del eje Tt, para ello calculamos los valores de Sigma2
Sigma2=(3./(U)).*(exp(U-1)).*cos((2.*V)-(pi/2));
subplot(1,3,2);
surf(X,Y,Sigma2);
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');
ylabel('Eje Y');
colorbar;
view(2);
title('Tensiones normales en la dirección del eje Tt')


%Tensiones normales en la dirección del eje Z, para ello calculamos los valores de Sigma3
Sigma3=(1./(U)).*(exp(U-1)).*cos((2.*V)-(pi/2));
subplot(1,3,3);
surf(X,Y,Sigma3);
axis equal
axis([-3,3,-1,3]);
xlabel('Eje X');
ylabel('Eje Y');
colorbar
view(2)
title('Tensiones normales en la dirección del eje Z');

8.6 Representación Gráfica de las Tensiones tangenciales

%Apartado 9
clear;clc;
%Parametrizamos la superficie
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);
%Creamos el mallado y pasamos a coordenas cilindricas
[U,V]=meshgrid(ro,tt);
X=U.*cos(V);
Y=U.*sin(V);
%Tensiones tangenciales respecto al plano ortogonal al eje Ro
Sigma=(1/2).*((U-1)./U).*(exp(U-1)).*(sin((2.*V)-(pi/2))).^2;
surf(X,Y,Sigma);
axis equal
axis([-3,3,-1,3]);
colorbar
view(2)
xlabel('Eje X');
ylabel('Eje Y');
title('Tensiones tangenciales');

8.7 Representación Gráfica de la Tensión de Von Mises

%Apartado 10
clear;clc;
%Parametrizamos la superficie
h=2/10;
ro=1:h:2;
tt=linspace(0.464,2.678,20);
%Creamos el mallado y pasamos a coordenas cilindricas
[U,V]=meshgrid(ro,tt);
X=U.*cos(V);
Y=U.*sin(V);
Z=zeros(size(X));
clf
%Procedemos a sacar los autovalores
Sig1=(3.*exp(U-1)).*cos(2.*V-(pi/2))./U; 
Sig2=(7.*exp(U-1)).*cos(2.*V-(pi/2))./U; 
Sig3=(3.*exp(U-1)).*cos(2.*V-(pi/2))./U; 
%Calculamos la tension de Von Mises usando la formula.
SigVM=sqrt((Sig1-Sig2).^2 +(Sig2-Sig3).^2 +(Sig3-Sig1).^2/2);
%Creamos la grafica de la superficie
surf(X,Y,SigVM)
axis([-3,3,-1,3]);
view(2)
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
title('Tension de Von Mises')