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

Trabajo realizado por estudiantes
Título	Temperatura en una placa semicircular en 2-D (Grupo 19A)
Asignatura	Teoría de Campos
Curso	2020-21
Autores	Guillermo Pásaro Martínez / María del Pilar Magariño Mchaar / Julio Montero Domínguez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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.

1 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].

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

2 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.

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

3 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:

[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+2}}{\vec{j}}[/math]

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:

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

4 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:

|[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]

5 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.

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

6 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.

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

7 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.

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

8 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.

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))

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.

9 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]

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')

10 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]

11 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:

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)

12 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.

%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

13 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:

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