The Finite Element Method Applied to a Problem of Blood Flow in Vessels

We use the finite element method to solve a convection-diffusion equation when convection is dominating, a problem which describes the behavior of the concentration of a solute in a blood vessel. A new technique for computing the discrete problem is used.


Introduction
Due to the fact that various chemicals, such as oxygen, carbon dioxide, or lipids, are transported by blood to and from the other tissues, including the skin, the study of the blood flow in the human vascular system is of great interest in medicine and medical engineering.
The skin is the largest organ of the human body, having very important psychosocial implications. It represents the only system completely displayed at the body surface, offering essential information regarding the homeostasis of the internal organs and thus the general senescence process. Senescence involves a complex of factors, of which vascularisation plays an important role. The perfusion degree of the tissues depends on the microscopical structure of the vascular walls, as well as the metabolical and biochemical changes associated with age. It is well known that the collagen glycation processes that occur in the vascular walls, and are involved in the senescence process, are associated with vascularisation deficits, that are specific to the ageing process generally, and to the associated pathology specifically (atherosclerosis, Alzheimer, metabolical diseases, rheumatoid arthritis). These modifications can be observed at the skin level by means of images techniques (see [1,2]).
These are only few reasons for which the study of the possible solutes in blood vessels is so important. In this sense, many studies were made and different mathematical models were given, depending on various factors, such as the health state of the patient (the existence of pathologies like atherosclerosis, etc.). For instance, in [3], an implementation of adaptive anisotropic meshes for this class of problems by developing an a posteriori error analysis for a simpler situation, namely, a single steady advection-diffusion-reaction equation with a given convective field, is given.

The Model Problem
In order to study the transport problem of a solute in a vessel, we consider the following partial derivative equation: x, y ∈ Γ out , c = ξg 3 , x, y ∈ Γ.
Here, c = c(x, y) represents the concentration of the solute in the blood. The first term in (1) describes the diffusion, the later contains two more terms: one for the diffusion along the x-direction: n 1 (∂c/∂x), and another for the diffusion along the y-direction: n 2 (∂c/∂y). m represents the diffusivity of the solute, n a given velocity field, α is a reaction coefficient, fa possible forcing term for the solution concentration due, for example, to chemical reactions. The ratio between coefficients m and n determines the predominance of the convection or the diffusion in the physical process. If: then the convection is dominant. computing the numerical solution of a convection-diffusion problem of the form (1) becomes increasingly difficult (converge slow or not at all) as the ratio (2) increases (i.e., convection is dominant in the process), and this is the case in the model problem (1) (the ratio here is 10 4 , see Section 5). The domain Ω chosen for this problem represents a part of an artery (considered to be rigid) and is shown in Figure 1.
The borders of Ω are as follows.
(i) Γ in , which corresponds to the point in where the solution is injected in the blood, thus the concentration of the solution is known here: u = g 1 , (x, y) ∈ Γ in .
(ii) Γ out , which corresponds to the point where the concentration of the solution can be measured: u = g 2 , (x, y) ∈ Γ out .
(iii) Γ: in these points the solution is in contact with the vessel walls, ξ represents the wall permeability, and g 3 is a function that decreases from g 1 to g 2 .

The Discretization of the Problem
In this paper, we use the finite element method (see [4,5]). In order to do this, as in [6][7][8], the domain [a, b] × [c, d] is divided in rectangular subdomains. For the discretizaton of the problem we use the finite element method. The domain is divided in rectangular subdomains, having the step h x on the Ox-direction and h y on Oy. The solution of the systems generated through discretization is obtained by Gauss full elimination method. The first level on which the solution is computed is l 0 , then this particular solution is used for obtaining the solution on higher order levels. The grid on the l level is divided by the one from the l 0 level in subdomains. On each of these, the corresponding system of linear equations will be solved, and the solutions from the subdomains are used in order to generate the solution on the l level.
The partial differential equations will be replaced by a liniar system of equations through the discretization method.
In order to achieve this, and keeping the notations used in [7], we choose the grid steps h lx = (b − a)/2 l+1 and h ly = (d − c)/2 l+1 , l being the number of the level. The corresponding number of interior grid points is n l = 2 l+1 − 1 on each direction. The largest possible step corresponds to the level denoted by l = 0 on which the grid has a single point: The grid on the level l will contain the points (x i , y j ), i, j = 1, 2, . . . , n l , and will be denoted by G l . The value of the exact solution in the point (x i , y j ) is denoted by c i, j .

Finite Element Discretization.
According to [4], in order to apply the finite element discretization, some transformations of the given equation have to be made. So, the equation to be discretized is multiplied by a test Using Green's formula, the equation above becomes as follows: The functions c and v are approximated using some continu- being the number of interior points of the grid on level l), through the following relations: where Replacing these approximations in (4), the system obtained is: where Computational , the values of the integral (8) for the problem (1) are given in the following 4 × 4 matrix: Thus the differential equation (1) is approximated in a grid point (x i , y j ), i, j = 1, . . . , n l by the following system of linear equations (see [9]):

Solving Method
The systems generated in Section 3 can be written on any level l. Each system contains n 2 l unknowns. The solution is exactly computed on a level l 0 , for example, on l 0 = 2 or l 0 = 3 using Gauss elimination method. The exact solution on the level l 0 , for the problem is approximated by c i , i ∈ {1, 2, . . . , n 2 l0 } (Figure 2), wich only contains an error term due to the discretization. In order to solve problem on the level l, the grid already obtained will be further divided. Thus, each domain from the grid, Ω k , will be splitted into n i subdomains, where n i = 2 li+1 − 1 and l i = l − l 0 − 1.
On each subdomain Ω k , the discretization of the differential equation leads to a system whose matrix has the same form as the one on l 0 level. But on the level l 0 the boundary values were given in the hypothesis. For the systems on the level l to be precisely solved on Ω k , one has to determine as accurate as possible the n i values on each of the sides of the domain Ω k . Two possible ways to accomplish this are given in the following subsections. [7], the value of the approximation on level l is denoted by c (l) . On the borders of Ω k , these values are defined through the following relations: . . . , n, j = 1, . . . , n, N = n i + 1; . . . , n, j = 0, . . . , n, k = 1, . . . , n i . (13)

Stellar Prolongation.
In [9], a new type of prolongation which we called "stellar prolongation" because the nodes h a 1 a 2 a n a 2n 0 x 0 h 0 a n 2 +1 a n 2 +2 a n 2 +n a n+1 a n+2 Figure 4 involved in computation are in the shape of a star, is presented. We use this technique in what follows, too.
In order to determine more accurately the values of the solution on the borders of Ω k , instead of pondered arithmetic mean prolongation one can use the solutions of the systems obtained by discretizing the initial equation in the grid points corresponding to the values a i and b i , i = 1, 2, . . . , n 2 + n from Figures 3 and 4.
The components of the constant terms vector are in this case: (17) c f r is a function which is zero inside the domain Ω on wich the system is solved and is equal to the border values on δΩ, and h is the grid step on l 0 level. The values b k , k = 1, 2, . . . , n(n + 1) are obtained from a system whose matrix is also of the form (14), but in which: l 6 l 3 0 · · · 0 l 9 l 6 l 3 · · · 0 0 l 9 l 6 · · · 0 . . . . . .

Computational and Mathematical Methods in Medicine
In the matrix A: x = x 0 and y = x 0 for the first line of blocks in (10), on Ω A and Ω B , x = 1 and y = 1 on Ω C and Ω D , while x = 1, y = 1 for the last line of blocks on Ω A and Ω B , and x = 1 − x 0 , y = 1 − x 0 on Ω C and Ω D . For the remainder of the lines: x = 1, y = 1 (see Figure 5).
For the matrix B: x = 1/x 0 and y = 1 for the first line of blocks on Ω A and Ω D , x = 1 and y = 1 on Ω B and Ω C . The last line has: x = 1, y = 1 on Ω A and Ω D , and on Ω B and Ω C : x = 1/(1 − x 0 ), y = 1. For the other lines: x = 1, y = 1 (see Figure 6).
The system obtained by discretizing the problem on every subdomain Ω iN+ j , i = 0, . . . , n, j = 1, . . . , n has n 2 i equations and unknowns and will be solved using the Gauss full elimination method. The solutions a and b computed above will be used as boundary conditions on this domain (Figure 7).
Reuniting the solutions computed on the grid corresponding to the level l 0 and the ones from every subdomain, one gets the final solution on the work level l.

Numerical Results and Conclusions
The method presented above allows the computing of the concentration of the solvent in blood at any point of the grid for different grid steps, and can be done for various dimensions of the vessel or values of the coefficients. The importance of such a method comes from the fact that the evaluation of cutaneous circulation can be a predictive parameter for the age-related pathology.