Boundary condition in FEM method

Categories: Math

Following the previous post, I'm going to discuss about the boundary conditions in FEM method. For FEM method, there are two types of boundary conditions: 1. Natural boundary condition 2. Essential boundary condition

The difference between these two kinds of boundary conditions are well summarized in the following two references: - https://caendkoelsch.wordpress.com/2018/06/09/what-is-the-difference-between-essential-and-natural-boundary-conditions-in-fem/ - https://www.researchgate.net/post/What_is_the_difference_between_essential_boundary_conditions_and_natural_boundary_conditions - https://math.stackexchange.com/questions/2089253/how-to-identify-natural-and-essential-boundary-conditions-of-this-differential-e - http://textofvideo.nptel.ac.in/105106051/lec3.pdf

To clearly show the difference, I combined the two references:

Essential boundary conditions Natural boundary conditions
They are imposed explicitly during solving These conditions are added during the formulation of FEM problem and are automatically, “naturally” satisfied without any external conditions
They are satisfied exactly by the construction of additional trial function or modification of the linear system to solve They are satisfied up to the order of the shape function
Some examples would be
Displacement in stress analysis Bending moment or shear forces in stress analysis
Temperature in thermal analysis Adiabatic boundary in heat conduction analysis

The part that I'm mainly interested in is the natural boundary condition. For fluid dynamic problem, we usually solve a PDE like the following:

$$\frac{\partial U}{\partial t} + \nabla \cdot \vec{F} = S$$

When we impose a flux boundary condition for $\vec{F}$, we clearly get a Neumann boundary condition, i.e. Natural boundary condition for FEM method. Taking the example of 1D steady state heat transfer equation:

$$\frac{d}{dx}\left( -k \frac{dT}{dx} \right) = f$$

If we want to impose the heat flux at the boundary, we have:

$$q_L = -\left. k \frac{dT}{dx} \right|_{x=0}$$

However, the problem comes! In FVM, we can exactly make the boundary condition satisfied. In FEM, the method doesn't achieve the goal. Then the question becomes how much error it gives, or what it applies to the boundary.

We first derive the weak form the original equation:

$$\begin{gathered} (-kT'',v) = (f,v)\\ -kT'_R\cdot v_R + kT'_L\cdot v_L - (-kT',v') = (f,v) \end{gathered}$$

With Dirac delta function, we can write the above one in the abstract variational form:

\begin{aligned} a(T,v) & = L(v) \\ \textrm{where:}\; a(T,v)& = (kT', v') \\ L(v) &= (f + kT'_R\cdot \delta(x-x_R) - kT'_L \cdot \delta (x-x_L), v) \end{aligned}

To investigate what FEM imposes for the boundary, we need to find out the equation with the term $T_L$. Let's consider the simplest P1 shape function, hat function as the following:

\begin{aligned} {\varphi}_i(x) = \left\lbrace\begin{array}{ll} 0, & x < x_{i-1},\\ (x - x_{i-1})/(x_{i} - x_{i-1}), & x_{i-1} \leq x < x_{i},\\ 1 - (x - x_{i})/(x_{i+1} - x_{i}), & x_{i} \leq x < x_{i+1},\\ 0, & x\geq x_{i+1}{\thinspace .} \end{array} \right. \end{aligned}

The figure above shows the shape functions and the function we want to approximate.

When we take the test function $v$ as the first test function, we get:

$$\int_{x=0}^{x=h} \left(-k \dfrac{dT}{dx}\right) \cdot \frac{1}{h} dx= \int_{x=0}^{x=h} f(x)\cdot \left(- \frac{x}{h} \right)dx -kT'_L$$

Physically, the above equation tells that the average heat flux equals to the summation of momentum of source term and the heat flux on the left boundary.

If we take the limit $h \rightarrow 0$, then the first term on the right diminishes and the left hand side becomes the heat flux on the left boundary.

🔔Conclusion, the natural boundary condition in FEM method is only matched in the "average" sense. Unless the boundary element has a length equal to zero, or you have a shape function exactly match the analytical solution, we can never match the flux on the boundary