In this document, we demonstrate how to solve a two-dimensional Biharmonic equation using triangular finite elements within
First, let consider the Biharmonic equation which is mathematically expressed as follow
In this study, we consider the homogeneous Biharmonic equation with the exact solution
The problem will be implemented in the rectangular domain , and we apply Dirichlet boundary conditions such that the exact solution is satisfied.
In order to formulate the finite element implementation of the Biharmonic equation, the weak formulation must be derived. The reader is referred to another document for detailed descriptions of weak formulations.
After introducing the weight or test function and changing the differential equation (1) into the integral form, we have
Integrating by parts twice and using of the divergence theorem give
where denotes the domain boundary and the outward normal derivative. Since the test functions satisfy homogeneous boundary conditions, the integrals for the last two terms of equation (3) vanish.
Therefore, an alternative version of the weak form of (1) is given by
It is noteworthy that the Galerkin method is employed in this study so that the choice of function spaces for the solution, , and the test function, , is the same. Therefore, both the solution and the test function can be approximated by the same set of shape function, as
where is the dimension of discrete approximating to function space.
Substituting the derivatives of (6) into (4), we have
We have that (7) must hold for any value of the coefficients , so the coefficients must satisfy the equations
After substituting the derivatives of (5) into (8), the weak form of (1) can be expressed as
Since we solve the Biharmonic equation with the Dirichlet boundary conditions in this study, the boundary terms in (9) vanish. This is a consequence of every test functions have to satisfy the homogeneous conditions at boundaries. Now, it can be seen that the weak formulation of the equation contains the second-order derivatives of the approximation, of the solution. Hence, in order to make the integral exists, the approximation and its first-order derivative are required to be continuous. Therefore, the -interpolations are required.
BellElement is introduced in this document to provide the continuously differentiable shape functions to interpolate the solution of a Biharmonic equation. This kind of element arises in order to satisfy continuity of the solutions of a fourth-order problem.
Specifically, we shall provide
oomph-libusing the Biharmonic equation as a case study.
BellElement provides a discretisation of the variational principle (9) with two-dimensional, three-node triangular finite elements. In this element, the subparametric idea is employed. This means that the degree of geometric shape functions (shape function approximated a geometry) is defined to be lower than the degree of polynomials for parametric shape functions (shape function approximated a variable). In this study, a geometry is approximated by the linear Lagrange polynomials while the Bell shape functions are employed to interpolate variables.
The two-dimensional linear Lagrange shape functions defined over a triangle are parametrised by the local coordinates and are defined as
Therefore, a position is approximated by
where are the nodal positions in the direction at node .
Similar to a geometric approximation, a variable can be approximated as
where are the element's 2D global coordinates. The function are the two-dimensional Bell shape functions parametrised by area coordinates which we will describe in the next two subsections.
The area coordinate system or Barycentric coordinates are defined on a triangle which is expressed in terms of its three vertices. We consider a triangle defined by three vertices as
The barycentric coordinates of point in the triangle are simply defined to be ratios of triangular areas as,
Area is the total area of the triangle and can be computed as
Also, are the sub-area of the triangle described in the figure. They can be computed as
where denotes the and coordinates of the point These ratios of areas form dimensionless coordinates in the plane defined by points
A barycentric combination of three points takes the form: where The third area coordinate can be expressed as
The three vertices have the barycentric coordinates as
The shape functions of the Bell element are defined over a triangle and derived from the quintic polynomials with 18 degrees of freedom. The element comprises three nodes with six degrees of freedom at each corner node including the value, and all first and second derivatives; The graphical description of the Bell element is illustrated below.
The shape functions of the Bell element are parametrised by area coordinates, . The subscripts of the shape functions mean that they are evaluated at node with degree of freedom where denotes degree of freedom that corresponds to the unknown value, correspond to the first derivatives with respect to the first and second global coordinates, correspond to the second derivatives with respect to the first and second global coordinates, and corresponds to the mixed derivative.
Unlike the Lagrangian interpolation functions, the Bell shape functions are derived by interpolating nodal derivatives as well as nodal displacements. Also, they have to satisfy the properties that at node and at other nodes. Furthermore, all of the first and the second derivatives of at all nodes. Also, at node and at other nodes. Similarly for the first derivative with respect to , we have that at node and at other nodes. These shape functions associated with the first derivatives are zero at all nodes for , , , . Moreover, , , and . Furthermore, these shape functions associated with the second derivatives are zero at all nodes for , , and .
The Bell shape functions are mathematically expressed as
These are shape functions associated with the first node with its 6 degrees of freedom. The shape functions at another two nodes can be obtained by performing a cyclic permutation of the Barycentric coordinates in (13). The parameters appearing in the equations can be obtained from
where are components of a vertex
Note that all the first and second derivatives can be obtained by the help of the chain rule as
and the derivatives of area coordinates with respect to global coordinates can be obtained from considering (12) with the substitutions of They are expressed as follows:
Here are graphical Bell shape functions at the first node, evaluated at . The figures show that the Bell shape function at the first node have to satisfy the properties that at node and at other nodes. Furthermore, all of the Bell shape functions corresponded to all first and second derivatives at all nodes.
Note that the Bell shape functions and also the first and second derivatives are differentiated with respect to global coordinates. Therefore, the Jacobian of mapping between global to local coordinates is no longer needed to derive the derivatives of the shape functions for this type of element.
BellElement is constructed with the subparametric idea, there are two sets of approximating functions defined in the element. The linear Lagrange polynomials defined in (10) will be inherited from
TElementShape<2,2> in order to approximate the geometry. These functions can be accessible by
In order to approximate variables in the
BellElement, the Bell shape functions defined in (13) will be overloaded from
BellElementShape::Bshape(...). This function will compute the basis functions at local coordinate and it provides -continuity.
As the global coordinates are required in the derivation of the shape functions of the
BellElementShape<>, the physical coordinates of vertices have to be passed as an argument when shape functions are overloaded. The interface of the Bell shape functions obtained from
BellElementShape<> at local coordinate is
BellElementShape<>::Bshape(s,psi,position). Similarly, the first and the second-order derivatives can be obtained at local coordinate as
Note that the vector
position is the collection of coordinates on vertices of Bell triangles. The order of nodal positions defined in the vector
position is correspond with the order of those on physical elements. Also, the vertex nodes 0,1, and 2 situate anticlockwise.
Now, to obtained the basis functions that employed to approximate variables in
BellElement, they can be accessible via
and, the first and second-order derivatives of the Bell shape functions can be accessible via
Before the numerical results obtained from solving the 2D Biharmonic equation with the Bell finite element will be illustrated, the boundary conditions have to be specified. Since the Bell finite element is employed, the boundary specifications have to correspond with degrees of freedom defined for the
BellElement in the previous section.
Since we consider the Biharmonic equation with the Dirichlet boundary conditions, we have that the physical conditions that allow to be pinned on boundaries are the value and the normal derivatives. However, there is no normal derivative defined as a degree of freedom to be specified for the Bell shape function. Therefore, the first-order derivatives have to be imposed instead and the boundary specification can be worked out from the normal derivative on that boundary.
Furthermore, the values of the unknown and the first-order derivatives that have to be specified on the boundaries of the domain for the Bell element can be obtained from
For the rest of degrees of freedom associated with the second-order derivatives, they are derived conditions that can be taken care by the natural boundary conditions. In this study, we set them to be pinned in order to reduce the number of degrees of freedom in the problem. The values of the second-order derivatives that have to be specified on the boundaries of the domain for the Bell element can be worked out from (14) and are expressed as
The following is a table illustrating the performance with the associated computational time of the Bell triangular elements. The accuracy of the solutions will base on the -norm error which is mathematically described as
where denote the exact and the finite element solutions, respectively.
Table 1: -norm error of the solution obtained from the Biharmonic implementation using the Bell element with various numbers of elements.
|Element size (h)||Number of elements||Number of dofs||Error||Time (sec)|
The figures below illustrate, respectively, the finite element solution and absolute error obtained from solving the Biharmonic equation with the
BellElement on structured meshes. This implementation employs 16 triangular elements in the mesh with the number of degrees of freedom equals to 18.
Next, plotting the logarithm of -error with that of the element size shows that the obtained rate of convergence of the Bell element is quartic (see figure 1.4).
Physical_Variables is where the source function and the exact solution are defined. The source function can be specified via
Physical_Variables::source_function() while the exact solutions are defined via
Physical_Variables::get_exact_u(). Note that the six exact solutions correspond to the six degrees of freedom defined on each node.
The driver code is very simple and short. It is where the problem is defined. In this study, the problem is constructed using the structured mesh with triangular elements in 2D. A number of nodes in an element has to be specific as a template parameter in the problem set up. This is crucial in order to take care of element nodes which we will describe in another tutorial. Following the usual self-test, we call the function
BiharmonicProblem::newton_solve() to compute the solution of the problem.
BiharmonicProblem::doc_solution() will output the solution afterwards.
The problem class has four member functions, illustrated as follows:
action_before_newton_solve(): Update the problem specifications before solve. Boundary conditions maybe set here.
action_after_newton_solve(): Update the problem specifications after solve.
doc_solution(): Pass the number of the case considered, so that output files can be distinguished.
From the above mentioned functions, only the problem constructor is non-trivial.
In the present problem, the function
Problem::actions_after_newton_solve() is not required, so it remains empty. Also, the class includes a private data member which store a pointer to a source function.
The problem constructor starts by overloading the function
Problem::mesh_pt() and set to the specific mesh used in this problem. In this tutorial, we implement the problem with 2D triangular structured mesh which is internally created by
SimpleRectangularTriMesh<>. The generated output will be used to build
Next, the boundary conditions of the problem will be taken care. We pin the nodal values on all boundaries and apply the Dirichlet boundary conditions. Note that the boundary identities are labelled anticlockwise where the first boundary is associated with the side
We then loop over the elements and set the pointer to the physical parameters, the function pointer to the source function which is the function on the right-hand side of the Biharmonic equation.
We finish the constructor by assigning the equation numbering scheme.
action_before_newton_solve(), the problem specifications will be updated before performing the newton solve. The boundary values will be (re)set from the exact solutions.
The post-processing routine writes the computed results to an output file.
A pdf version of this document is available.