# Software profesional en Acústica 2020-21 (M2i)

*This notebook was adapted from Chapter 1 of [The FEniCS Tutorial Volume I](https://fenicsproject.org/pub/tutorial/sphinx1/) by Hans Petter Langtangen and Anders Logg, released under CC Attribution 4.0 license. It has been created by Xiangmin Jiao (University of Stony Brook University) and it is available in the repository [Unifem/FEniCS-note](https://github.com/unifem/fenics-notes).*

# The equations of vibrations in linear elasticity

Vibrations in linear elasticity is the study of how solid objects are responding to a mechanical vibration and become 
internally stressed due to prescribed time-harmonic loading conditions. It is an important problem
in modern engineering. Its corresponding PDE is a generalization of the
Helmholtz equation, and it is among one of the most popular PDEs in 
engineering. We now study its variational formulation and how to solve
this problem using FEniCS in 2D.

## PDE problem

The time-harmonic equation governing elastic vibrations of a body $\Omega$ can be written as

\begin{equation}
\label{ftut-elast-varform-equilibrium}\tag{1}
-\omega^2\rho\boldsymbol{u} -\boldsymbol{\nabla}\cdot\boldsymbol{\sigma} = \boldsymbol{f}\hbox{ in }\Omega,
\end{equation}

where $\boldsymbol{\sigma}$ is the *stress tensor*, and $\boldsymbol{f}$ is the *body force per unit
volume*, $\rho$ is the *mass density*, and $\omega$ the angular frequency. For isotropic materials, the stress tensor is further related to the deformation by 
the following two equations:
\begin{align}
\boldsymbol{\sigma} &= \lambda\,\hbox{tr}\,(\boldsymbol{\varepsilon}) \boldsymbol{I} + 2\mu\boldsymbol{\varepsilon},
\label{ftut-elast-varform-stresstrain}\tag{2}\\
\boldsymbol{\varepsilon} &= \frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^{\top}\right),
\label{ftut-elast-varform-strainu}\tag{3}
\end{align}
where $\boldsymbol{\varepsilon}$ is the *symmetric strain-rate tensor* (symmetric gradient), 
and $\boldsymbol{u}$ is the *displacement vector field*, $\boldsymbol{I}$ denotes the *identity tensor*, 
$\mathrm{tr}$ denotes the *trace operator* on a tensor, and $\lambda$ and $\mu$ 
are material properties known as *Lamé's elasticity parameters*.

We can combine (\ref{ftut-elast-varform-stresstrain}) and
(\ref{ftut-elast-varform-strainu}) to obtain
\begin{equation}
\label{ftut-elast-varform-stressu}\tag{4}
\boldsymbol{\boldsymbol{\sigma}} = \lambda(\boldsymbol{\nabla}\cdot \boldsymbol{u})\boldsymbol{I} + \mu(\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^{\top})
\end{equation}

Note that
(\ref{ftut-elast-varform-equilibrium})--(\ref{ftut-elast-varform-strainu})
can easily be transformed to a single vector PDE for $\boldsymbol{u}$, which is the
governing PDE for the unknown $\boldsymbol{u}$ (Navier's equation).  In the
derivation of the variational formulation, however, it is convenient
to keep the equations split as above.

## Variational formulation

The variational formulation of
([1](#mjx-eqn-ftut-elast-varform-equilibrium)--[3](#mjx-eqn-ftut-elast-varform-strainu))
consists of forming the inner product of
[1](#mjx-eqn-ftut-elast-varform-equilibrium) and a *vector* test function
$\boldsymbol{v}\in \hat{V}$, where $\hat{V}$ is a vector-valued test function space, and
integrating over the domain $\Omega$:

\begin{equation} 
-\omega^2\int_\Omega \rho\boldsymbol{u}\cdot \boldsymbol{v}\ \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},\label{eq:elast:varform:variational}\tag{5}
\end{equation}
where $\mathrm{d}\boldsymbol{x}$ denotes the volume measure.

Since $\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}$ contains second-order derivatives of the primary
unknown $\boldsymbol{u}$, we integrate this term by parts:

\begin{equation}
-\omega^2\int_\Omega \rho\boldsymbol{u}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{x}-\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},
\label{eq:elast:varform:integbypart}\tag{6}
\end{equation}

where the colon operator is the inner product between tensors (summed
pairwise product of all elements), $\boldsymbol{n}$ is the outward unit normal
at the boundary, and $\mathrm{d}\boldsymbol{s}$ is a measure in surface area.

The quantity $\boldsymbol{\sigma}\cdot \boldsymbol{n}$ is known as the
*traction* or stress vector at the boundary, and is often prescribed
as a boundary condition on a part $\partial_T\Omega$ of the boundary 
as $\boldsymbol{\sigma}\cdot \boldsymbol{n} = \boldsymbol{T}$, whereas the 
remaining part of the boundary would have the value of the displacements prescribed.

Substituting (\ref{eq:elast:varform:integbypart}) into (\ref{eq:elast:varform:variational}), 
we thus obtain
\begin{equation}
-\omega^2\int_\Omega \rho\boldsymbol{u}\cdot \boldsymbol{v}\ \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}.
\label{eq:elast:varform:varform_original}\tag{7}
\end{equation}

Inserting the expression ([4](#mjx-eqn-ftut-elast-varform-stressu)) for
$\boldsymbol{\sigma}$ gives the variational form with $\boldsymbol{u}$ as unknown. 

### Symmetrizing $\boldsymbol{\nabla} \boldsymbol{v}$

One can show that the inner product of a symmetric tensor $\boldsymbol{A}$ and an
anti-symmetric tensor $\boldsymbol{B}$ vanishes. Since $\boldsymbol{\sigma}$ is a
symmetric tensor, if we express $\boldsymbol{\nabla} \boldsymbol{v}$ as a sum
of its symmetric and anti-symmetric parts, only the symmetric part will
survive in the product $\boldsymbol{\sigma} :\boldsymbol{\nabla} \boldsymbol{v}$. Thus replacing $\boldsymbol{\nabla} \boldsymbol{u}$ by the symmetric gradient gives rise to the slightly different variational form

\begin{equation}
-\omega^2\int_\Omega \rho\boldsymbol{u}\cdot \boldsymbol{v}\ \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}
\label{ftut-elast-varform-sigma_inner_gradv}\tag{8}
\end{equation}

where $\boldsymbol{\epsilon}(\boldsymbol{v})$ is the symmetric part of $\boldsymbol{\nabla} \boldsymbol{v}$:
\begin{equation}
\boldsymbol{\epsilon}(\boldsymbol{v}) = \frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{v} + (\boldsymbol{\nabla} \boldsymbol{v})^{\top}\right).
\label{ftut-elast-varform-sigma_inner_eps}\tag{9}
\end{equation}

The formulation (\ref{ftut-elast-varform-sigma_inner_gradv}) is what naturally
arises from the minimization of elastic potential energy, and is a more
popular formulation than ([7](#mjx-eqn-eq:elast:varform:varform_original})).

Here, $\boldsymbol{\epsilon}$ is a useful operator. The symmetric strain-rate tensor $\boldsymbol{\varepsilon}$ in ([3](#mjx-eqn-ftut-elast-varform-strainu)) is equal to $\boldsymbol{\epsilon}(\boldsymbol{u})$. 

### Enforcing boundary conditions

Now let us consider how to enforce boundary conditions. 
For Dirichlet boundaries, we will enforce boundary-conditions strongly.
For these points, no test functions are associated with the Dirichlet nodes.

For traction boundary conditions, we will enforce the boundary condition
weakly using the variational form ([8](#mjx-eqn-ftut-elast-varform-sigma_inner_gradv)).
Similar to the Helmholtz equation, we require their corresponding test
functions $\boldsymbol{v}$ vanish along $\partial \Omega$ for interior points.
Then, the boundary integral above has no effects for points on
$\partial\Omega\setminus\partial\Omega_T$.

### Summary of variational form
In summary, the variational problem is to find $\boldsymbol{u}$ in a vector function space $\hat{V}$ such that
\begin{equation}
a(\boldsymbol{u},\boldsymbol{v}) = L(\boldsymbol{v})\quad\forall \boldsymbol{v}\in\hat{V},\tag{10}
\end{equation}
where 
\begin{align}
a(\boldsymbol{u},\boldsymbol{v}) &= -\omega^2\int_\Omega \rho\boldsymbol{u}\cdot \boldsymbol{v}\ \mathrm{d}\boldsymbol{x}
+ \int_\Omega\sigma(\boldsymbol{u}) :\varepsilon(\boldsymbol{v})\ \mathrm{d}\boldsymbol{x},\tag{11}\\
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}\tag{12},
\end{align}
and 
\begin{equation}
\boldsymbol{\sigma}(\boldsymbol{u}) = \lambda(\boldsymbol{\nabla}\cdot \boldsymbol{u})\boldsymbol{I} + \mu(\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^{\top}).\tag{13}\\
\end{equation}

## FEniCS implementation

To demonstrate the implementation, we will model a clamped beam deformed under a time-harmonic surface force on the opposite free cross-section surface. This can be modeled by setting $\boldsymbol{T}=\boldsymbol{f}=(0,0,1)$ on that boundary part $x=L$. The beam is
box-shaped with length $L$ and has a square cross section of width $W$. We
set $\boldsymbol{u}=(0,0,0)$ at the clamped end, $x=0$. The rest of the lateral boundary is
traction free; that is, we set $\boldsymbol{T} = 0$. Therefore,
$$L(\boldsymbol{v}) = \int_{\{x=L\}} \boldsymbol{f}\cdot \boldsymbol{v} \mathrm{d}\boldsymbol{s}$$
for this problem.

### Import packages

We start by importing fenics and enforcing matplotlib. In addition, we import `mshr` for mesh generation.


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

### Generate the mesh and function spaces

Our action startes by generating meshes and defining function spaces.


```python
# Create mesh and define function space
length = 1; width = 0.2

mesh = RectangleMesh(Point(0, 0), Point(length, width), 50, 10, "right/left")
V = VectorFunctionSpace(mesh, 'P', 1)
```

### Define the variational problem

The primary unknown is now a vector field $\boldsymbol{u}$ and not a scalar field,
so we need to work with a vector function space. We will use 
piecewise-linear basis functions for all the components.



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

With `u = TrialFunction(V)` we get `u` as a vector-valued finite element
function with three components for this 2D problem.

Next, we define the stress tensor and $a$. The gradient and divergence operators 
now have a prefix `nabla_`.
This is not strictly necessary in the present problem, but is
recommended in general for vector PDEs arising from continuum mechanics,
if you interpret $\boldsymbol{\nabla}$ as a vector in the PDE notation. See
the notes on `grad(u)` vs. `nabla_grad(u)` below.


```python
from ufl import nabla_div

# Define strain and stress
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

# Physical constants
rho = 1.
omega = 2*np.pi*3
beta = 1.25
lambda_ = beta
mu = 1
d = u.geometric_dimension()  # space dimension

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define a
a = -omega**2*rho*inner(u, v)*dx + inner(sigma(u), epsilon(v))*dx
```

To define the partition of the boundary in the computational domain, each face of the finite element mesh is marked with a number flag or *marker*:


```python
# Initialize subdomain and boundary markers
tol = 1e-3
clamped_boundary = CompiledSubDomain('on_boundary and near(x[0],L,tol)', L=0., tol=tol)
traction_boundary = CompiledSubDomain('on_boundary and near(x[0],L,tol)', L=length, tol=tol)

# Initialize mesh function for boundary
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)
clamped_boundary.mark(boundary_markers, 1) 
traction_boundary.mark(boundary_markers, 2)

# Define the boundary measure with the boundary markers
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
```

To define $L$, $\boldsymbol{f}=(0, 0, 1)$ is a constant vector, instead of a scalar.
Such a vector constant is specified as `Constant(0, 0, 1)` in FEniCS.


```python
# Define L
f = Constant((0, 1))
L = dot(f, v)*ds(2)
```

### Define boundary conditions

We only specify the Dirichlet boundary condition. For the boundary condition 
$u=(0, 0, 0)$, we must set a vector value to zero, not just a scalar. 
We specify the vector constant as `Constant((0, 0, 0))`.


```python
# Define boundary condition
bc = DirichletBC(V, Constant((0, 0)), boundary_markers, 1)
```

### Solve the variational problem

Finally, we can solve the problem.


```python
# Compute solution
u = Function(V)
solve(a == L, u, bc)
```

### Plot the solution

Any component of the solution can be computed or even the modulus of the vector-field solution


```python
# Plot the Finite Element approximation
def plot_solution(u):
    '''plot solution of FEM-based simulation'''
    fig = plt.figure(figsize=(10,10))
    fig = plot(u)
    plt.xlabel(r'$x$ / m')
    plt.ylabel(r'$y$ / m')
    plt.colorbar(fig, fraction=0.038, pad=0.04);
    
# Plot vector field   
plot_solution(u)

# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
Q = FunctionSpace(mesh, 'P', 1)
u_magnitude = project(u_magnitude, Q)
plot_solution(u_magnitude)
```


```python
# Plot second component of the solution
plot_solution(u[0])
plot_solution(u[1])
```

## Note on `grad(u)` vs. `nabla_grad(u)`

For scalar functions, $\boldsymbol{\nabla} u$ has a clear meaning as the vector
$$\boldsymbol{\nabla} u =\left(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y},
\frac{\partial u}{\partial z}\right).$$

However, for a vector-valued function $\boldsymbol{u}$, the meaning is less clear.
Some sources define $\boldsymbol{\nabla}\boldsymbol{u}$ as the matrix with elements
$\partial u_j / \partial x_i$, especially in continuum mechanics.
Some other sources prefer $\partial u_i / \partial x_j$, such as in the 
*Jacobian matrix*.

In our notes, $\boldsymbol{\nabla}$ always denotes the matrix 
with elements $\partial u_j / \partial x_i$, to be consistent with most PDEs.
The corresponding operator in FEniCS is `nabla_grad(u)`.
For the Jacobian matrix, its corresponding operator in FEniCS is `grad(u)`.

76b49e6dd2a20a0e1039ab3af4603d64cdc77874
ipynb
Jupyter Notebook
notebooks/FEniCS_elasticity.ipynb
maprieto/SoftwareProfesionalAcustica
d6b9bd0d0e182c0c71a878a891c7420ced1b263d
MIT
notebooks/FEniCS_elasticity.ipynb
maprieto/SoftwareProfesionalAcustica
d6b9bd0d0e182c0c71a878a891c7420ced1b263d
MIT
notebooks/FEniCS_elasticity.ipynb
maprieto/SoftwareProfesionalAcustica
d6b9bd0d0e182c0c71a878a891c7420ced1b263d
MIT
Qwen/Qwen-72B
1. YES
2. YES

