## 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 \]

## Principle of Virtual Work

## 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.