# Ecuación de la elasticidad lineal

## Ecuación diferencial parcial

Escribimos el problema de deformación elastica de un cuerpo con dominio geométrico $\Omega$ como
\begin{equation}
\rho \ddot{\boldsymbol{u}} + \rho \eta \dot{\boldsymbol{u}} - \boldsymbol{\nabla}\cdot\boldsymbol{\sigma} = \boldsymbol{f}\hbox{ in }\Omega,
\end{equation}

definiendo $\boldsymbol{\sigma}$ como el *tensor de esfuerzos*, $\boldsymbol{f}$ como la *fuerza de cuerpo por unidad de volumen* y $\boldsymbol{u}$ como el *vector de desplazamientos*. Además, las derivadas temporales se denotan aquí con un punto sobre la variable.

Restringiendo el problema a materiales isotrópicos, podemos definir el tensor de esfuerzos y su relación con las deformaciones como
\begin{align}
\boldsymbol{\sigma} &= \lambda\,\hbox{tr}\,(\boldsymbol{\varepsilon}) \boldsymbol{I} + 2\mu\boldsymbol{\varepsilon}, \\
\boldsymbol{\varepsilon} &= \frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^{\top}\right),
\end{align}
donde $\boldsymbol{\varepsilon}$ es el *tensor simetrico de la tasa de deformación*, $\boldsymbol{I}$ denota el *tensor identidad*, $\mathrm{tr}$ denota la *traza* en un tensor, $\lambda$ y $\mu$ 
son los *parametros de Lamé* (propiedades del material) y $\rho$ es la densidad.

Si combinamos las definiciones de esfuerzo y deformación, podemos escribir el esfuerzo como
\begin{equation}
\boldsymbol{\boldsymbol{\sigma}} = \lambda(\boldsymbol{\nabla}\cdot \boldsymbol{u})\boldsymbol{I} + \mu(\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^{\top})
\end{equation}

## Formulación variacional

La formulación variacional del problema de deformación elastica la definimos como la integral sobre el dominio $\Omega$ del producto interno entre la ecuación y una función vectorial de test $\boldsymbol{v}\in \hat{V}$, donde $\hat{V}$ es el espacio de la funcion de test.

\begin{equation} 
\int_\Omega \rho \ddot{\boldsymbol{u}} \mathrm{d}\boldsymbol{x} + \int_\Omega \rho \eta \dot{\boldsymbol{u}} \mathrm{d}\boldsymbol{x} -\int_\Omega (\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}) \cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{x} = \int_\Omega \boldsymbol{f}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{x},
\end{equation}

Como $\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}$ contiene derivadas de segundo orden en $\boldsymbol{u}$, debilitamos esta parte de la ecuación integrando por partes:

\begin{equation}
-\int_\Omega (\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}) \cdot \boldsymbol{v} \ \mathrm{d}\boldsymbol{x}
= \int_\Omega \boldsymbol{\sigma} : \boldsymbol{\nabla} \boldsymbol{v} \ \mathrm{d}\boldsymbol{x} - \int_{\partial\Omega} (\boldsymbol{\sigma}\cdot \boldsymbol{n})\cdot \boldsymbol{v} \ \mathrm{d}\boldsymbol{s},
\end{equation}
donde $\boldsymbol{n}$ es vector normal en el contorno.

$\boldsymbol{\sigma}\cdot \boldsymbol{n}$ es conocido como la *tracción* (el esfuerzo en el contorno) y puede ser usado para prescribir una condicion de contorno de tipo Neumann o Robin de forma $\boldsymbol{\sigma}\cdot \boldsymbol{n} = \boldsymbol{T}$.

Combinando las 2 ecuaciones anteriores podemos escribir el problema variacional como

\begin{equation}
\int_\Omega \rho \ddot{\boldsymbol{u}} \mathrm{d}\boldsymbol{x} + \int_\Omega \rho \eta \dot{\boldsymbol{u}} \mathrm{d}\boldsymbol{x} + \int_\Omega \boldsymbol{\sigma} : \boldsymbol{\nabla} \boldsymbol{v}\ \mathrm{d}\boldsymbol{x} =
\int_\Omega \boldsymbol{f}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{x} + \int_{\partial\Omega} \boldsymbol{T}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{s}.
\end{equation}

Donde usando la definición de $\boldsymbol{\sigma}$ podemos escribir la ecuación variacional con $\boldsymbol{u}$ como incognita.

### Parte simetrica de  $\boldsymbol{\nabla} \boldsymbol{v}$

Ya que el producto de un tensor simetrico con un anti-simetrico es cero, y sabiendo que el tensor de esfuerzos $\boldsymbol{\sigma}$ es simetrico, podemos reemplazar en la ecuacion anterior $\boldsymbol{\nabla} \boldsymbol{v}$ por su parte simétrica $\boldsymbol{\epsilon}(\boldsymbol{v})$, dando como resultado

\begin{equation}
\int_\Omega \rho \ddot{\boldsymbol{u}} \mathrm{d}\boldsymbol{x} + \int_\Omega \rho \eta \dot{\boldsymbol{u}} \mathrm{d}\boldsymbol{x} + \int_\Omega \boldsymbol{\sigma} : \boldsymbol{\epsilon}(\boldsymbol{v})\ \mathrm{d}\boldsymbol{x} = \int_\Omega \boldsymbol{f}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{x} + \int_{\partial\Omega} \boldsymbol{T}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{s}
\end{equation}

### Estableciendo las condiciones de contorno

Ahora consideremos cómo hacer cumplir las condiciones de contorno.
Para los contornos con condiciones de Dirichlet, aplicaremos las condiciones de manera fuerte. Para estos puntos, las funciones de test asociadas con los nodos de Dirichlet se desvanecen.

Para condiciones de tracción en los contornos, aplicaremos la forma variacional.
Similarmente a lo realizado con la ecuación de Poisson, requerimos que el espacio funcional de test correspondiente a 
las funciones $ \boldsymbol{v} $ desaparezcan a lo largo de $ \partial \Omega $ para los puntos interiores.
Entonces, la integral de contorno anterior no tiene efectos para los puntos en
$ \partial \Omega \setminus \partial \Omega_T $.

### Forma variacional 
En resumen, el problema variacional es el siguiente: encontrar $\boldsymbol{u}$ en un espacio funcional vectorial $\hat{V}$, tal que
\begin{equation}
m(\boldsymbol{u},\boldsymbol{v}) + c(\boldsymbol{u},\boldsymbol{v}) + a(\boldsymbol{u},\boldsymbol{v}) = L(\boldsymbol{v})\quad\forall \boldsymbol{v}\in\hat{V},
\end{equation}
donde
\begin{align}
m(\boldsymbol{u},\boldsymbol{v}) &= \int_\Omega \rho \ddot{\boldsymbol{u}} \mathrm{d}\boldsymbol{x},\\
c(\boldsymbol{u},\boldsymbol{v}) &= \int_\Omega \rho \eta \dot{\boldsymbol{u}} \mathrm{d}\boldsymbol{x},\\
a(\boldsymbol{u},\boldsymbol{v}) &= \int_\Omega\sigma(\boldsymbol{u}) :\varepsilon(\boldsymbol{v})\ \mathrm{d}\boldsymbol{x},\\
L(\boldsymbol{v}) &= \int_\Omega \boldsymbol{f}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{x} + \int_{\partial\Omega_T} \boldsymbol{T}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{s},
\end{align}

## Implementación en FEniCS

Modelamos una barra 3D empotrada en un extremo y con una carga aplicada en el otro. Adicionalmente incluimos el peso de la barra.

La carga y el peso pueden ser facilmente modeladas sumando al lado derecho de la ecuacion el termino $\boldsymbol{f}=(0,-\rho g + s,0)$, en donde $g$ es la magnitud de la aceleracion por gravedad y $s(t,x)$ es una distribución de carga cualquiera. La barra es un paralelepipedo de longitud $L$ y de sección transversal $W \times H$. Fijamos el desplazamiento en la parte empotrada ($x=0$) como $\boldsymbol{u}=(0,0,0)$, y el resto de la barra la dejamos libre de tracciones $\boldsymbol{T} = 0$. Entonces, el lado derecho de la forma variacional es de la forma
$$L(\boldsymbol{v}) = \int_\Omega \boldsymbol{f}\cdot \boldsymbol{v} \mathrm{d}\boldsymbol{x}$$



### Importar paquetes

Importamos los paquetes de solucion (**dolfin**), de mallado (**mshr**), de cálculo numérico (**numpy**) para realizar operaciones numericas de proceso y posproceso, así como **matplotlib** para ver los resultados directamente durante la ejecución del programa


```python
from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm
```

### Dominio geométrico, propiedades, malla y subdominios

Definimos el dominio computacional $\Omega$, definimos las propiedades, generamos la malla y creamos subdominios para aplicar la condición de contorno Dirichlet y la carga $s(t)$.


```python
# Dominio geometrico y malla
L = 1; W = 0.1; H =0.04

mesh = BoxMesh(Point(0, 0, 0), Point(L, W, H), 60, 10, 5)
x = SpatialCoordinate(mesh)

# Propiedades del material
# Constantes de Lame
E  = 1000.0
nu = 0.3
mu    = Constant(E / (2.0*(1.0 + nu)))
lmbda = Constant(E*nu / ((1.0 + nu)*(1.0 - 2.0*nu)))
# Densidad
rho = Constant(1.0)
# Coeficiente de amortiguacion
eta = Constant(0.1)

#Subdominios
# Izquierda
def left(x, on_boundary):
    return near(x[0], 0.) and on_boundary     
# Derecha
def right(x, on_boundary):
    return near(x[0], 1.) and on_boundary
```

### Generación de los espacios funcionales sobre la malla

Una vez creada la malla, podemos crear un espacio funcional de elementos finitos $V$ sobre dicho dominio geométrico discreto:


```python
# Espacio de funciones para desplazamiento, velocidad y aceleraciones
V = VectorFunctionSpace(mesh, "Lagrange", 2)
# Espacio de funciones para esfuerzos
Vsig = TensorFunctionSpace(mesh, "DG", 0)
```

En este caso, tenemos un espacio funcional que es vectorial. Esto se puede entender como un arreglo de espacios funcionales escalares: uno para cada componente dimensional de la deformación del sólido elástico.

Usamos un elemento de tipo Lagrangeano (Lagrange) con un orden de interpolacion de primer orden. Adicionalmente creamos un espacio funcional discontinuo (DG) con el que representaremos los esfuerzos en el postproceso.

### Definición del problema variacional

La incógnita principal de este problema es un campo vectorial $ \boldsymbol{u} $ y no un campo escalar. Por lo que necesitamos trabajar con un espacio funcional vectorial.


```python
u = TrialFunction(V)
v = TestFunction(V)
```

Con `u = TrialFunction(V)` definimos `u` como un espacio funcional de elementos finitos vectorial con tres componentes: uno para cada dimensión de este problema tridimensional.

### Formas variacionales como funciones

Definimos las formas variacionales $m(\boldsymbol{u},\boldsymbol{v})$, $c(\boldsymbol{u},\boldsymbol{v})$, $a(\boldsymbol{u},\boldsymbol{v})$ y $L(\boldsymbol{v})$, la funcion para calcular los esfuerzos $\boldsymbol{\sigma}(\boldsymbol{u})$ (`sigma`) y una función de proyecccion local para representar el tensor de esfuerzos en el postproceso.

Todas estas funciones las definimos como funciones de python, usando las funciones intrinsecas de FEniCS (`inner`, `grad`, `sym`).


```python
# Tensor de esfuerzos
def sigma(u):
    return 2.0*mu*sym(grad(u)) + lmbda*tr(sym(grad(u)))*Identity(len(u))

# Matriz de masa
def mmat(u,v):
    return rho*inner(u, v)*dx

# Matriz de rigidez elastica
def kmat(u, v):
    return inner(sigma(u), sym(grad(v)))*dx

# Amortiguacion de Rayleigh
def cmat(u, u_):
    return eta*mmat(u, v)

# Proyeccion local
def local_project(v, V, u=None):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dx
    b_proj = inner(v, v_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return
```

### Definición de condiciones de contorno

Especificamos la condición de contorno de Dirichlet $u=(0,0,0)$ en el subdominio izquierdo (usando la funcion `DirichletBC`) y definimos el area donde aplicamos la carga $s(t)$ (lado derecho).


```python
# Condicion de contorno izquierda (Dirichlet)
zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, left)

# Condicion de contorno derecha
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
force_boundary = AutoSubDomain(right)
force_boundary.mark(boundary_subdomains, 3)
dss = ds(subdomain_data=boundary_subdomains)
```

### Definicion de las fuerzas de cuerpo

Definimos la fuerza de cuerpo $\boldsymbol{f}$ compuesta por el peso de la barra $\rho g$ y la función temporal en el lado derecho $s(t)$


```python
# Definicion de las cargas
k = 0.5
s0 = 1.
s = Expression(("0", "s0*sin(k*pi*t)","0"), t=0, k=k, s0=s0, degree=0)
g = Constant((0.,-1,0.))
```

### Integrador temporal

Usamos el método de Newmark como integrador temporal para las derivadas temporales del problema. Creamos funciones para actualizar los campos de desplazamiento, velocidad y aceleracion en cada paso de tiempo. Tambien definimos campos vectoriales `u_old`, `v_old` y `a_old`, que vamos a usar en la integracion temporal.

Tambien definimos el tiempo total de simulación y el paso de tiempo.


```python
T  = 10.0 # tiempo total
dt = 0.5 # paso de tiempo
Nt = int(T/dt)
ti = np.linspace(0, T, Nt)

# Campos para el paso de tiempo anterior
u_old = Function(V, name="Desplazamiento")
v_old = Function(V, name="Velocidad")
a_old = Function(V, name="Aceleracion")

# Metodo de Newmark (punto medio)
gamma   = Constant(0.5)
beta    = Constant((gamma+0.5)**2/4.)

# Aceleracion
# a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
def update_a(u, u_old, v_old, a_old, ufl=True):
    if ufl:
        dt_ = dt
        beta_ = beta
    else:
        dt_ = float(dt)
        beta_ = float(beta)
    return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old

# Velocidad
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
def update_v(a, u_old, v_old, a_old, ufl=True):
    if ufl:
        dt_ = dt
        gamma_ = gamma
    else:
        dt_ = float(dt)
        gamma_ = float(gamma)
    return v_old + dt_*((1-gamma_)*a_old + gamma_*a)

# Actualizar campos
def update_fields(u, u_old, v_old, a_old):
    """Update fields at the end of each time step."""

    # Vectores de referencia
    u_vec, u0_vec  = u.vector(), u_old.vector()
    v0_vec, a0_vec = v_old.vector(), a_old.vector()

    # Actualizacion de velocidad y aceleracion
    a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
    v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)

    # Actualizacion de desplazamiento
    v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec
    u_old.vector()[:] = u.vector()
```

### Definicion del problema variacional

Usando todos los ingredientes anteriores, escribimos el problema variacional y exprezamos el problema en forma matricial $\boldsymbol{K} \boldsymbol{u} = \boldsymbol{b}$


```python
# Ecuacion de elasticidad en forma residual
a_new = update_a(u, u_old, v_old, a_old, ufl=True)
v_new = update_v(a_new, u_old, v_old, a_old, ufl=True)
res = mmat(a_new,v) + cmat(v_new,v) \
       + kmat(u,v) - rho*dot(s,v)*dss(3) - dot(rho*g,v)*dx
# Ensamble del sistema de ecuaciones
B_form = lhs(res)
L_form = rhs(res)
K, b = assemble_system(B_form, L_form, bc)
```

    Calling FFC just-in-time (JIT) compiler, this may take some time.
    Calling FFC just-in-time (JIT) compiler, this may take some time.


### Arrays de postproceso e impresion de datos

Creamos arrays de numpy para guardar algunos resultados y abrimos el archivo de **XDMF** donde vamos a almacenar los datos de postproceso


```python
# Inicializacion de variables de salida
u_tip = np.zeros((Nt,)) # desplazamiento en lado derecho de la barra
E_damp = 0
energia = np.zeros((Nt, 2))

# Creacion de archivos de salida (ParaView)
xdmf_file = XDMFFile("bar.xdmf")
xdmf_file.parameters["flush_output"] = True
xdmf_file.parameters["functions_share_mesh"] = True
xdmf_file.parameters["rewrite_function_mesh"] = False
```

## Solución del problema

En esta sección finalmente resolvemos el problema. Iniciamos definiendo nuestro campo vectorial para la solución `u` y uno donde calcularemos el tensor de esfuerzos `sig`.
Usamos el ciclo `for` para iterar sobre todos los pasos de tiempo, y en este bucle evaluamos la fuerza $s(t)$ en el paso de tiempo actual, ensamblamos la forma $L(v)$ y solucionamos el problema.
Finalmente, actualizamos los campos de desplazamiento, velocidad, aceleración y esfuerzo, y guardamos los datos necesarios para el postproceso.


```python
u = Function(V, name="Desplazamiento")
sig = Function(Vsig, name="sigma")

# Propiedades del solver
solver = LUSolver(K, "mumps")
solver.parameters["symmetric"] = True

# Bucle temporal
for i in tqdm(range(Nt),desc='Time loop'):

    t = ti[i]

    # Fuerza evaluada en t_{n+1}=t_{n+1}-dt
    s.t = t

    # Solucionar desplazamientos
    b = assemble(L_form)
    bc.apply(b)
    solver.solve(K, u.vector(), b)

    # Actualizar campos
    update_fields(u, u_old, v_old, a_old)

    # Calcular esfuerzos
    local_project(sigma(u), Vsig, sig)

    # Guardar solucion (vizualizacion)
    xdmf_file.write(u_old, t)
    xdmf_file.write(v_old, t)
    xdmf_file.write(a_old, t)
    xdmf_file.write(sig, t)

    # Guardar solucion (desplazamiento y energia)
    u_tip[i] = u(1., 0.05, 0.02)[1]
    E_elas = assemble(0.5*kmat(u_old, u_old))
    E_kin = assemble(0.5*mmat(v_old, v_old))
    energia[i, :] = np.array([E_elas, E_kin])
```

    Time loop: 100%|██████████| 20/20 [01:31<00:00,  4.59s/it]


## Post-proceso


### Visualización y postprocesamiento interno de la solución

Usando el array `u_tip`, podemos graficar el desplazamiento en el lado derecho de la barra


```python
# Imprimir desplazamiento en x=(1,0.05,0.02)
plt.figure()
plt.plot(ti, u_tip)
plt.xlabel("Time")
plt.ylabel("Desplazamiento")
plt.savefig('images/tip.png',format='png',dpi=400)
plt.show()
```

Usando las soluciones de desplazamiento `u_old` y velocidad `v_old`, podemos definir la energia elastica como $E_e = \int_\Omega \frac{1}{2} \boldsymbol{\sigma}(\boldsymbol{u}):\boldsymbol{\epsilon}(\boldsymbol{u}) dx$ y la energia cinetica como $E_k = \int_\Omega \frac{1}{2} \rho \dot{\boldsymbol{u}} \cdot\dot{ \boldsymbol{u}} dx$


```python
# Imprimir energia
plt.figure()
plt.plot(ti, energia)
plt.legend(("elastica", "cinetica", "amortiguacion"))
plt.xlabel("Time")
plt.ylabel("Energia")
plt.savefig('images/energy.png',format='png',dpi=400)
plt.show()
```

**Reconocimiento**: Este cuaderno fue adaptado de [The FEniCS Tutorial Volume I] (https://fenicsproject.org/pub/tutorial/sphinx1/) elaborado por Hans Petter Langtangen y Anders Logg, publicado bajo licencia CC Attribution 4.0.

7a12e06209db6ee7af99740862ee095da0d03df0
ipynb
Jupyter Notebook
elasticidad_lineal.ipynb
ciaid-colombia/Taller-FEniCS
629b26c4a756d1017e6b409e5ed700498a9edacd
MIT
elasticidad_lineal.ipynb
ciaid-colombia/Taller-FEniCS
629b26c4a756d1017e6b409e5ed700498a9edacd
MIT
elasticidad_lineal.ipynb
ciaid-colombia/Taller-FEniCS
629b26c4a756d1017e6b409e5ed700498a9edacd
MIT
2021-04-24T14:08:38.000Z
2021-04-24T14:08:38.000Z
Qwen/Qwen-72B
1. YES
2. YES