MateWiki - Contribuciones del usuario [es]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-11T09:11:21Z

Julio Montero Domínguez: /* Divergencia del campo de desplazamientos */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
[[Archivo:1.jpeg|right|300px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).

En la gráfica podemos observar que la temperatura máxima la encontramos en la parte superior de nuestro anillo, en el punto (0,2) y la temperatura mínima en la base de nuestro semianillo.

[[Archivo:2.jpeg|right|400px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
[[Archivo:Grupo19A 2020-21 Gradiente temperatura.jpeg|right|]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}

==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
La divergencia la usamos para observar cómo cambia el volumen local de las distintas áreas de la placa cuando se aplica una fuerza que ejerce un desplazamiento, en este caso, el desplazamiento es el vector <math>\vec u</math>.

Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \rho \vec u_\rho}{\partial \rho} +\frac{1}{\rho} \frac{\partial \rho \vec u_\theta}{\partial \theta} +\frac{1}{\rho} \frac{\partial \rho \vec u_z}{\partial z}
=\frac{1}{\rho}0 + \frac{1}{\rho} \rho (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2}) cos(\theta) + 0
</math>

Y realizando estos cálculos obtenemos el resultado final para la divergencia:

<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>

A continuación obtendremos analíticamente los puntos para los que la divergencia es nula, máxima y mínima, aunque esto ya lo podemos observar en la representación gráfica que hemos hecho con matlab.

La divergencia es nula o bien cuando el <math>\ cos(\theta)= 0</math> (es decir para los valores de <math>\ \theta</math> :<math>\ \theta=\frac{pi}{2}</math> y <math>\ \theta=\frac{3pi}{2}</math>) o si <math>\ (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})=0</math> es decir para el valor <math>\ \rho=1</math>

La divergencia es máxima para <math>\ \theta=2pi</math> ,<math>\ \theta=0</math> (que es cuando el coseno es máximo) y además para los <math>\ \rho</math> que cumplan la siguiente ecuación:
<math> 34\rho ^3-3\rho ^2-10=0</math> Resultado de máximizar (derivando e igualando a cero):<math>\ (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>

La divergencia mínima se obtiene para el valor de <math>\ \theta=pi</math> y al igual que para la divergencia máxima para los valores de <math>\ \rho</math> que cumplan la ecuación:
<math> 34\rho ^3-3\rho ^2-10=0</math>

En la gráfica se puede observar como en la parte derecha de la placa la divergencia es máxima, es decir, el volumen local se expande. En cambio, en la parte izquierda del semianillo la divergencia es negativa, por lo que el volumen local se contrae.
[[Archivo:divergencia7.png|right|550px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==
El rotacional se utiliza para ver cómo afecta el campo de fuerzas a un objeto en su movimiento, pero, en especial, a su movimiento angular, es decir, el rotacional nos indica cómo girará un cuerpo en un determinado punto afectado por el campo de fuerzas.

Como hemos hecho con anterioridad, calcularemos el rotacional del campo como se expone abajo para representarlo posteriormente en un gráfico con indicaciones a calor.

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))
}}

[[Archivo: Captura rotacional.PNG|1000px]]

Como se puede apreciar en la imagen, existe un mayor rotacional en la parte superior central. En cambio, en la parte inferior el rotacional es próximo a 0, es decir, la placa casi no rota sobre si misma.

==Tensión de deformaciones==
Nos definen el tensor :<math> e \vec u </math> como <math> e \vec u = \frac{(\nabla \vec u + \nabla \vec u ^t)}{2} </math>

Además nos permiten escribir el tensor de tensiones a través de la fórmula: <math> \sigma= ẞ \nabla \cdot \vec u I +2ẟe</math>

Primero calcularemos <math> \nabla \vec u </math>= ( <math> senθ(\frac{2}{15}-\frac{1}{20\rho}+\frac{1}{20\rho3}\vec g_θ</math>,<math>-\rho senθ \frac{4\rho-3}{60}+\frac{1}{60\rho^2}\vec g_\rho </math>,<math>cosθ(\frac{4\rho -3}{60}+\frac{1}{60\rho^2})\vec g_θ </math>)

Para simplificar la notación hemos llamado <math>b</math> a <math>\frac{4\rho-3}{60}+\frac{1}{60\rho^2}</math> y hemos llamado <math>a</math> a <math> senθ (\frac{2}{15}-\frac{1}{20\rho}+\frac{1}{20\rho3})</math>

Lo primero que haremos sera calcular e para luego poder sustituirlo en la fórmula que nos da las tensiones, para ello necesitaremos el gradiente y su traspuesto:
<math>\nabla \vec u= \begin{pmatrix}0&-\rho senθb&0\\a&cosθ b&0\\0&0&0\end{pmatrix} </math> <math>\nabla \vec u ^t= \begin{pmatrix}0& a& 0\\-\rho senθb & cosθ b & 0 \\0&0&0 \end{pmatrix} </math>

Una vez hallado esto tenemos que:

<math> e(u) = \frac{\begin{pmatrix} 0 & a-\rho senθ b & 0 \\ a-\rho senθ b & 2cosθ b & 0 \\ 0 & 0 & 0 \end{pmatrix}}{2} </math>

Ya podemos sustituir en la fórmula de la tensíon para unos coeficientes de lamé ẞ=ẟ=1 :

<math> \sigma= \begin{pmatrix} cosθ b & 0 & 0 \\ 0 & cosθ b & 0 \\ 0 & 0 & cosθ b \end{pmatrix} + \begin{pmatrix} 0 & a-\rho senθ b & 0 \\ a-\rho senθ & 2cosθ b & 0 \\ 0 & 0 & 0 \end{pmatrix}= \begin{pmatrix} cosθ b & a-\rho senθ b & 0 \\ a-\rho senθ b & 3cosθ b & 0 \\ 0 & 0 & cosθ b \end{pmatrix} </math>
Operando obtenemos que:

<math> \sigma= \begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix} </math>

[[Archivo:demasiadastensiones.jpeg|centre|1200px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);

%Trazado de la tensión normal grho*tensión*grho
t1=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta);
subplot(1,3,1)
surf(MX,MY,t1)
title('grho*tensión*grho')

%Trazado de la tensiónn normal (grho/rho)tensión(grho/rho)
t2=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta).*(1./mrho.^2);
subplot(1,3,2)
surf(MX,MY,t2)
axis equal
title('(grho/rho).*tensión(grho/rho)')

%Trazado de la tensiónn normal k*tensión*k
t3=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta);
subplot(1,3,3)
surf(MX,MY,t3)
axis equal
title('k*tensión*k')
}}

==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==

El círculo de Mohr es un método gráfico que nos ayuda a determinar el estado tensional en los distintos puntos de un cuerpo.

Gracias a este círculo podemos saber cuáles son las tensiones principales de nuestro problema.

[[Archivo:5.png|right]]
{{matlab|codigo=
%círculo de mohr
sigx=75;
sigy=45;
sigz=50;
tauxy=50;
R=sqrt(((sigx-sigy)/2)^2+(tauxy)^2)
C=(sigx+sigy)/2;
sig1=C+R
sig2=C-R
sigma=linspace(sig2,sig1,1001);
tau1=sqrt(R^2-(sigma-C).^2);
tau2=-tau1
plot(sigma,tau1,'k',sigma,tau2,'k')
xlabel('Sigma')
ylabel('tau')
axis('square')
hold on
plot([sig1,sig2],[0,0],'k')
plot([sigx,sigy],[-tauxy,tauxy],'r')
grid on
x1= sig2:0.15:sig1+0.25
y1= (sigz/(R-C))*x1-(sigz/(R-C))*C
x2=sig2:0.15:sig1
y2=0
x3=-C
y3=-R:0.15:R
plot(X,Y,x1,y1,x2,y2,x3,y3)
hold off
}}

==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-11T09:08:40Z

Julio Montero Domínguez: /* Divergencia del campo de desplazamientos */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
[[Archivo:1.jpeg|right|300px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).

En la gráfica podemos observar que la temperatura máxima la encontramos en la parte superior de nuestro anillo, en el punto (0,2) y la temperatura mínima en la base de nuestro semianillo.

[[Archivo:2.jpeg|right|400px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
[[Archivo:Grupo19A 2020-21 Gradiente temperatura.jpeg|right|]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}

==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
La divergencia la usamos para observar cómo cambia el volumen local de las distintas áreas de la placa cuando se aplica una fuerza que ejerce un desplazamiento, en este caso, el desplazamiento es el vector <math>\vec u</math>.

Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \rho \vec u_\rho}{\partial \rho} +\frac{1}{\rho} \frac{\partial \rho \vec u_\theta}{\partial \theta} +\frac{1}{\rho} \frac{\partial \rho \vec u_z}{\partial z}
=\frac{1}{\rho}0 + \frac{1}{\rho} \rho (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2}) cos(\theta) + 0
</math>

Y realizando estos cálculos obtenemos el resultado final para la divergencia:

<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>

A continuación obtendremos analíticamente los puntos para los que la divergencia es nula, máxima y mínima, aunque esto ya lo podemos observar en la representación gráfica que hemos hecho con matlab.

La divergencia es nula o bien cuando el <math>\ cos(\theta)= 0</math> (es decir para los valores de <math>\ \theta</math> :<math>\ \theta=\frac{pi}{2}</math> y <math>\ \theta=\frac{3pi}{2}</math>) o si <math>\ (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})=0</math> es decir para el valor <math>\ \rho=1</math>

La divergencia es máxima para <math>\ \theta=2pi</math> ,<math>\ \theta=0</math> (que es cuando el coseno es máximo) y además para los <math>\ \rho</math> que cumplan la siguiente ecuación:
<math> 34\rho ^3-3\rho ^2-10=0</math> Resultado de máximizar (derivando e igualando a cero):<math>\ (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>

La divergencia mínima se obtiene para el valor de <math>\ \theta=pi</math> y al igual que para la divergencia máxima para los valores de <math>\ \rho</math> que cumplan la ecuación:
<math> 34\rho ^3-3\rho ^2-10=0</math>

En la gráfica podemos ver como la placa se deforma en cada punto va cambiando debido al desplazamiento. Este es el sentido que tiene la divergencia.
[[Archivo:divergencia7.png|right|550px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==
El rotacional se utiliza para ver cómo afecta el campo de fuerzas a un objeto en su movimiento, pero, en especial, a su movimiento angular, es decir, el rotacional nos indica cómo girará un cuerpo en un determinado punto afectado por el campo de fuerzas.

Como hemos hecho con anterioridad, calcularemos el rotacional del campo como se expone abajo para representarlo posteriormente en un gráfico con indicaciones a calor.

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))
}}

[[Archivo: Captura rotacional.PNG|1000px]]

Como se puede apreciar en la imagen, existe un mayor rotacional en la parte superior central. En cambio, en la parte inferior el rotacional es próximo a 0, es decir, la placa casi no rota sobre si misma.

==Tensión de deformaciones==
Nos definen el tensor :<math> e \vec u </math> como <math> e \vec u = \frac{(\nabla \vec u + \nabla \vec u ^t)}{2} </math>

Además nos permiten escribir el tensor de tensiones a través de la fórmula: <math> \sigma= ẞ \nabla \cdot \vec u I +2ẟe</math>

Primero calcularemos <math> \nabla \vec u </math>= ( <math> senθ(\frac{2}{15}-\frac{1}{20\rho}+\frac{1}{20\rho3}\vec g_θ</math>,<math>-\rho senθ \frac{4\rho-3}{60}+\frac{1}{60\rho^2}\vec g_\rho </math>,<math>cosθ(\frac{4\rho -3}{60}+\frac{1}{60\rho^2})\vec g_θ </math>)

Para simplificar la notación hemos llamado <math>b</math> a <math>\frac{4\rho-3}{60}+\frac{1}{60\rho^2}</math> y hemos llamado <math>a</math> a <math> senθ (\frac{2}{15}-\frac{1}{20\rho}+\frac{1}{20\rho3})</math>

Lo primero que haremos sera calcular e para luego poder sustituirlo en la fórmula que nos da las tensiones, para ello necesitaremos el gradiente y su traspuesto:
<math>\nabla \vec u= \begin{pmatrix}0&-\rho senθb&0\\a&cosθ b&0\\0&0&0\end{pmatrix} </math> <math>\nabla \vec u ^t= \begin{pmatrix}0& a& 0\\-\rho senθb & cosθ b & 0 \\0&0&0 \end{pmatrix} </math>

Una vez hallado esto tenemos que:

<math> e(u) = \frac{\begin{pmatrix} 0 & a-\rho senθ b & 0 \\ a-\rho senθ b & 2cosθ b & 0 \\ 0 & 0 & 0 \end{pmatrix}}{2} </math>

Ya podemos sustituir en la fórmula de la tensíon para unos coeficientes de lamé ẞ=ẟ=1 :

<math> \sigma= \begin{pmatrix} cosθ b & 0 & 0 \\ 0 & cosθ b & 0 \\ 0 & 0 & cosθ b \end{pmatrix} + \begin{pmatrix} 0 & a-\rho senθ b & 0 \\ a-\rho senθ & 2cosθ b & 0 \\ 0 & 0 & 0 \end{pmatrix}= \begin{pmatrix} cosθ b & a-\rho senθ b & 0 \\ a-\rho senθ b & 3cosθ b & 0 \\ 0 & 0 & cosθ b \end{pmatrix} </math>
Operando obtenemos que:

<math> \sigma= \begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix} </math>

[[Archivo:demasiadastensiones.jpeg|centre|1200px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);

%Trazado de la tensión normal grho*tensión*grho
t1=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta);
subplot(1,3,1)
surf(MX,MY,t1)
title('grho*tensión*grho')

%Trazado de la tensiónn normal (grho/rho)tensión(grho/rho)
t2=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta).*(1./mrho.^2);
subplot(1,3,2)
surf(MX,MY,t2)
axis equal
title('(grho/rho).*tensión(grho/rho)')

%Trazado de la tensiónn normal k*tensión*k
t3=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta);
subplot(1,3,3)
surf(MX,MY,t3)
axis equal
title('k*tensión*k')
}}

==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==

El círculo de Mohr es un método gráfico que nos ayuda a determinar el estado tensional en los distintos puntos de un cuerpo.

Gracias a este círculo podemos saber cuáles son las tensiones principales de nuestro problema.

[[Archivo:5.png|right]]
{{matlab|codigo=
%círculo de mohr
sigx=75;
sigy=45;
sigz=50;
tauxy=50;
R=sqrt(((sigx-sigy)/2)^2+(tauxy)^2)
C=(sigx+sigy)/2;
sig1=C+R
sig2=C-R
sigma=linspace(sig2,sig1,1001);
tau1=sqrt(R^2-(sigma-C).^2);
tau2=-tau1
plot(sigma,tau1,'k',sigma,tau2,'k')
xlabel('Sigma')
ylabel('tau')
axis('square')
hold on
plot([sig1,sig2],[0,0],'k')
plot([sigx,sigy],[-tauxy,tauxy],'r')
grid on
x1= sig2:0.15:sig1+0.25
y1= (sigz/(R-C))*x1-(sigz/(R-C))*C
x2=sig2:0.15:sig1
y2=0
x3=-C
y3=-R:0.15:R
plot(X,Y,x1,y1,x2,y2,x3,y3)
hold off
}}

==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T23:34:06Z

Julio Montero Domínguez: /* Círculo de Mohr */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
[[Archivo:1.jpeg|right|300px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).

En la gráfica podemos observar que la temperatura máxima la encontramos en la parte superior de nuestro anillo, en el punto (0,2) y la temperatura mínima en la base de nuestro semianillo.

[[Archivo:2.jpeg|right|400px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
[[Archivo:Grupo19A 2020-21 Gradiente temperatura.jpeg|right|]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}

==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
La divergencia la usamos para observar cómo de rápido cierto punto varía posición con respecto a una fuerza que ejerce un desplazamiento, en este caso, el desplazamiento es el vector <math>\vec u</math>.

Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \rho \vec u_\rho}{\partial \rho} +\frac{1}{\rho} \frac{\partial \rho \vec u_\theta}{\partial \theta} +\frac{1}{\rho} \frac{\partial \rho \vec u_z}{\partial z}
=\frac{1}{\rho}0 + \frac{1}{\rho} \rho (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2}) cos(\theta) + 0
</math>

Y realizando estos cálculos obtenemos el resultado final para la divergencia:

<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>

A continuación obtendremos analíticamente los puntos para los que la divergencia es nula, máxima y mínima, aunque esto ya lo podemos observar en la representación gráfica que hemos hecho con matlab.

La divergencia es nula o bien cuando el <math>\ cos(\theta)= 0</math> (es decir para los valores de <math>\ \theta</math> :<math>\ \theta=\frac{pi}{2}</math> y <math>\ \theta=\frac{3pi}{2}</math>) o si <math>\ (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})=0</math> es decir para el valor <math>\ \rho=1</math>

La divergencia es máxima para <math>\ \theta=2pi</math> ,<math>\ \theta=0</math> (que es cuando el coseno es máximo) y además para los <math>\ \rho</math> que cumplan la siguiente ecuación:
<math> 34\rho ^3-3\rho ^2-10=0</math> Resultado de máximizar (derivando e igualando a cero):<math>\ (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>

La divergencia mínima se obtiene para el valor de <math>\ \theta=pi</math> y al igual que para la divergencia máxima para los valores de <math>\ \rho</math> que cumplan la ecuación:
<math> 34\rho ^3-3\rho ^2-10=0</math>

En la gráfica podemos ver como la placa se deforma en cada punto va cambiando debido al desplazamiento. Este es el sentido que tiene la divergencia.
[[Archivo:divergencia7.png|right|550px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==
El rotacional se utiliza para ver cómo afecta el campo de fuerzas a un objeto en su movimiento, pero, en especial, a su movimiento angular, es decir, el rotacional nos indica cómo girará un cuerpo en un determinado punto afectado por el campo de fuerzas.

Como hemos hecho con anterioridad, calcularemos el rotacional del campo como se expone abajo para representarlo posteriormente en un gráfico con indicaciones a calor.

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))
}}

[[Archivo: Captura rotacional.PNG|1000px]]

Como se puede apreciar en la imagen, existe un mayor rotacional en la parte superior central. En cambio, en la parte inferior el rotacional es próximo a 0, es decir, la placa casi no rota sobre si misma.

==Tensión de deformaciones==
Nos definen el tensor :<math> e \vec u </math> como <math> e \vec u = \frac{(\nabla \vec u + \nabla \vec u ^t)}{2} </math>

Además nos permiten escribir el tensor de tensiones a través de la fórmula: <math> \sigma= ẞ \nabla \cdot \vec u I +2ẟe</math>

Primero calcularemos <math> \nabla \vec u </math>= ( <math> senθ(\frac{2}{15}-\frac{1}{20\rho}+\frac{1}{20\rho3}\vec g_θ</math>,<math>-\rho senθ \frac{4\rho-3}{60}+\frac{1}{60\rho^2}\vec g_\rho </math>,<math>cosθ(\frac{4\rho -3}{60}+\frac{1}{60\rho^2})\vec g_θ </math>)

Para simplificar la notación hemos llamado <math>b</math> a <math>\frac{4\rho-3}{60}+\frac{1}{60\rho^2}</math> y hemos llamado <math>a</math> a <math> senθ (\frac{2}{15}-\frac{1}{20\rho}+\frac{1}{20\rho3})</math>

Lo primero que haremos sera calcular e para luego poder sustituirlo en la fórmula que nos da las tensiones, para ello necesitaremos el gradiente y su traspuesto:
<math>\nabla \vec u= \begin{pmatrix}0&-\rho senθb&0\\a&cosθ b&0\\0&0&0\end{pmatrix} </math> <math>\nabla \vec u ^t= \begin{pmatrix}0& a& 0\\-\rho senθb & cosθ b & 0 \\0&0&0 \end{pmatrix} </math>

Una vez hallado esto tenemos que:

<math> e(u) = \frac{\begin{pmatrix} 0 & a-\rho senθ b & 0 \\ a-\rho senθ b & 2cosθ b & 0 \\ 0 & 0 & 0 \end{pmatrix}}{2} </math>

Ya podemos sustituir en la fórmula de la tensíon para unos coeficientes de lamé ẞ=ẟ=1 :

<math> \sigma= \begin{pmatrix} cosθ b & 0 & 0 \\ 0 & cosθ b & 0 \\ 0 & 0 & cosθ b \end{pmatrix} + \begin{pmatrix} 0 & a-\rho senθ b & 0 \\ a-\rho senθ & 2cosθ b & 0 \\ 0 & 0 & 0 \end{pmatrix}= \begin{pmatrix} cosθ b & a-\rho senθ b & 0 \\ a-\rho senθ b & 3cosθ b & 0 \\ 0 & 0 & cosθ b \end{pmatrix} </math>
Operando obtenemos que:

<math> \sigma= \begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix} </math>

[[Archivo:demasiadastensiones.jpeg|centre|1200px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);

%Trazado de la tensión normal grho*tensión*grho
t1=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta);
subplot(1,3,1)
surf(MX,MY,t1)
title('grho*tensión*grho')

%Trazado de la tensiónn normal (grho/rho)tensión(grho/rho)
t2=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta).*(1./mrho.^2);
subplot(1,3,2)
surf(MX,MY,t2)
axis equal
title('(grho/rho).*tensión(grho/rho)')

%Trazado de la tensiónn normal k*tensión*k
t3=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta);
subplot(1,3,3)
surf(MX,MY,t3)
axis equal
title('k*tensión*k')
}}

==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==

El círculo de Mohr es un método gráfico que nos ayuda a determinar el estado tensional en los distintos puntos de un cuerpo.

Gracias a este círculo podemos saber cuáles son las tensiones principales de nuestro problema.

[[Archivo:5.png|right]]
{{matlab|codigo=
%círculo de mohr
sigx=75;
sigy=45;
sigz=50;
tauxy=50;
R=sqrt(((sigx-sigy)/2)^2+(tauxy)^2)
C=(sigx+sigy)/2;
sig1=C+R
sig2=C-R
sigma=linspace(sig2,sig1,1001);
tau1=sqrt(R^2-(sigma-C).^2);
tau2=-tau1
plot(sigma,tau1,'k',sigma,tau2,'k')
xlabel('Sigma')
ylabel('tau')
axis('square')
hold on
plot([sig1,sig2],[0,0],'k')
plot([sigx,sigy],[-tauxy,tauxy],'r')
grid on
x1= sig2:0.15:sig1+0.25
y1= (sigz/(R-C))*x1-(sigz/(R-C))*C
x2=sig2:0.15:sig1
y2=0
x3=-C
y3=-R:0.15:R
plot(X,Y,x1,y1,x2,y2,x3,y3)
hold off
}}

==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T23:33:52Z

Julio Montero Domínguez: /* Círculo de Mohr */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
[[Archivo:1.jpeg|right|300px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).

En la gráfica podemos observar que la temperatura máxima la encontramos en la parte superior de nuestro anillo, en el punto (0,2) y la temperatura mínima en la base de nuestro semianillo.

[[Archivo:2.jpeg|right|400px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
[[Archivo:Grupo19A 2020-21 Gradiente temperatura.jpeg|right|]]
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}

==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
La divergencia la usamos para observar cómo de rápido cierto punto varía posición con respecto a una fuerza que ejerce un desplazamiento, en este caso, el desplazamiento es el vector <math>\vec u</math>.

Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \rho \vec u_\rho}{\partial \rho} +\frac{1}{\rho} \frac{\partial \rho \vec u_\theta}{\partial \theta} +\frac{1}{\rho} \frac{\partial \rho \vec u_z}{\partial z}
=\frac{1}{\rho}0 + \frac{1}{\rho} \rho (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2}) cos(\theta) + 0
</math>

Y realizando estos cálculos obtenemos el resultado final para la divergencia:

<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>

A continuación obtendremos analíticamente los puntos para los que la divergencia es nula, máxima y mínima, aunque esto ya lo podemos observar en la representación gráfica que hemos hecho con matlab.

La divergencia es nula o bien cuando el <math>\ cos(\theta)= 0</math> (es decir para los valores de <math>\ \theta</math> :<math>\ \theta=\frac{pi}{2}</math> y <math>\ \theta=\frac{3pi}{2}</math>) o si <math>\ (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})=0</math> es decir para el valor <math>\ \rho=1</math>

La divergencia es máxima para <math>\ \theta=2pi</math> ,<math>\ \theta=0</math> (que es cuando el coseno es máximo) y además para los <math>\ \rho</math> que cumplan la siguiente ecuación:
<math> 34\rho ^3-3\rho ^2-10=0</math> Resultado de máximizar (derivando e igualando a cero):<math>\ (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>

La divergencia mínima se obtiene para el valor de <math>\ \theta=pi</math> y al igual que para la divergencia máxima para los valores de <math>\ \rho</math> que cumplan la ecuación:
<math> 34\rho ^3-3\rho ^2-10=0</math>

En la gráfica podemos ver como la placa se deforma en cada punto va cambiando debido al desplazamiento. Este es el sentido que tiene la divergencia.
[[Archivo:divergencia7.png|right|550px]]
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==
El rotacional se utiliza para ver cómo afecta el campo de fuerzas a un objeto en su movimiento, pero, en especial, a su movimiento angular, es decir, el rotacional nos indica cómo girará un cuerpo en un determinado punto afectado por el campo de fuerzas.

Como hemos hecho con anterioridad, calcularemos el rotacional del campo como se expone abajo para representarlo posteriormente en un gráfico con indicaciones a calor.

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))
}}

[[Archivo: Captura rotacional.PNG|1000px]]

Como se puede apreciar en la imagen, existe un mayor rotacional en la parte superior central. En cambio, en la parte inferior el rotacional es próximo a 0, es decir, la placa casi no rota sobre si misma.

==Tensión de deformaciones==
Nos definen el tensor :<math> e \vec u </math> como <math> e \vec u = \frac{(\nabla \vec u + \nabla \vec u ^t)}{2} </math>

Además nos permiten escribir el tensor de tensiones a través de la fórmula: <math> \sigma= ẞ \nabla \cdot \vec u I +2ẟe</math>

Primero calcularemos <math> \nabla \vec u </math>= ( <math> senθ(\frac{2}{15}-\frac{1}{20\rho}+\frac{1}{20\rho3}\vec g_θ</math>,<math>-\rho senθ \frac{4\rho-3}{60}+\frac{1}{60\rho^2}\vec g_\rho </math>,<math>cosθ(\frac{4\rho -3}{60}+\frac{1}{60\rho^2})\vec g_θ </math>)

Para simplificar la notación hemos llamado <math>b</math> a <math>\frac{4\rho-3}{60}+\frac{1}{60\rho^2}</math> y hemos llamado <math>a</math> a <math> senθ (\frac{2}{15}-\frac{1}{20\rho}+\frac{1}{20\rho3})</math>

Lo primero que haremos sera calcular e para luego poder sustituirlo en la fórmula que nos da las tensiones, para ello necesitaremos el gradiente y su traspuesto:
<math>\nabla \vec u= \begin{pmatrix}0&-\rho senθb&0\\a&cosθ b&0\\0&0&0\end{pmatrix} </math> <math>\nabla \vec u ^t= \begin{pmatrix}0& a& 0\\-\rho senθb & cosθ b & 0 \\0&0&0 \end{pmatrix} </math>

Una vez hallado esto tenemos que:

<math> e(u) = \frac{\begin{pmatrix} 0 & a-\rho senθ b & 0 \\ a-\rho senθ b & 2cosθ b & 0 \\ 0 & 0 & 0 \end{pmatrix}}{2} </math>

Ya podemos sustituir en la fórmula de la tensíon para unos coeficientes de lamé ẞ=ẟ=1 :

<math> \sigma= \begin{pmatrix} cosθ b & 0 & 0 \\ 0 & cosθ b & 0 \\ 0 & 0 & cosθ b \end{pmatrix} + \begin{pmatrix} 0 & a-\rho senθ b & 0 \\ a-\rho senθ & 2cosθ b & 0 \\ 0 & 0 & 0 \end{pmatrix}= \begin{pmatrix} cosθ b & a-\rho senθ b & 0 \\ a-\rho senθ b & 3cosθ b & 0 \\ 0 & 0 & cosθ b \end{pmatrix} </math>
Operando obtenemos que:

<math> \sigma= \begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix} </math>

[[Archivo:demasiadastensiones.jpeg|centre|1200px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);

%Trazado de la tensión normal grho*tensión*grho
t1=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta);
subplot(1,3,1)
surf(MX,MY,t1)
title('grho*tensión*grho')

%Trazado de la tensiónn normal (grho/rho)tensión(grho/rho)
t2=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta).*(1./mrho.^2);
subplot(1,3,2)
surf(MX,MY,t2)
axis equal
title('(grho/rho).*tensión(grho/rho)')

%Trazado de la tensiónn normal k*tensión*k
t3=((4.*mrho.^3-3.*mrho.^2-1)./60.*mrho.^2).*cos(mtheta);
subplot(1,3,3)
surf(MX,MY,t3)
axis equal
title('k*tensión*k')
}}

==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==

El círculo de Mohr es un método gráfico que nos ayuda a determinar el estado tensional en los distintos puntos de un cuerpo
Gracias a este círculo podemos saber cuáles son las tensiones principales de nuestro problema.

[[Archivo:5.png|right]]
{{matlab|codigo=
%círculo de mohr
sigx=75;
sigy=45;
sigz=50;
tauxy=50;
R=sqrt(((sigx-sigy)/2)^2+(tauxy)^2)
C=(sigx+sigy)/2;
sig1=C+R
sig2=C-R
sigma=linspace(sig2,sig1,1001);
tau1=sqrt(R^2-(sigma-C).^2);
tau2=-tau1
plot(sigma,tau1,'k',sigma,tau2,'k')
xlabel('Sigma')
ylabel('tau')
axis('square')
hold on
plot([sig1,sig2],[0,0],'k')
plot([sigx,sigy],[-tauxy,tauxy],'r')
grid on
x1= sig2:0.15:sig1+0.25
y1= (sigz/(R-C))*x1-(sigz/(R-C))*C
x2=sig2:0.15:sig1
y2=0
x3=-C
y3=-R:0.15:R
plot(X,Y,x1,y1,x2,y2,x3,y3)
hold off
}}

==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T20:28:48Z

Julio Montero Domínguez: /* Rotacional del campo de desplazamientos */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}
==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
Calculamos la divergencia para ver como se produce el desplazamiento de la placa en cada uno de sus puntos. Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \rho \vec u_\rho}{\partial \rho} +\frac{1}{\rho} \frac{\partial \rho \vec u_\theta}{\partial \theta} +\frac{1}{\rho} \frac{\partial \rho \vec u_z}{\partial z}
=\frac{1}{\rho}0 + \frac{1}{\rho} \rho (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2}) cos(\theta) + 0
</math>
Y realizando estos cálculos obtenemos el resultado final para la divergencia:

[[Archivo:divergencia7.png|right|550px]]

<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==
Ahora calcularemos el módulo del rotacional en todos los puntos de la placa:

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))
}}

[[Archivo: Captura rotacional.PNG|1000px]]

Como se puede apreciar en la imagen, existe un mayor rotacional en la parte superior central. En cambio, en la parte inferior el rotacional es próximo a 0, es decir, la placa casi no rota sobre si misma.

==Tensión de deformaciones==
==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==
==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T20:23:29Z

Julio Montero Domínguez: /* Rotacional del campo de desplazamientos */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}
==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
Calculamos la divergencia para ver como se produce el desplazamiento de la placa en cada uno de sus puntos. Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \rho \vec u_\rho}{\partial \rho} +\frac{1}{\rho} \frac{\partial \rho \vec u_\theta}{\partial \theta} +\frac{1}{\rho} \frac{\partial \rho \vec u_z}{\partial z}
=\frac{1}{\rho}0 + \frac{1}{\rho} \rho (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2}) cos(\theta) + 0
</math>
Y realizando estos cálculos obtenemos el resultado final para la divergencia:
<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3-3\rho ^2-1}{60\rho ^2})</math>
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==
Ahora calcularemos el módulo del rotacional en todos los puntos de la placa:

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))
}}

[[Archivo: Captura rotacional.PNG|1000px]]

==Tensión de deformaciones==
==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==
==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T20:21:17Z

Julio Montero Domínguez: /* Rotacional del campo de desplazamientos */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}
==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
Calculamos la divergencia para ver como se produce el desplazamiento de la placa en cada uno de sus puntos. Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \rho \vec u_\rho}{\partial \rho} +\frac{1}{\rho} \frac{\partial \rho \vec u_\theta}{\partial \theta} +\frac{1}{\rho} \frac{\partial \rho \vec u_z}{\partial z}
=\frac{1}{\rho}0 + \frac{1}{\rho} \rho (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2})
</math>

<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2})</math>
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==
Ahora calcularemos el módulo del rotacional en todos los puntos de la placa:

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))
}}

[[Archivo: Captura rotacional.PNG]]

==Tensión de deformaciones==
==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==
==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Archivo:Captura rotacional.PNG

2020-12-07T20:17:00Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T20:14:01Z

Julio Montero Domínguez: /* Rotacional del campo de desplazamientos */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}
==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
Calculamos la divergencia para ver como se produce el desplazamiento de la placa en cada uno de sus puntos. Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \rho \vec u_\rho}{\partial \rho} +\frac{1}{\rho} \frac{\partial \rho \vec u_\theta}{\partial \theta} +\frac{1}{\rho} \frac{\partial \rho \vec u_z}{\partial z}
=\frac{\partial \frac{yx}{10}}{\partial x} + \frac{\partial \frac{-x^2}{10}}{\partial y} + \frac{\partial 0}{\partial z}
=\frac{y}{10} + 0 + 0
=\frac{y}{10}
</math>

<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2})</math>
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==
Ahora calcularemos el módulo del rotacional en todos los puntos de la placa:

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))
}}

==Tensión de deformaciones==
==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==
==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T20:09:49Z

Julio Montero Domínguez: /* Rotacional del campo de desplazamientos */

{{ TrabajoED | Temperatura en una placa semicircular en 2-D (Grupo 19A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC20/21|2020-21]] |
Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez }}

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen comprendido entre los radios 1 y 2, en el plano y>=0.
:*La temperatura \(T(x,y)\), dada por el campo escalar <math>T(x,y)=log(y^2+2)</math>
:*Los desplazamientos <math>\vec u(x,y)</math> producidos por la acción de una fuerza <math>\vec{F}</math> determinada. 
De esta forma, si definimos <math>\vec r_{0}(x,y)</math> el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:
:<math>\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:
:<math>\vec u(\rho,\theta)=sin(\theta)*f(\rho) \vec g_\theta</math>
donde \(f(\rho)\) es una cierta función.

==Mallado==
Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar; primero, realizamos una parametrización de las coordenadas para trabajar en coordenadas cilíndricas, utilizando como variable <math>ρ</math> y <math>θ</math>.

Tomamos como ejes \((x,y) ∈ [−3,3] × [-1,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, <math>h=\frac{1}{10}</math> para las variables <math>ρ</math> y <math>θ</math>.
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Campo escalar: Temperatura==
Una vez tenemos la base, el mallado, podemos pasar a representar los campos. Empezamos con el campo escalar de la temperatura sobre la placa, con función <math>T(x,y)=log(y^2+2)</math>, representaremos tanto sus curvas de nivel así como la superficie que genera (considerando el eje z como eje de las temperaturas).
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
}}

==Gradiente de la Temperatura, campo escalar==
Utilizamos el gradiente para observar como un campo escalar varía a lo largo de la superficie donde tiene efecto. Es un vector que desarrollamos de las derivadas direccionales del campo escalar cuya dirección nos indica la mayor variación posible desde un determinado punto, que coincide con la ortogonal a la línea de nivel de dicho punto. 
En primer lugar, se procede al cálculo del gradiente:

:<big><math>\nabla{T(x,y)}=\frac{\partial{T}}{\partial{x}}{\vec{i}}+\frac{\partial{T}}{\partial{y}}{\vec{j}}={0}{\vec{i}}+{\frac{2y}{y^2}}{\vec{j}}</math></big>

Una vez calculado, podemos pasar a dibujarlo con Matlab, donde se puede ver fácilmente que los vectores del gradiente son perpendiculares a las líneas de nivel:
{{matlab|codigo=
h=0.1 %paso de muestreo
rho=1:h:2 %variable rho
theta=0:h:pi %variable theta
[mrho,mtheta]=meshgrid(rho,theta) %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
mx=mrho.*cos(mtheta)
my=mrho.*sin(mtheta)
z=zeros(size(mx)) %la z es cero por ser una placa sin espesor
mesh(mx,my,z) %representamos de tal forma que solo se vean las líneas
%Pintamos las curvas de nivel correspondientes a la temperatura
t=log(my.^2)+2
pcolor(mx,my,t)
shading flat %para quitar las líneas de cuadrícula y que no se superpongan con las de nivel
hold on
contour(mx,my,t,50,'k') %el color de las líneas de nivel sale en negro
axis([-3,3,-1,3]) %modificamos los ejes
%representamos el campo vectorial del gradiente de la temperatura
u=(0.*mx)
v=((2.*my)./((my.*my)+2))
w=(0.*mx)
%dibujamos
quiver3(mx,my,z,u,v,w)
axis([-3,3,-1,3]) %volvemos a asignar los ejes para que tanto vectores como la superficie mantengan la misma escala
hold off
}}
==Campo de desplazamientos==
Ahora calcularemos el campo de desplazamientos <math>\vec u(\rho,\theta)</math> ocasionado por la acción de una fuerza.

Para ello sabemos que:

|<big><math>\nabla </math> x <math>\vec u(\rho,\theta)</math>| = <math>\frac{2\rho -1}{10} sin(\theta)</math>

Los puntos en <math>ρ = 1 </math> no sufren desplazamiento

Primero calculamos el rotacional de <math>\vec u(\rho,\theta)</math> en función de \(f(\rho)\):

El resultado que nos sale es: <math>\frac{1}{\rho} sin(\theta) (2\rho f(\rho) +\rho^2 (f'(\rho))</math>

Ese resultado lo igualamos a <math>\frac{2\rho -1}{10} sin(\theta)</math>

El resultado es una ecuación diferencial. Resolviendola nos sale que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math>

Para saber el valor de la constante, sabemos que cuando <math>ρ</math> = 1 no hay desplazamiento, es decir, si <math>ρ</math> = 1 ;<math>\vec u(\rho,\theta)=0 </math> .

La única manera de que no haya desplazamiento es que:

\(f(\rho)\) = <math>\frac{4\rho -3}{60}</math> + <math>\frac{k}{\rho ^2}</math> = 0

Por lo tanto, el campo de desplazamiento es:

<math>\vec u(\rho,\theta) = sin(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2}) \vec g_\theta</math>

==Campo de vectores==

Utilizando el campo de desplazamiento calculado previamente, pasamos a su representación. Al ser un campo vectorial, realizaremos un programa que nos permita calcular cada vector dependiendo de las coordenadas del mallado, a la vez que representar dichos vectores sobre la placa base. Los puntos de radio mayor poseen un módulo mayor que aquellos que se encuentran en el radio interior, lo que indica que su deformación será mayor.

[[Archivo:apar5.jpg|right|500px]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
UX=inline('-x.(sin(y).^2).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
UY=inline('x.sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=UX(mrho,mtheta);
V=UY(mrho,mtheta);
hold on
mesh(MX,MY,z); %representamos el mallado
quiver(MX,MY,U,V,'r'); %representamos los vectores del campo
axis([-3,3,-1,3]) %ajustamos los ejes
hold off
}}

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

Procedemos a ver como son los puntos desplazados, realizando varios gráficos; primero uno original de la placa, luego uno de los puntos ya deformados y por último uno solapado sobre el otro. Como se puede apreciar, los puntos mas exteriores sufren más deformación que aquellos que se encuentran más adentro.
[[Archivo:21 Desplazamiento.JPG|1000px|centro]]

{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);
MY=mrho.*sin(mtheta);
z=zeros(size(MX)); %la z es cero por ser una placa sin espesor
%Describimos los componentes del vector u
RUX=inline('x.cos(y)-x.(sin(y).^2).*((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
RUY=inline('x.sin(y)+x.*sin(y).*cos(y).((4.*x.^3-3.*x.^2-1)./60.*x.^2)','x','y');
U=RUX(mrho,mtheta);
V=RUY(mrho,mtheta);
subplot(1,3,1)
%Representamos la placa antes del desplazamiento:
surf(MX,MY,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Original')
subplot(1,3,2)
%Representamos la placa después del movimiento:
surf(U,V,z);
view(2)
axis equal
axis([-3,3,-1,4]) %ajustamos los ejes
title('Desplazado')
subplot(1,3,3)
%Hacemos una comparativa
hold on
surf(MX,MY,z);
surf(U,V,z);
axis equal
view(2)
axis([-3,3,-1,4]) %ajustamos los ejes
title('Comparativa')
hold off
}}

==Divergencia del campo de desplazamientos==
Calculamos la divergencia para ver como se produce el desplazamiento de la placa en cada uno de sus puntos. Para calcular la divergencia puesto que estamos trabajando en coordenadas cilíndricas utilizaremos la siguiente fórmula:
:<math>
\nabla \cdot \vec u
=\frac{1}{\rho}\frac{\partial \vec u_\rho}{\partial \rho} + \frac{\partial \vec u_2}{\partial y} + \frac{\partial \vec u_3}{\partial z}
=\frac{\partial \frac{yx}{10}}{\partial x} + \frac{\partial \frac{-x^2}{10}}{\partial y} + \frac{\partial 0}{\partial z}
=\frac{y}{10} + 0 + 0
=\frac{y}{10}
</math>

<math>\ div\vec{u} = cos(\theta) (\frac{4\rho ^3+3\rho ^2-1}{60\rho ^2})</math>
{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%representación de la divergencia
duv=((4.*mro.^3-3.*mro.^2-1)./60.*mro.^2).*cos(mteta)
maxdiv=max(max(duv))
mindiv=min(min(duv))
surf(mx,my,duv)
axis equal
axis([-3,3,-1,3])
colorbar
}}

==Rotacional del campo de desplazamientos==

{{matlab|codigo=
h=0.1 %paso de muestreo
ro=1:h:2 %variable ro
teta=0:h:pi %variable teta
[mro,mteta]=meshgrid(ro,teta)%cruzamos
%pasamos a cartesianas para hacer la representación
%multiplicando componente por componente
mx=mro.*cos(mteta)
my=mro.*sin(mteta)
z=zeros(size(mx))%la z es cero por ser una placa
mesh(mx,my,z)%representamos de tal forma que solo se vean las líneas
axis([-3,3,-1,3])%modificamos los ejes
%dibujamos el rotacional
rot=(sin(mteta).*((1/15)-(1./(30.*mro.^2))).*mro)+2.*sin(mteta).*(((4.*mro-3)./60)-(1./(60.*mro.^2)))
surf(mx,my,rot)
colorbar;
axis([-3,3,-1,3])
maxrot=max(max(rot))

==Tensión de deformaciones==
==Tensiones tangenciales==
Realizamos los cálculos necesarios para obtener las tensiones normales al vector <math>\vec{g_p}</math>, es decir:

:<math>|\sigma\cdot\vec{g_p}-(\vec{g_p}\cdot\sigma\cdot\vec{g_p})\cdot\vec{g_p}|=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}</math>

:<math>-(\begin{pmatrix} 1, & 0, & 0 \end{pmatrix}\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} & \frac{4ρ^3-3ρ^2-1}{20ρ^2}*cos(θ) & 0 \\ 0 & 0 & \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-(\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ), & sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}, & 0 \end{pmatrix}*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})*\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=</math>

:<math>=|\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}-\begin{pmatrix} \frac{4ρ^3-3ρ^2-1}{60ρ^2}*cos(θ) \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} 0 \\ sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4} \\ 0 \end{pmatrix}|=sen(θ)*\frac{-4ρ^5+3ρ^4+4ρ^3+ρ^2-20}{60*ρ^4}</math>

==Tensión de Von Mises==
La tensión de Von Mises es un campo escalar que utilizamos para comprobar el comportamiento de un determinado material ante un esfuerzo y distinguir un comportamiento plástico de uno elástico además de comprobar de dónde proviene el fallo. Viene dado por los autovalores de la matriz de las tensiones:
:<math>\sigma_{VM}=\sqrt{\frac{(\sigma_{1}-\sigma_{2})^{2}+(\sigma_{2}-\sigma_{3})^{2}+(\sigma_{3}-\sigma_{1})^{2}}{2}}</math> 
Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:
[[Archivo:21 Tension de Von Mises.JPG|1000px|centro]]
{{matlab|codigo=
h=0.1; %paso de muestreo
rho=1:h:2; %variable rho
theta=0:h:pi; %variable theta
[mrho,mtheta]=meshgrid(rho,theta); %cruzamos
%Pasamos a cartesianas para hacer la representación multiplicando componente por componente
MX=mrho.*cos(mtheta);lx=length(rho);
MY=mrho.*sin(mtheta);ly=length(theta);

a=inline('(1./(sqrt(1+(y./x).^2))).*((4.*(x.^2+y.^2).^(3/2)-3.*(x.^2+y.^2)-1)./(60.*(x.^2+y.^2)))','x','y');
b=inline('((y./x)./(sqrt(1+(y./x).^2))).*((-4.*(x.^2+y.^2).^(5/2)+3.*(x.^2+y.^2).^2+4.*(x.^2+y.^2).^(3/2)+(x.^2+y.^2)+20)./(60.*(x.^2+y.^2).^2))','x','y');
%Calculamos la tensión
for i=1:ly
for j=1:lx
v(1,1)=a(MX(i,j),MY(i,j));
v(1,2)=b(MX(i,j),MY(i,j));
v(1,3)=0;
v(2,1)=b(MX(i,j),MY(i,j));
v(2,2)=3.*a(MX(i,j),MY(i,j));
v(2,3)=0;
v(3,1)=0;
v(3,2)=0;
v(3,3)=a(MX(i,j),MY(i,j));
[W,D]=eig(v);
Valor(i,j)=sqrt(((D(1)-D(2))^2+((D(2)-D(3))^2+((D(3)-D(1))^2))*1/2));
end
end
%Dibujo de la superficie de tensiones de Von Mises
subplot(1,3,1)
surf(MX,MY,Valor)
axis equal
title ('Vista de la superficie')
%Dibujo cenital
subplot(1,3,2)
pcolor(MX,MY,Valor)
axis([-3,3,0,2.5])
axis square

title ('Vista cenital de la superficie')
%Calculo de las tensiones máximas
subplot(1,3,3)
hold on
t=Valor(1,1:lx);
%Dibujo de la tension de forma plana
plot(MX,t,'b');
maxT=max(t);
%Resaltamos los puntos de mayor tensión
for k=1:length(t);
if t(k)==maxT
plot(MX(k),maxT,'xr','markersize',10)
end
end
axis([-3 3 -0 2]);
title('Representación del máximo')
axis([-2,2,0,0.1])
axis square

hold off
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxT)
}}

==Círculo de Mohr==
==Masa de la placa==
Por último, vamos a calcular la masa de la placa sabiendo que su densidad es: <math> 𝑑(𝑥, 𝑦) = 1 + 𝑥𝑦 log(1 + 𝑥 + 𝑦^2) </math>

Para ello parametrizaremos la superficie la placa en coordenadas cartesianas de manera que tenemos:

<math>x= u*cos(v) </math>

<math>y= u*sin(v) </math>

<math>z=0 </math>

En donde:

<math>u∈[1,2] </math>

<math>v∈[0,π] </math>

Después calcularemos el vector de posición <math> \vec r(x,y) </math>:

<math>\vec r(x,y)= u*cos(v){\vec{i}}+ u*sin(v){\vec{j}} </math>

A continuación derivamos este vector respecto de u y de v:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} = cos(v){\vec{i}} + sin(v){\vec{j}}</math>

<math>\frac{\partial{\vec r(x,y) }}{\partial{v}} = -u*sin(v){\vec{i}} + u*cos(v){\vec{j}} </math>

Ahora calculamos el siguiente producto vectorial:

<math>\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}} = u {\vec{k}}</math>

El módulo de este producto vectorial es el área de la superficie:

<math>|\frac{\partial{\vec r(x,y) }}{\partial{u}} </math> x <math> \frac{\partial{\vec r(x,y) }}{\partial{v}}| = u </math>

La integral de la densidad por el área es la masa de la placa:

Para ello, antes vamos a pasar la densidad a coordenadas cilindricas:

<math> d(u,v) = 1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v)) </math>

<math>\int_{A} d(u,v) dA=\int_{0}^{π}\int_{1}^{2} d(u,v) dudv=\int_{0}^{π}\int_{1}^{2} u(1+u^2cos(v)sin(v)log(1+u*cos(v)+u^2 sin^2(v))) dudv </math>

Como es una integral bastante compleja, utilizamos matlab:

{{matlab|codigo=

f=@(u,v) u.*(1+((u.^(2)).*cos(v).*sin(v)).*(log(1+(u.*cos(v))+((u.^(2)).*(sin(v).^(2))))));
res=integral2(f,0,pi,1,2);
fprintf('El valor de la integral es: %2.3f\n',res)
}}

Con este código obtenemos que la masa es m= 8,1757u

[[Categoría:Teoría de Campos]]
[[Categoría:TC20/21]]

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T20:01:22Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:58:32Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:54:43Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:36:44Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:35:44Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:34:52Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:29:11Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:18:44Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:13:37Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:12:25Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:10:06Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:01:53Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T19:00:12Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:59:16Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:46:21Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:38:50Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:37:41Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:34:05Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:29:27Z

Julio Montero Domínguez: /* Campo de desplazamientos */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:25:05Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:23:36Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:17:43Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:15:14Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:13:51Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:13:06Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:11:53Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:11:11Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:10:00Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:08:56Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:07:49Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:04:47Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T18:03:30Z

Julio Montero Domínguez: /* Campo de vectores */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T17:20:15Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T17:19:45Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T17:19:30Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T17:05:53Z

Julio Montero Domínguez: /* Masa de la placa */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T17:02:21Z

Julio Montero Domínguez: /* Campo de vectores */

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T17:00:36Z

Julio Montero Domínguez: /* Campo de vectores */

Archivo:Capturilla.PNG

2020-12-07T16:59:32Z

Julio Montero Domínguez:

Temperatura en una placa semicircular en 2-D (Grupo 19A)

2020-12-07T16:05:39Z

Julio Montero Domínguez: /* Campo de desplazamientos */