Grupo B12 Trabajo 5: Estudio de tensiones sobre una placa mediante diferenciación automática

Trabajo realizado por estudiantes
Título	Estudio de tensiones sobre una placa mediante diferenciación automática. Grupo B12
Asignatura	Teoría de Campos
Curso	2021-2022
Autores	Víctor Sillero González Miguel Ureña Pliego Rodrigo Vázquez Pérez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El trabajo correspondiente a este grupo es el número 5, en particular consiste en el estudio de temperaturas, deformaciones y tensiones sobre una placa compuesta de medio anillo circular. La peculiaridad de este trabajo radica en que se ha realizado mediante diferenciación automática gracias a una librería realizada por los componentes del grupo. Dicho programa ha sido realizado en el leguaje de programación "Julia", y por ende todo el desarrollo ha sido realizado mediante dicho lenguaje.

Todos los gráficos se han creado a partir de la librería Makie de Julia. Se ha creado una función con los estilos que queríamos usar que toma como input campos y funciones así como una lista de puntos y algunos parámetros de estilo y dibuja la gráfica correspondiente. Al ser tanto el programa para las gráficas como la librería de diferenciación automática demasiado largos para ser publicados en esta página se adjunta un link donde se pueden descargar y ver. Aunque todos las gráficas mostradas se han producido por métodos computacionales, también se han realizado los cálculos a mano para poder mostrar las ecuaciones de los campos y tensores calculados.

La librería utilizada a lo largo de todo el trabajo queda adjunta en el anexo.

En este caso los datos aportados por el enunciado han sido los siguientes:

• Placa plana que ocupa la mitad de un anillo circular centrado en el origen y comprendido entre los radios 1 y 2, en el plano y ≥ 0.
• Función temperatura en la placa: [math] T(x,y) = x^2 + (y-1)^2 [/math].
• Campo de deformaciones: [math]\vec u(ρ,θ) = \frac{ρ}{10} sin(θ) + \frac{ρ}{5} cos(θ) [/math].

Portada programa Julia

1 Introducción

Este trabajo tiene como objetivo el estudio de las tensiones y fuerzas que producen desplazamientos y deformaciones en una placa, así como su campo de temperaturas a partir de un campo de desplazamientos y una función escalar de temperaturas. Además, se va a estudiar el cálculo de derivadas y operadores diferenciales mediante métodos computacionales.

1.1 Diferenciación automática

1.1.1 Métodos habituales: Aproximación numérica

Normalmente cuando se quiere obtener el valor de la derivada de una función en un punto mediante un ordenador se realiza una aproximación a partir de la definición de derivada. [math] (f(x+h) - f(x)) / h [/math] dando un valor concreto lo menor posible para h. Un ejemplo del método con la función [math] f(x) = x^2 [/math] y [math] h = 10^{-10} [/math] en [math] x = 3 [/math]

h = 10^-10
f(x) = x^2
derivada = (f(3+h) - f(3)) / h

El resultado dado no es exactamente 6 pero se consiguen 6 decimales de precisión por lo que es una buena aproximación. Este algoritmo realiza 3 operaciones. Evaluar f en x + h, evaluar f en x y dividir todo ello entre h.

Un problema del método es si se emplea una h demasiado pequeña puede dar problemas.

Si se supera el número más pequeño que se puede expresar con el tipo de dato utilizado (Float64 en nuestro caso) el resultado mostrado será 0. Por lo tanto un manejo correcto de los valores de h es esencial y complica el método bastante.

1.1.2 Números duales

El método empleado en este trabajo es exacto. No se realiza ninguna aproximación. Para ello es necesario definir un nuevo tipo de número distinto de los reales o los complejos, los números duales.

Los números duales son números del tipo a + bε siendo a y b números reales o duales.

Denominamos parte real al número a y parte derivada al número b.

Además se define la siguiente propiedad. ε es distinto de 0 y [math] ε^2 = 0 [/math]. ε por lo tanto no es un número real.

Las operaciones aritméticas básicas se definen de manera sencilla.

• Suma: [math] (a + bε) + (c + dε) = (a + c) + (b + d)ε [/math]
• Producto [math] (a + bε) * (c + dε) = ac + (bc + ad)ε + (bd)^2 * ε^2[/math] , como [math] ε^2 = 0 [/math] el resultado final es [math] (a + bε) * (c + dε) = ac + (bc + ad)ε[/math]

1.1.3 Series de Taylor y números duales

Sea f una función cualquiera y a un punto del dominio de f la serie de Taylor de f en el punto a es

[math] f(x) = f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^{2} + ... = \sum_{i=0}^\infty \frac{f^{i}(a)}{i!}(x-a)^{i} [/math] Si se evalúa la serie de Taylor de f en el valor dual x = a + ε se obtiene

[math] f(a + ε) = f(a) + f'(a)ε + \frac{f''(a)}{2!}ε^{2} + ... [/math] y como [math] ε^{2} = 0 [/math] el resultado final es

[math] f(a + ε) = f(a) + f'(a)ε [/math]

La parte real de [math] f(a + ε) [/math] es f evaluado en a, es decir f(a) y la parte derivada de f(a + ε) es la derivada de f evaluada en a, es decir f'(a).

1.1.4 Cálculo de derivadas por números duales

Gracias a las propiedades comentadas en la sección anterior podemos establecer un nuevo método esta vez exacto para calcular derivadas.

[math] f'(a) = parte Derivada[f(a + ε)] [/math]

Este método consiste en únicamente una operación siendo no solamente exacto sino además más rápido que el método aproximado.

1.1.5 Programación

En primer lugar se define un nuevo tipo de dato Dual, con las propiedades de parte real x y parte derivada ε.

mutable struct Dual{T<:Real} <: Real 
    x::T
    ϵ::T
end

Ahora se definen las operaciones básicas sobre el tipo de dato dual. Aquí se muestran algunas de ellas.

#negación 
Base.:-(a::Dual) = Dual(-1 * a.x, -1* a.ϵ)

#suma
Base.:+(a::Dual, b::Dual) = Dual(a.x + b.x, a.ϵ .+ b.ϵ)

#diferencia
Base.:-(a::Dual, b::Dual) = Dual(a.x - b.x, a.ϵ - b.ϵ)

#producto
Base.:*(a::Dual, b::Dual) = Dual(a.x * b.x, b.x * a.ϵ + a.x * b.ϵ)

#reciproco 
recip = a::Dual -> Dual(1. / a.x, (-1* a.ϵ) / (a.x ^2))

#cociente
Base.:/(a::Dual, b::Dual) = a * recip(b)

#exponencial, potencia y raíz_cuadrada
Base.:exp(a::Dual) = exp(a.x) * Dual(1.,a.ϵ)

#logaritmo 
Base.:log(a::Dual) = Dual(log(a.x), a.ϵ / a.x)

#trigonometria
Base.:sin(d::Dual) = Dual(sin(d.x), d.ϵ * cos(d.x))
Base.:cos(d::Dual) = Dual(cos(d.x), - d.ϵ * sin(d.x))

Y ya podríamos calcular derivadas. Por ejemplo creando la función f(x) = x^2 y evaluándola en el punto Dual(4, 1) se obtiene el número Dual(16, 8), 8 es la derivada de x^2 en el punto 4.

f(x) = x^2 
f(Dual(4., 1.))

A partir de este código se programa una librería capaz de calcular gradientes, divergencias y rotacionales de cualquier función. Además es capaz de trabajar en coordenadas cilíndricas esféricas y cartesianas.

Por ejemplo el tensor gradiente de una función u en un punto en la librería se calcula con el comando:

u.grad(punto)

Existen algoritmos más avanzados que hacen este proceso aún más rápido pero para cubrir nuestras necesidades la velocidad aportada por esta librería es suficiente.

Para realizar el trabajo se define la placa, la función de temperaturas y el campo de desplazamientos y el programa crea todas las gráficas necesarias, por lo que nuestro programa es capaz de resolver cualquiera de los trabajos que se han planteado de manera instantánea.

λ = 1
μ = 1

rango = 0.1

ρ_s = range(1, 2, 11) # discretización de variables
θ_s = range(0, pi, 32)
z_s = zeros(size(ρ_s))

#campo de temperaturas
t = x -> x[1]^2 +(x[2] - 1)^2 #cartesianas

T = funcion_escalar(t, tipo = 1)

u = funcion_vectorial(x -> [(x[1] / 10) * sin(x[2]), (x[1] / 5) * cos(x[2]), 0], tipo = 2) # cilíndricas

2 Dibujo del mallado

En primer lugar es necesario dibujar la placa que se va a estudiar. La placa se define en coordenadas cilíndricas (base física cilíndrica) con la siguiente parametrización: [math] ρ = u -\gt [1,2] ; θ = v -\gt [0, π]; z = 0[/math]. En el caso de que la placa fuese más compleja, nuestro programa también aceptaría definirla mediante una función que tomando un punto como entrada, devuelva 1 o 0 dependiendo de si dicho punto se encuentra en la placa o no.

Los ejes tomados para dibujar las gráficas han sido [−2, 2] × [0, 2] y se ha dibujado un mallado con el comando "wireframe". Por otra parte, el muestreo utilizado ha sido [math]h= \frac{1}{10}[/math].

Placa representada

3 Campo de las temperaturas en la placa

3.1 Dibujo de las curvas de nivel

Una vez definida la función temperatura T(x,y) en cartesianas, una parametrización para la placa y habiendo elegido el número de curvas a representar, (20), mediante la función "contour" se han dibujado las curvas de nivel.

Curvas de nivel de la Temperatura

3.2 Máximo de temperatura

Para obtener el máximo de la temperatura con precisión la función "dibujar2D" que empleamos para dibujar las gráficas devuelve el mayor valor e indica dicho punto mediante un aspa roja.

3.3 Dibujo del gradiente como campo vectorial ([math] \nabla [/math]T)

De la función de la temperatura dada, se ha obtenido el gradiente: [math] 2x \vec i + 2 (y - 1) \vec j [/math].

Dibujado junto a las curvas de nivel, se observa fácilmente la perpendicularidad entre dichas curvas y los vectores representados. Esto tiene sentido puesto que el gradiente indica la dirección de mayor crecimiento mientras que su perpendicular indica el crecimiento nulo. Esta última dirección es, por tanto, tangente a las curvas de nivel.

Gradiente

4 Tensiones y fuerzas

4.1 Posición inicial y final de la placa

Como punto de partida se dibuja la posición inicial y final de la placa. Los puntos de la posición final se han calculado sumando a cada punto de la placa inicial su vector desplazamiento dado por la función [math] \vec u [/math]. Además, restando a la placa final la media de los desplazamientos en todos los puntos se ha obtenido una gráfica que muestra únicamente las deformaciones que la placa experimenta. De esta forma, permite ver que la placa final está achatada y deja de ser circular.

Comparación de la placa inicial y final con deformación de la placa sin desplazamiento

El objetivo final es conocer el campo de fuerzas que ha producido este desplazamiento. Para ello a partir del campo de desplazamientos se calculan las tensiones que experimenta la placa y a partir de éstas, el campo de fuerzas.

4.2 Representación del campo vectorial de desplazamientos

Tras definir la función [math] \vec u [/math] en el código, y mediante el uso de la función "arrows" (gracias a que el campo no dispone de componente vertical), se obtiene la siguiente gráfica:

Campo de desplazamientos

4.3 Divergencia y rotacional

Calcular la divergencia y el rotacional del campo de desplazamientos es de gran utilidad, ya que estos aportan información muy valiosa sobre las características de la deformación que experimenta la placa.

El rotacional mide cuánto rotan los puntos de la placa. Su módulo cuantifica la rotación, mientras que su dirección y sentido dan el eje de rotación. La divergencia en cambio, mide cuánto cambia el área en un entrono de cada punto.

La deformación de un entorno de un punto de la placa, no su desplazamiento, es por lo tanto la rotación y el cambio de tamaño de ese entorno.

Para tener una mejor intuición de lo que esto significa, se ha dibujado el módulo del rotacional y la divergencia de una rotación de ángulo π/2, una homotecia de la placa del factor 2 y un desplazamiento sin deformación de dirección 2 [math] \vec i [/math]. En el desplazamiento como no hay deformación, el rotacional y la divergencia deben valer cero. El rotacional es mayor que cero e igual en todos los puntos y la divergencia es cero. Esto se debe a que todos los puntos rotan por igual (π/2) y su área no cambia. En la homotecia por el factor 2, el rotacional debe ser nulo y la divergencia 4, ya que cada lado se duplica produciendo una cuadruplicación del área.

Explicación divergencia y rotacional

4.3.1 Divergencia del campo desplazamientos

Calculando la divergencia de forma manual y con diferenciación automática, se ha concluido que esta es nula en cualquier punto: [math] ∇ ⋅ \vec u = \frac{1}{ρ} [\frac{δ}{δρ} (ρ (\frac{ρ}{10} sen θ)) + \frac{δ}{δθ} (\frac{ρ}{5} cos θ) + \frac{δ}{δz} (ρ 0)] = 0 [/math].

Lo que indica la divergencia es la variación del diferencial de área. En este caso concreto en que es nula, se puede apreciar en la gráfica del punto anterior como aunque los cuadrados del mallado de la placa puedan ser deformados mantienen el mismo tamaño, es decir, ni se dilatan ni se contraen.

Divergencia del campo

4.3.2 Visualización del módulo del rotacional del campo desplazamientos

Mediante el cálculo manual se ha obtenido el siguiente rotacional: [math] V × \vec u = \frac{1}{ρ} \left| \begin{matrix} \vec e_ρ & ρ \vec e_θ & \vec e_z \\ \frac{δ}{δρ} & \frac{δ}{δθ} & \frac{δ}{δz} \\ \frac{ρ}{10} sen θ & ρ \frac{ρ}{5} cos θ & 0 \end{matrix} \right| = \frac{3}{10} cos(θ) \vec e_θ [/math].

Por tanto su módulo es: [math] |V × \vec u| = |\frac{3}{10} cos(θ)| [/math].

El resultado de representarlo sobre el plano XY es el siguiente:

Rotacional del campo desplazamiento

El rotacional indica el giro experimentado por la placa en su deformación, lo que se puede apreciar en la representación anterior y teniendo en cuenta el punto 4.3. Así, los extremos de la placa, que experimentan un giro, tienen el mayor valor del rotacional. Cabe decir que estos valores tendrán signo contrario, pues los giros son opuestos. Esto no se tiene en consideración para la representación gráfica, pues en ella buscamos el módulo del rotacional, que siempre es positivo.

4.4 Tensor de tensiones

Sabiendo que se encontra en un medio lineal, isótropo y homogéneo (así lo indica el enunciado), la función del tensor de tensiones viene descrita de la siguiente forma: [math] σ = λ ∇ ⋅ \vec u 1 + 2 μ ε [/math].

En dicha fórmula conocemos que los Coeficientes de Lamé (λ y μ) que dependen de las propiedades elásticas de cada material. Tomando en este caso que ambos son 1.

La parte simétrica del tensor gradiente del campo de deplazamientos, queda también definida de la siguiente manera: [math] ε(\vec u) = \frac{∇ \vec u + ∇ \vec u ^t}{2} [/math].

De forma manual se ha obtenido que el gradiente del campo es [math] ∇ \vec u(ρ,θ)=\begin{pmatrix} \frac{1}{10} sen θ & - \frac{ρ}{10} cos θ & 0 \\ \frac{1}{5} cos θ & - \frac{ρ}{10} sen θ & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math].

Por tanto, la matriz traspuesta del gradiente del campo es: [math] ∇ \vec u^t(ρ,θ)=\begin{pmatrix} \frac{1}{10} sen θ & \frac{1}{5} cos θ & 0 \\ - \frac{ρ}{10} cos θ & - \frac{ρ}{10} sen θ & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math].

Resolviendo la parte simétrica del tensor gradiente de forma manual, se obtiene como resultado: [math] ε_{ij}(ρ,θ)=\begin{pmatrix} \frac{1}{10} sen θ & cos θ (\frac{1}{10} - \frac{ρ}{20} ) & 0 \\ cos θ (\frac{1}{10} - \frac{ρ}{20}) & - \frac{1}{10} sen θ & 0 \\ 0 & 0 &0 \end{pmatrix}[/math].

De esta forma, se concluye que el tensor de tensiones viene representado de la siguiente forma: [math] σ_{ij}(ρ,θ)=\begin{pmatrix} \frac{1}{5} sen θ & cos θ (\frac{1}{5} - \frac{ρ}{10} ) & 0 \\ cos θ (\frac{1}{5} - \frac{ρ}{10}) & - \frac{1}{10} sen θ & 0 \\ 0 & 0 &0 \end{pmatrix}[/math].

4.4.1 Representación de las tensiones normales

Representación de las tensiones normales

4.4.1.1 Tensión normal en el eje [math] \vec e_ρ [/math]

Conociendo que las tensiones en el eje [math] \vec e_ρ, [/math] son [math] \vec e_ρ ⋅ σ ⋅ \vec e_ρ [/math], calculándolo de forma manual se obtiene el siguiente resultado: [math] \frac {1}{5} sen θ \vec e_ρ [/math].

4.4.1.2 Tensión normal en el eje [math] \vec e_θ [/math]

Conociendo que las tensiones en el eje [math] \vec e_θ, [/math] son [math] \vec e_θ ⋅ σ ⋅ \vec e_θ [/math], calculándolo de forma manual se obtiene el siguiente resultado: [math] - \frac {1}{10} sen θ \vec e_θ [/math].

4.4.1.3 Tensión normal en el eje [math] \vec e_z [/math]

Conociendo que las tensiones en el eje [math] \vec e_z, [/math] son [math] \vec e_z ⋅ σ ⋅ \vec e_z [/math], calculándolo de forma manual se obtiene el siguiente resultado: [math] 0 \vec e_z [/math].

4.4.2 Tensiones tangenciales respecto al plano ortogonal a [math] \vec e_ρ [/math]

En este caso, simplemente se ha aplicado la fórmula correspondiente que es la siguiente: [math] |σ ⋅ \vec e_ρ - \vec e_ρ ⋅ σ ⋅ \vec e_ρ| [/math].

Desarrollando la fórmula: [math] \begin{pmatrix} \frac{1}{5} sen θ \\ cos θ (\frac{1}{5} - \frac{ρ}{10}) \\ 0 \end{pmatrix} - \begin{pmatrix} \frac {1}{5} sen θ \\ 0 \\ 0 \end{pmatrix} = cos θ (\frac {1}{5} - \frac {ρ}{10}) [/math].

Comparación Tensiones tangenciales

4.5 Tensión de Von Mises

La tensión de Von Mises viene determinada por la fórmula: [math] σ_{VM} = \sqrt{\frac{(\vec σ_1 - \vec σ_2)^2 + (\vec σ_2 - \vec σ_3)^2 + (\vec σ_3 - \vec σ_1)^2}{2}} [/math].

Puesto que cada punto de la placa lleva asociado un tensor de la tensión, siendo [math] \vec σ_1, \vec σ_2, \vec σ_3; [/math] los autovalores de dicha matriz, se debe ir punto a punto calculándolos, por tanto es necesario un bucle que rellene una matriz con los valores de la tensión Von Mises para cada punto.

La tensión de Von Mises muestra los puntos donde se produciría la rotura de la placa. En este caso rompería la parte superior de la placa.

La representación es la siguiente:

Tensión Von Mises

4.6 Fuerza causante del desplazamiento

Conociendo que la placa viene deformada por la función vectorial dada, se puede obtener la fuerza que motiva el desplazamiento. Esa viene definida por: [math] \vec F = - \nabla ⋅ \vec u [/math].

El resultado de la fuerza obtenido manualmente es el descrito a continuación: [math] - \frac {1}{10} sen θ \vec e_ρ + \frac {1}{5} (\frac {1}{ρ} - 2) cos θ \vec e_θ [/math]. La componente según [math] \vec e_z [/math] es nula.

Calculado por diferenciación automática queda la siguiente gráfica:

Campo de fuerzas

Si de dibuja el campo de fuerzas y la las posiciones inicial y final no queda del todo claro que éste sea el campo que motiva el desplazamiento y la deformación de la placa.

Campo de Fuerzas y desplazamientos

Si se dibuja la resultante de las fuerzas el desplazamiento de la placa queda explicado, mientras que dibujando únicamente la placa final deformada, sin considerar el desplazamiento y el campo de fuerzas, se pueden apreciar claramente que deformaciones como el achatamiento de la forma circular de la placa son causadas por el campo de fuerzas.

Campo de fuerzas, desplazamiento y deformación de la placa

5 Anexo

En este enlace se puede descargar el código que se ha empleado. Para que funcione es necesario tener instalado Julia y las librerías Makie y LinearAlgebra. https://drive.google.com/drive/folders/1aeOkuIzscmac7BiWYsRidg_zSQZc8gfB?usp=sharing