Abstract

In this study, linear elasticity equation is solved on different geometries and under various boundary conditions.  Classical elasticity problems like plane stress, uniaxial stress problem are implemented on regarding geometry. Also, behaviors of relatively complex structures are computed as well. The key of this project is familarize the work flow and regarding command using FEniCS to solve partial differencial equation (PDE). In general, FEniCS is easy to use, and freedom of controling the computing process is allowed at certain degree through input script. In this study, some more complex structures suffered from memory limit of local machine. The ability of FEniCS will be fully expanded on high-performance clusters.

Introduction

FEniCS is an open source computing platform for solving partial differential equations (PDEs). It enables users to quickly translate scientific models into efficient finite element code \cite{project}.  Once user specify the the variational form, domain, boundary condition and the function space, FEniCS autometically generates the global stiffness matrix and force vector and solve for the unknown displacement vector. In this study, the linear elasticity equation is applied to different geometry and boundary condition. 
The main constituents of a finite element method for the solution of a boundary-value problem are \cite{BELYTSCHKO_2008} :
  1. The variational or weak statement of the problem.
  2. The approximate solution of the variational equations through the use of finite element finctions.
We first start with the methmatical statement of the problem to be solve. The governing equation for small elastic deformation of a domain \(\Omega\)  in differencial form (or strong form) can be written as \cite{solvers} :
\(-\nabla\cdot\sigma\ =\ f\ in\ \Omega\),
\(\sigma\ =\ \lambda tr\left(\varepsilon\right)I\ +2\mu\varepsilon\),
\(\varepsilon\ =\ \frac{1}{2}\left(\nabla u\ +\left(\nabla u\right)^T\right)\),
where \(\sigma\) is the stress tensor, \(f\) is the body force per unit volume, \(\lambda\) and \(\mu\) are Lamé's elasticity parameters for material in domain \(\Omega\)\(I\) is the identity temsor, \(tr\) is the trace operator on a tensor, \(\varepsilon\) is the symmetric strain-rate temsor, and \(u\) is the displacement  vector field. 
Combine the second and the third equation, we have:
\(\sigma\ =\ \lambda\left(\nabla\cdot u\right)I\ +\ \mu\left(\nabla u\ +\ \left(\nabla u\right)^T\right)\)
Next, we convert the equation into variational form (or weak form):
\(-\int_{_{\Omega}}^{ }\left(\nabla\cdot\sigma\right)\cdot vdx\ =\ \int_{_{\Omega}}^{ }f\cdot vdx\)
\(-\int_{\Omega}^{ }\left(\nabla\cdot\sigma\right)\cdot vdx\ =\ \) \(\int_{\Omega}^{ }\sigma\ :\ \nabla vdx\ \)  \(-\int_{\partial\Omega}^{ }\left(\sigma\cdot n\right)\cdot vds\)
The colon operator is the inner product between tensors, \(n\) is the outward unit normal at the boundary. \(\sigma\cdot n\) represants the value of traction vector \(T\) at the boundary and often prescribed as boundary condition. And the value of the displacement is assumed as Dirichlet boundary condition.
We have:
\(\int_{\Omega}^{ }\sigma\ :\ \nabla vdx\ =\ \int_{\Omega}^{ }f\cdot vdx\ +\ \int_{\partial\Omega_T}^{ }T\cdot vds\)
The variational formulation can be finally summarized as:
\(a\left(u,v\right)\ =\ L\left(v\right)\ \ \ \forall v\ \in\ V\),
where:
\(a\left(u,v\right)\ =\ \int_{\Omega}^{ }\sigma\left(u\right)\ :\ \nabla vdx\),
\(\sigma\left(u\right)\ =\ \lambda\left(\nabla\cdot u\right)I\ +\ \mu\left(\nabla u\ +\ \left(\nabla u\right)^T\right)\),
\(L\left(v\right)\ =\ \int_{\Omega}^{ }f\cdot vdx\ +\ \int_{\partial\Omega_T}^{ }T\cdot vds\)
Once the variational form is defined, solving for displacement \(u\) is relatively striaight forward using FEniCS. As soon as \(u\) is computed, we can compute von Mises stress defined as:
\(\sigma_M=\sqrt{\frac{3}{2}s\ :\ s}\) 
where \(s\) is the deviatoric stress tensor given as:
\(s\ =\ \sigma-\frac{1}{3}tr\left(\sigma\right)I\)
The von Mises stress is often used in determining whether an isotropic meterial will yield when subjected to a complex loading condition \cite{stress}.
In this study, the stress is applied through specify the displacement of boundaries as Dirichlet boundary conditions. And the problem is scaled and physical parameters are dimensionsless and on scale of 1.

FEniCS Implementation Steps

A typical FEniCS Python code has the following structure: 
  1. Import essencial libraries (e.g. fenics, dolfin, mshr)
  2. Define parameters in the equation
  3. Define domain, create mesh and define function space
  4. Define boundary condition
  5. Define variational problem
  6. Compute solution 
  7. Save results for visulization
The actrually Python scripts are not included in the report. Please follow this link to check the scripts.

Results