<jupyter_start><jupyter_text>Elastodinámica parte 2Nota: 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 20 de <NAME>, El método de los elementos finitos, Editorial Reverté, (2018). Se puede acceder por la plataforma eLibro de UTN.- Sacamos el código del sitio de <NAME>, [seguir enlace](https://comet-fenics.readthedocs.io/en/latest/).  IntroducciónEn esta sección trataremos el problema elastodinámico con el enfoque 2, lo repetimos aquí2. Agregando una viscosidad a conocida como **amortiguamiento de Rayleigh**. Este método es el abordaje más comunmente utilizado en estructuras.Las ecuaciones son similares, repetiremos aquí y mostraremos la diferencias. Ecuación de balance de momentoLa formulación dura será la misma:$$ \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 .$$ Relación constitutiva entre $\sigma$ y $\varepsilon$Consideramos la siguiente relación constitutiva:$$\sigma = 2\mu\varepsilon+\left(\lambda tr\left(\varepsilon\right)\right)I\tag{2}.$$Notemos que aquí **no** aparece un coeficiente de viscosidad que aparecía en el tutorial anterior. En el siguiente punto veremos como incorporarlo. Formulación variacionalLa formulación variacional de la Ec. (1) que utilizamos fue (ya aplicamos 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{3}.$$Esta ecuación se puede pensar como:$$m\left(\ddot{u},w\right)+k\left(u,w \right)=L\left(w\right)\tag{4}$$donde $m$ es la forma bilineal asociada a la matriz de *masa*, $k$ es la forma bilineal asociada a la matriz de *rigidéz*. Son bilineales porque dependen tanto de $u$ como de $w$. Luego podemos discretizar el espacio y esta ecuación se podrá escribir como:$$w^{T}M\ddot{u}+w^{T}K u=w^{T}F$$donde M, K, y F son matrices.Aquí se puede introducir la viscocidad en la Ec. (4), con una ecuación que involucre la velocidad $\dot{u}$.$$m\left(\ddot{u},w\right)+c\left(\dot{u},w\right)+k\left(u,w \right)=L\left(w\right)\tag{5}.$$Al discretizar en espacio aparece $C$, se conoce poco de esta matríz y porque escencialmente no se conoce la viscosidad. Una aproximación que se utiliza mucho es la que se denomina matriz de **amortiguamiento de Rayleigh**:$$C = \eta_{M} M +\eta_{K}K\tag{6}$$es una combinación de la matriz de masa y de rigidéz. $\eta_{M}$ y $\eta_{K}$ son números positivos (en la literatura también se los conocen como $\alpha$ y $\beta$). Al primero se lo llama parámetro de amortiguamiento viscoso y a $\eta_{K}$ parámetro de amortiguamiento estructural, depende de la máxima tensión de vibración y microscópicamente de la fricción interna. Discretización del tiempoPara la discretización en el tiempo comentaremos otro método interesante es el $\alpha$-generalizado. Este método consiste en resolver la evolución dinámica entre $t_{n}$ y $t_{n+1}$ como sigue:$$m\left( \ddot{u}_{n+1-\alpha_{m}},w\right)+c\left( \dot{u}_{n+1-\alpha_{f}},w\right)+k\left(u_{n+1-\alpha_{f}},w\right)=L\left(t_{n+1-\alpha_{f}},w\right)\tag{7}$$donde $x_{n+1-\alpha}=\left(1-\alpha\right)x_{n+1}+\alpha x_{n}$. Luego se utilizarán las relaciones de Newmark Ecs. (4) y (5) del [tutorial anterior](https://github.com/rirastorza/Intro2FEM/blob/master/Problemas_dependientes_del_tiempo/Problema_elastodinamico.ipynb). Aquí sólo computaremos en primera instancia la aceleración del paso $n+1$ mediante:$$\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{8}.$$sale de despejar la aceleración en $n+1$ de Ec. (4) del tutorial anterior. Luego, actualizamos la velocidad con la Ec. (5) del tutorial anterior, la repetimos:$$\dot{u}_{n+1}= \dot{u}_{n}+\Delta t \left[\gamma \ddot{u}_{n+1}+\left(1-\gamma\right)\ddot{u}_{n}\right]\tag{9}.$$**La idea general es combinar las Ecs. (7), (8), y (9) para obener una relación que permita calcular el desplazamiento en el instante $n+1$, es decir, $u_{n+1}$ en función de $u_{n}$, $\dot{u}_{n}$, y $\ddot{u}_{n}$.** Haciendo esto nos queda el siguiente sistema:$$ \bar{K}\{u_{n+1}\} = \{F(t_{n+1-\alpha_{f}})\} - \alpha_{f}K\{u_n\} - C(c_{1}\{u_n\}+c_2\{\dot{u}_n\}+c_3\{\ddot{u}_n\})-M(m_1\{u_n\}+m_2\{\dot{u}_n\}+m_3\{\ddot{u}_n\})\tag{10}$$donde:- $\bar{K} = K+c_1C+m_1M$- $c_1 = \dfrac{\gamma(1-\alpha_f)}{\beta\Delta t}$- $c_2 = 1-\gamma(1-\alpha_f)/\beta$- $c_3 = \Delta t(1-\alpha_f)(1-\dfrac{\gamma}{2\beta})$- $m_1 = \dfrac{(1-\alpha_m)}{\beta\Delta t^2}$- $m_2 = \dfrac{(1-\alpha_m)}{\beta\Delta t}$- $m_3 = 1-\dfrac{1-\alpha_m}{2\beta}.$A la matriz $\bar{K}$ se la llama matriz de **rigidez modificada**. Note además que la parte derecha de la Ec. (10) depende del desplazamiento y sus derivadas evaluadas en el instante de tiempo $t_{n}$. Estamos obteniendo el valor del desplazamiento en el tiempo $t_{n+1}$ y luego, con las Ecs. (8) y (9) actualizaremos la velocidad y aceleración, respectivamente.  CódigoTodo el código siguiente está implementado en el *ejemplo17.py* subido a la carpeta de ejemplos. El comienzo es similar al ejemplo16.py, no lo repetiremos aquí. Cambian ligeramente las funciones de actualización, que ahora podrán utilizar variables simbólicas (UFL). Además, definiremos de manera diferente las funciones $m$, $c$, $k$, y $L$, que tendrán que ver con las matrices de **masa, de amortiguamiento, de rigidéz y de carga**.Aquí mostramos las definiciones:<jupyter_code>#Tensor de tensiones
def sigma(r):
    return 2.0*mu*sym(grad(r)) + lmbda*tr(sym(grad(r)))*Identity(len(r))

#Masa
def m(u, w):
    return rho*inner(u, w)*dx

#Rigidez
def k(u, w):
    return inner(sigma(u), sym(grad(w)))*dx

#Amortiguamiento de Rayleigh
def c(u, w):
    return eta_m*m(u, w) + eta_k*k(u, w)

#Fuerzas en el borde derecho
def Wext(w):
    return dot(w, p)*dss(3)<jupyter_output><empty_output><jupyter_text>También debemos definir una función para computar los tiempos intermedios (entre muestra y muestra) necesarios para el método $\alpha$-generalizado ($x_{n+1-\alpha}$).<jupyter_code>def avg(x_old, x_new, alpha):
    return alpha*x_old + (1-alpha)*x_new<jupyter_output><empty_output><jupyter_text>También se puede formular directamente el sistema a resolver, involucrando la matriz de rigidez modificada y los coeficientes introducidos anteriormente. Dado que la matriz del sistema a resolver es la misma para cada paso de tiempo (paso de tiempo constante), no es necesario factorizar el sistema en cada paso. Se puede al principio y solo realizar el ensamblaje de la sustitución del lado derecho variable para obtener la solución de manera mucho más eficiente. Esto se hace definiendo un objeto LUSolver y pidiendo reutilizar la factorización de la matriz:<jupyter_code>res = m(avg(a0, a_new, alpha_m), w) + c(avg(v0, v_new, alpha_f), w) + k(avg(u0, du, alpha_f), w) - Wext(w)
a_form = lhs(res)
L_form = rhs(res)

K, res = assemble_system(a_form, L_form, bc)
solver = LUSolver(K, "mumps")
solver.parameters["symmetric"] = True<jupyter_output><empty_output><jupyter_text>Noten que se utiliza la función [assemble_system](https://fenicsproject.org/docs/dolfin/2017.2.0/python/programmers-reference/fem/assembling/assemble_system.html), que nos da como resultado la matriz de rigidez (la llamamos $K$ en el código, depende de $u$ y $w$) y el vector de cargas del sistema modificado (la llamamos $res$ en el código, depende solo de $w$).Ahora comenzamos el loop de los pasos de tiempo. Debe remarcarse que con el enfoque $\alpha$-generalizado cada paso de tiempo calculado no será el pasado ($n$) o el presente ($n+1$), será un intermedio $t_{n+1-\alpha_f}=t_{n+1}-\alpha_fdt$:<jupyter_code>for (i, dt) in enumerate(np.diff(time)):

    t = time[i+1]
    print("Time: ", t)
    p.t = t-float(alpha_f*dt)
    res = assemble(L_form)#Actualizo la forma lineal (es variante en el tiempo).
    bc.apply(res)#Acá lo aplico como condición de borde nueva en cada paso.
    
    solver.solve(K, u.vector(), res)#Acá resuelvo

    
    actualiza_campos(u, u0, v0, a0)# Actualizo las cantidades
    xdmf_file.write(u, t)# Las puedo guardar en un archivo si quiero.
    
    local_project(sigma(u), Vsig, sig)# Calculo la tensión
    xdmf_file.write(sig, t)

    p.t = t #Setea el tiempo en el paso siguiente.
    
    #Calculo de energías
    if MPI.comm_world.size == 1:
        u_tip[i+1] = u(1., 0.05, 0.)[1]
    E_elas = assemble(0.5*k(u0, u0))
    E_kin = assemble(0.5*m(v0, v0))
    E_damp += dt*assemble(c(v0, v0))
    E_tot = E_elas+E_kin+E_damp 
    energies[i+1, :] = np.array([E_elas, E_kin, E_damp, E_tot])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15014855687606113, 'fraction_numerical': 0.014537351443123939, 'mean_word_length': 4.667468430547204, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 19, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '454629', 'license': 'CC0-1.0', 'n_tokens_mistral': 3613, 'n_tokens_neox': 3527, 'n_words': 1036, 'path': 'Problemas_dependientes_del_tiempo/Problema_elastodinamico2.ipynb', 'repo_name': 'gmg-utn/elementos_finitos', 'size': '9438'}