Basic FEA Theory – Part 1 – Principle of Virtual Work

Introduction

In two recent articles I presented code for short FE solvers written in Python (Full Finite Element Solver in 100 Lines of Python and Full Finite Element Solver in 200 Lines of Python). In this article I will start to explain the equations that are used in the Python code. It turns out that the foundation of the FE method is the Principle of Virtual Work, which although a bit abstract, is easy to derive. In a future article I will cover isoparametric shape functions, how to get the strains from the nodal displacements, how to get the stresses from the strains, and how to put everything together to get the global K U = R equation, which can be used to calculate the displacement field. Hopefully this will all be clear at the end of this article series 🙂.

Force Equilibrium

A stepping-stone for the derivation of the principle of virtual work is the force equilibrium equation, which is the same as Newton’s classic equation: $ F = ma$. For a solid 3D body this can be written: \[ \int_V b_i dV + \int_A t_i dA = \int_V \rho a_i dV, \] where $b_i$ is the i-th component of the body force, and $t_i$ is the i-th component of the surface traction. To simplify this equation we need to recall the definition of surface tractions: $ t_i = \sigma_{ji} n_j$, giving \[\int_V b_i dV + \int_A \sigma_{ji} n_j dA = \int_V \rho a_i dV. \] The next step is to use the divergence theorem which can be written \[ \int_V B_{i,i} dV = \int_A B_i n_i dA. \] Here, the divergence theorem can be used to convert the area integral to a volume integral, giving \[ \int_V (b_i + \sigma_{ji,j} – \rho a_i) dV = 0.\] This equation has to be true for any volume, so the integrand has to always be zero. Hence, if there is no acceleration, the equation for force equilibrium can be written: \[ \sigma_{ji,j} + b_i = 0 \]

Note: The equilibrium equation is valid for all materials!

Principle of Virtual Work

The principle of virtual work is a global (integral) form of the equilibrium equation. The derivation starts by integrating the product of the equilibrium equation with a virtual displacement field over the body: \[ -\int_V (\sigma_{ji,j} + \bar{b}_i) \delta u_i dV + \int_{A_t} (t_i – \bar{t}_i) \delta u_i dA = 0. \] In this equation both integrands are identically zero, and quantities with a bar denotes a prescribed value. It many seem odd to formulate this equation since it is clear that is will always be zero, but as I will show below, the equation can be manipulated into a different form that is quite useful. Now consider the following integral in more detail \[ \begin{align} \int_{A_t} t_i \delta u_i dA &= \int_A t_i \delta u_i dA – \int_{A_u} t_i \delta u_i dA,\\ &= \int_V (\sigma_{ji}\delta u_i)_{,j}\, dV – \int_{A_u} t_i \delta_i dA,\\ &= \int_V \sigma_{ji,j} \delta u_i dV + \int_V \sigma_{ji}\delta u_{i,j}\, dV – \int_{A_u} t_i \delta u_i dA. \end{align} \] Inserting this equation into the first equation in this section yields $$ \int_V \sigma_{ji} \delta u_{i,j} dV – \int_V \bar{b}_i \delta u_i dV – \int_{A_t} \bar{t}_i \delta u_i dA – \int_{A_u} t_i \delta u_i dA = 0.$$ But $\delta u_i=0$ on $A_u$ and $\sigma_{ij} \delta u_{i,j} = \sigma_{ij} \delta\varepsilon_{ij}$ giving \[ \int_V \sigma_{ij} \delta\varepsilon_{ij} dV – \int_V \bar{b}_i \delta u_i dV – \int_{A_t} \bar{t}_i \delta u_i dA = 0. \] Which also can be written $ \delta W_i + \delta W_e = 0 $, where $ \delta W_i $ is the internal virtual work, and $ \delta W_e $ is the external virtual work. The equation for the principle of virtual work can also be written extended to include external point forces ($f_i$): \[ \int_V \sigma_{ij} \delta \varepsilon_{ij} dV = \int_V b_i \delta u_i dV + \int_A t_i \delta u_i dA + \sum \delta u_i f_i, \] where I dropped the bar on the body force and surface tractions for simplicity.

Note: A deformable system is in equilibrium if the sum of the external virtual work and the internal virtual work is zero for virtual displacements that are kinematically admissible.

Discretize the Body Force Into Finite Elements

For finite elements calculations it is better to divide the integrals in the equation for the principle of virtual work into integrals over each finite element:

\[\displaystyle\sum_m \int_{V^{(m)}} \delta\boldsymbol{\varepsilon}^\top \boldsymbol{\sigma} dV = \sum_m\int_{V^{(m)}} \delta\mathbf{u}^\top \mathbf{b} dV + \sum_m \int_{A^{(m)}} \delta \mathbf{u}^\top \mathbf{t} dA + \sum \delta \mathbf{u}^\top \mathbf{f}\]

In this equation I also switched from index notation to vector notation. This equation looks complicated, and it may not be obvious that we have gained any particular insight here. After all, the principle of virtual work equation is typically not solved directly as it is written. The key step to solve this equation is to introduce isoparametric shape functions. Which will be covered in the next part in this series.