Mathematical Physics - Volume II - Numerical Methods

Chapter 3. Comparison of finite element method and finite difference method

82

3.6 Program for solving of elliptical two-dimensional problems us ing finite element method 3.6.1 Uvodne napomene A program which uses finite element method will be shown in this section. It consists of complete sub programs (with insignificant changes) and solution to a number of problems given in reference [5]. Let us first define the problem that needs to be solved: A real function u ( x , y ) should be determined in the points within region Ω ⊂ R 2 , such that it satisfies the linear partial differential equation : − ∇ · ( k ∇ u )+ bu = f . (3.150) Let us mention that this differential equation is a bit more general than (3.51). The difference lies in term bu . By introducing this term, it is possible to analyse problems such as continuously supported membrane, or chemical reactions which depend on the temperature. At points within a part of boundary Ω ( ∂ Ω 1 ), function u can have known values (essential boundary conditions): u = ˆ u , whereas in the case of points on the other part of the boundary ( ∂ Ω 2 ), this function can have known values of derivatives in the direction perpendicular to the boundary (natural boundary conditions): = pu − γ . This form of natural boundary conditions is also somewhat broader than (3.54), and allows the analysis of problems where the derivative in the direction of the normal depends on the value of function u at the boundaries. Region Ω can be divided into subregions such that k , b and f can have different values in each of them. k and b are constant in these subregions. As for f , we can assume that it has the following structure f ( x , y ) = ¯ f + ˆ f δ ( x − x i , y − y i ) where ¯ f and ˆ f are constants. Even though we allow the region Ω to be divided into subregions, the jump in quantity ∇ u · n has to be equal to zero in all points at the boundaries between individual subregions. A problem defined in this way can be translated into variation form in a way shown in the previous chapters. After the finite element approximation, it is transformed into a system of linear algebraic equations in terms of unknown values of function u in nodal points. A triangular isoparametric element with six nodes and rectangular isoparametric element with nine nodes will be used here. Let us describe these elements and conventions related to them, which are necessary to follow in order to ensure that the program will function properly. The following figure shows these elements in parametric ( ξ , η ) plane, with standard node and side numbering for each of them. − k ∂ u ∂ n