<jupyter_start><jupyter_text>Elastodinámica parte 1Nota: a este documento lo podremos ejecutar de manera dinámica si tienen instalado:- Python 3.5 o más nuevo instalado.- [Jupyter Notebook](https://jupyter.readthedocs.io/en/latest/install.html).- [FEniCS](https://fenicsproject.org/).La visualización del mismo es óptima utilizando Jupyter Notebook. Referencias- Capítulo 26 de <NAME>, Kent-<NAME>, <NAME>, [Automated Solution of Differential Equations by the Finite Element Method](https://www.springer.com/gp/book/9783642230981). The FEniCS Book (2012), Springer.- Sacamos varias cosas del sitio de <NAME>, [seguir enlace](https://comet-fenics.readthedocs.io/en/latest/).  IntroducciónEn esta sección trataremos el problema elastodinámico, es decir, a un sólido se le aplica una fuerza que puede ser variante en el tiempo. Esto hará que se pueda acelerar por la segunda ley de Newton. Entonces se debe aplicar la ecuación de balance de momento. Como ya sabemos, el sólido se acelerará pero ese movimiento deberá amortiguarse, esto se puede modelar de varias formas, lo comentaremos en dos partes:  1. Utilizando una **relación constitutiva** (relaciona la tensión $\sigma$ y la deformación $\varepsilon$) que contiene la viscosidad. Comenzaremos por este por su similitud con lo que venimos aprendiendo.2. Agregando una viscosidad a conocida como **amortiguamiento de Rayleigh**. Este método es el abordaje más comunmente utilizado en estructuras.Tal como nos pasó en la [sección anterior](https://github.com/rirastorza/Intro2FEM/blob/master/Problemas_dependientes_del_tiempo/Introduccion_problema_termico.ipynb), aparece una derivada, aunque en este caso es la derivada segunda del desplazamiento ($\ddot{u}$), que es nuestra incógnita. Para esto se puede utilizar un método implícito que es incondicionalmente estable: la familia de [**métodos Newmarks o métodos $\alpha$-generalizados**](https://link.springer.com/article/10.1007/s00466-001-0273-z). Ecuación de balance de momentoEn lo que sigue trabajaremos con el enfoque 1), la formulación dura será la siguiente:$$ \left \{ \begin{array}{l} \rho \ddot{u}-\nabla\cdot \sigma = b \ \  \text{en} \ \ \Omega\times \text{I} \\ u\left(t\right)=g\text{ en } \ \partial \Omega_{D}\times \text{I} \text{ (condición de borde Dirichlet) } \\ \sigma\left(t\right)\cdot n= h\text{ en } \ \partial \Omega_{N}\times \text{I} \text{ (condición de borde Newmann) } \\ u\left(x,0\right)=u_{0}\left(x\right) \ \  \text{en} \ \ \Omega \text{ (condición inicial de posición) } \\ \dot{u}\left(x,0\right)=v_{0}\left(x\right) \ \  \text{en} \ \ \Omega \text{ (condición inicial de velocidad) } \end{array} \tag{1}\right .$$La ecuación es similar a la que vimos en la sección de [elasticidad estática](https://github.com/rirastorza/Intro2FEM/blob/master/Elementos_finitos_en_2D/ElasticidadLineal2D.ipynb), con la diferencia que aparece el término nuevo $\rho \ddot{u}$ que justamente considera la aceleración por unidad de volumen por eso aparece multiplicada por la densidad $\rho$. A estas fuerzas que aparecen en los nodos del sólido se las denomina fuerzas inerciales (obviamente $u$ es el desplazamiento y $\sigma$ es la tensión). Note que también aparece los dominios $\Omega$ y $\partial \Omega$ multiplicados por $\text{I}$, esta notación se introduce para representar que además existe un recorrido en el tiempo $\text{I}=\left[0,T\right]$, donde $T$ es el tiempo final. Modelado modificando la relación constitutiva entre $\sigma$ y $\varepsilon$La Ec. (1) relaciona la tensión con el desplazamiento, ahora debemos también relacionarlos con la deformación. Si, según el punto 1) arriba mencionado, consideramos la siguiente relación constitutiva:$$\sigma = 2\mu\varepsilon+\left(\lambda tr\left(\varepsilon\right)+\eta tr\left(\dot{\varepsilon}\right)\right)I\tag{2}.$$y a su vez sabemos que $\varepsilon = \frac{1}{2}\left(\nabla u + \nabla u^{T}\right)$ y $\dot{\varepsilon}= \frac{1}{2}\left(\nabla \dot{u} + \nabla \dot{u}^{T}\right)=\frac{1}{2}\left(\nabla v + \nabla v^{T}\right)$, tenemos la relación entre los tres: $\sigma$, $\varepsilon$, y $u$ (también con la velocidad $v$). Además, aquí aparece un **coeficiente de viscosidad ($\eta$)** que tiene en cuenta la fricciones microscópicas del material, también la fricción con el aire de la pieza, entre otras cosas. En realidad lo debemos pensar como aquello que consuma energía para que la pieza se detenga, se amortigüe. Los parámetros $\mu$ y $\lambda$ son los coeficientes de Lamé que ya hemos visto. Formulación variacionalHacemos lo de siempre, multiplicamos por la función de prueba e integramos en el dominio $\Omega$. Vamos a considerar una función de prueba $w$ (para no confundirnos con la velocidad $v$). Obviamente se asume nula ($w=0$) en $\partial \Omega_{D}$. $$\int_{\Omega}\rho \ddot{u} \ w \ dx - \int_{\Omega} \left(\nabla \cdot \sigma\right) \ w \ dx - \int_{\Omega} bw \ dx = 0 \tag{3}$$Nuevamente ocurre lo de siempre, debemos bajar el órden de la derivada utilizando la fórmula de Green. $$\int_{\Omega}\rho \ddot{u} \ w \ dx + \int_{\Omega} \sigma : \underbrace{\nabla w}_{\varepsilon\left(w\right)} \ dx-\int_{\partial\Omega_{N}}\underbrace{h}_{\sigma \cdot n} w \ ds - \int_{\Omega} bw \ dx = 0 \tag{4}$$ Discretización del tiempoHasta ahora no hemos hablado nada sobre la discretización en el tiempo. En este caso aplicaremos la familia de métodos denominados Newmark, que son muy utilizados en dinámica de estructuras. Para esto, tomamos la aproximación de segundo órden del desplazamiento en $n+1$ considerando un paso de tiempo $\Delta t$:$$u_{n+1}=u_{n}+\Delta t \dot{u}_{n}+\frac{1}{2}\Delta t^{2}\left[2\beta \ddot{u}_{n+1}+\left(1-2\beta\right)\ddot{u}_{n}\right]\tag{4}$$y de la velocidad$$\dot{u}_{n+1}= \dot{u}_{n}+\Delta t \left[\gamma \ddot{u}_{n+1}+\left(1-\gamma\right)\ddot{u}_{n}\right]\tag{5}.$$Aparecen dos parámetros ($\beta$ y $\gamma$) que, dependiendo de cómo sea su elección, será un método **incondicionalmente estable**. Despejamos $\ddot{u}_{n+1}$ de la Ec. (4) y se obtiene:$$\ddot{u}_{n+1}= \frac{1}{\beta\Delta t^{2}}\left(u_{n+1}-u_{n}-\Delta t \dot{u}_{n}\right)-\left(\frac{1}{2\beta}-1\right)\ddot{u}_{n}\tag{6}$$y reemplazándola en Ec. (5) se tiene:$$\dot{u}_{n+1}=\frac{\gamma}{\beta\Delta t}\left(u_{n+1}-u_{n}\right)-\left(\frac{\gamma}{\beta}-1\right)\dot{u}_{n}-\Delta t \left(\frac{\gamma}{2\beta}-1\right)\ddot{u}_{n}$$si observamos esta última ecuación, nos damos cuenta que lo único que no conocemos es $u_{n+1}$, el resto de las variables corresponden al paso anterior $n$ de la iteración. Nota: existe también una generalización de esta familia de métodos, que se denomina métodos $\alpha$-generalizados, de los cuales hablaremos más adelante.  Barra empotrada en 3DComo ejemplo de aplicación usaremos la barra empotrada, pero ahora en 3D. [](barra3D.png)Las diferencias con el ejemplo anterior de elasticidad estática es que ahora, además de la fuerza peso de la barra, aparece la tensión en la dirección negativa de $z$ en la punta de la barra. CódigoTodo el código siguiente está implementado en el *ejemplo16.py* subido a la carpeta de ejemplos. Comenzamos como siempre, importando las librerías y definiendo los parámetros, en este caso, mecánicos ($\rho$, $E$, $\nu$ y lo nuevo, la viscosidad $\eta$) y de tiempo ($\Delta t$, $\beta$ y $\gamma$).<jupyter_code>from dolfin import *

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

E = 1000.0
nu = 0.3
mu = E/(2.0*(1.0+nu))
lmbda = E*nu/((1.0+nu)*(1.0-2.0*nu))
#Densidad de masa y coeficiente de viscosidad
rho = 1.0
eta = 500.0

#Parametros de tiempo
beta = 0.25
gamma = 0.5
dt = 0.08
t = 0.0
T = 4.0<jupyter_output><empty_output><jupyter_text>Luego definimos las funciones que nos definen los subdominios para las condiciones de borde.<jupyter_code>def izquierda(x, on_boundary):
    return near(x[0], 0.) and on_boundary

def derecha(x, on_boundary):
    return near(x[0], 1.) and on_boundary<jupyter_output><empty_output><jupyter_text>Ahora creamos una función para la discretización del tiempo que actualiza el desplazamiento, la velocidad y la aceleración según las Ecs. (5) y (6).<jupyter_code>def actualiza(u,u0,v0,a0,beta,gamma,dt):
    u_vec,u0_vec = u.vector(),u0.vector()
    v0_vec,a0_vec = v0.vector(),a0.vector()
    
    #Actualizar la aceleración y la velocidad
    a_vec = (1.0/(beta*dt*dt))*(u_vec-u0_vec-v0_vec*dt)-(1.0/(2*beta)-1.0)*a0_vec
    
    v_vec = v0_vec+dt*(gamma*a_vec+(1.0-gamma)*a0_vec)
    
    #Acualiza tn <-- tn+1
    v0.vector()[:],a0.vector()[:] = v_vec,a_vec
    u0.vector()[:] = u.vector()<jupyter_output><empty_output><jupyter_text>Luego creamos la malla, y los campos del paso actual y el siguiente.<jupyter_code>mesh = BoxMesh(Point(0., 0., 0.), Point(1., 0.1, 0.04), 60, 10, 5)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)

u1 = TrialFunction(V)
w = TestFunction(V)

#Campos del el paso previo (desplazamiento, velocidad, aceleración)
u = Function(V)
u0 = Function(V)
v0 = Function(V)
a0 = Function(V)<jupyter_output><empty_output><jupyter_text>Noten que aparece tres variables una que es una trial u1, y dos funciones u y u0.Ahora noten cómo se definen las condiciones de contorno, y por último la tensión aplicada en la parte derecha de la barra como una *Expression*.<jupyter_code>boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
force_boundary = AutoSubDomain(derecha)
force_boundary.mark(boundary_subdomains, 3)

dss = ds(subdomain_data=boundary_subdomains)

#Condición de borde del costado izquierdo
zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, izquierda)

#Carga aplicada
p0 = 1.
cutoff_Tc = T/5
#NOTAR COMO SE DEFINE LA FUNCIÓN!
h = Expression(("0", "t <= tc ? p0*t/tc : 0","0"), t=0, tc=cutoff_Tc, p0=p0, degree=0)<jupyter_output><empty_output><jupyter_text>Ahora comenzaremos a escribir la formulación variacional, aunque debemos definir: la velocidad (v1) y aceleración (a1) según la discretización en tiempo, $h$ aplicado en el borde derecho (marcado con la etiqueta 3 en el diferencial dss) y la relación constitutiva definida por la Ec. (2). Noten que como el medio es isotrópico y homogéneo se tomó: $\varepsilon = \nabla u$ y $\dot{\varepsilon} = \nabla \dot{u} = \nabla v$.<jupyter_code>#Velocidad y aceleración en tn+1
v1=(gamma/(beta*dt))*(u1-u0)-(gamma/beta-1.0)*v0-dt*(gamma/(2.0*beta)-1.0)*a0
a1= (1.0/(beta*dt**2.0))*(u1-u0-dt*v0)-(1.0/(2.0*beta)-1.0)*a0

#Trabajo de fuerzas externas
def Wext(w):
    return dot(w, h)*dss(3)

#Tensor de tensión
def sigma(u,v):
    return 2.0*mu*sym(grad(u))+(lmbda*tr(grad(u))+eta*tr(grad(v)))*Identity(u.geometric_dimension())<jupyter_output><empty_output><jupyter_text>Finalmente escribimos la formulación variacional tal como la Ec. (4).<jupyter_code>F = (rho*dot(a1,w)+inner(sigma(u1,v1),sym(grad(w))))*dx-Wext(w)
a = lhs(F)
L = rhs(F)<jupyter_output><empty_output><jupyter_text>Note que usamos la función dot, que indica el producto escalar. Lo que sigue ahora es iterar en todos los pasos de tiempo, guardando en cada paso las variable desplazamiento en un archivo pdv, por ejemplo. También debemos actualizar en cada paso los valores de desplazamiento, velocidad y aceleración con la función que creamos (actualiza).<jupyter_code>file = File("data/desplazamiento.pvd")
u_tip = []
while t <=T:
    t += dt
    h.t = t
    print("Time: ", t)
    #problem.solve()
    solve(a == L, u, bc)
    #u = problem.solve()
    u_tip.append(u(1., 0.05, 0.0)[1])
    actualiza(u,u0,v0,a0,beta,gamma,dt)
    file<<u<jupyter_output><empty_output><jupyter_text>Luego podemos también dibujar, por ejemplo, el desplazamiento de la punta de la barra.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
u_final = np.asarray(u_tip)
tiempo = np.arange(0,dt*(len(u_final)),dt)
plt.figure()
plt.plot(tiempo,u_tip)
plt.xlabel("Time")
plt.ylabel("Tip displacement")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13838054650288292, 'fraction_numerical': 0.02189354056990056, 'mean_word_length': 4.895566502463054, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 41, 'https://': 7, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '526404', 'license': 'CC0-1.0', 'n_tokens_mistral': 4486, 'n_tokens_neox': 4338, 'n_words': 1313, 'path': 'Problemas_dependientes_del_tiempo/Problema_elastodinamico.ipynb', 'repo_name': 'gmg-utn/elementos_finitos', 'size': '11989'}