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

Authors: Gabriela Nuţ; Ioana Chiorean; Maria Crişan
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## Abstract

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.


---

## Body

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


## 2. The Model Problem

In order to study the transport problem of a solute in a vessel, we consider the following partial derivative equation:(1)mΔc+n∇c=f,(x,y)∈Ω,c=g1,(x,y)∈Γin,c=g2,(x,y)∈Γout,c=ξg3,(x,y)∈Γ.n=(n1,n2).
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: n1(∂c/∂x), and another for the diffusion along the y-direction: n2(∂c/∂y). m represents the diffusivity of the solute, n a given velocity field, α is a reaction coefficient, f—a possible forcing term for the solution concentration due, for example, to chemical reactions.The ratio between coefficientsm and n determines the predominance of the convection or the diffusion in the physical process. If:(2)|n|m≫1,
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 104, 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.Figure 1The 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=g1, (x,y)∈Γin.(ii)
Γout, which corresponds to the point where the concentration of the solution can be measured: u=g2, (x,y)∈Γout.(iii)
Γ: in these points the solution is in contact with the vessel walls, ξ represents the wall permeability, and g3 is a function that decreases from g1 to g2.


## 3. 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–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 stephx on the Ox-direction and hy 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 l0, 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 l0 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 hlx=(b-a)/2l+1 and hly=(d-c)/2l+1, l being the number of the level. The corresponding number of interior grid points is nl=2l+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: ((b-a)/2,(d-c)/2). The grid on the level l will contain the points (xi,yj),i,j=1,2,…,nl, and will be denoted by Gl. The value of the exact solution in the point (xi,yj) is denoted by ci,j.


### 3.1. 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 function v∈U(Ω)={v∣v∈H1(Ω),v=g on Γ}, then is integrated on the domain Ω; H1(Ω)⊂C0 is the space of functions with square integrable derivatives on Ω,(3)-m∬ΩΔcvdxdy+∬Ωn∇cvdxdy+α∬Ωcvdxdy=∬Ωfvdxdy.
Using Green’s formula, the equation above becomes as follows:(4)m∬Ω∇c∇vdxdy-m∫δΩ∂c∂nvds+∬Ωn∇cvdxdy+α∬Ωcvdxdy=∬Ωfvdxdy,u,v∈U(Ω).
The functions c and v are approximated using some continuous functions, Φi (Φi(xj,yj)=δij,i,j=1,…,N, N=nl2 being the number of interior points of the grid on level l), through the following relations:(5)c≈∑i=1NuiΦi,v≈∑j=1NvjΦj,
where ci=c(xi,yi),i=1,…,N. Replacing these approximations in (4), the system obtained is:(6)∑j=1NKijcj=Fi,i=1,…,N,
where(7)Kij=∬Ω[m(∂Φi∂x∂Φj∂x+∂Φi∂y∂Φj∂y)+n1Φi∂Φj∂x+n2Φi∂Φj∂y+αΦiΦj]dxdy,Fi=∬ΩfΦidxdy.The restrictions ofK and F on a domain ΩA=[a,b]×[c,d] are as follows:(8)kijA=m∬ΩA(∂ΨiA∂x∂ΨjA∂x+∂ΨiA∂y∂ΨjA∂y)dxdy+n1∬ΩAΨiA∂ΨjA∂xdxdy+n2∬ΩAΨiA∂ΨjA∂ydxdy+α∬ΩAΨiAΨjAdxdy,i,j=1,…,4,(9)fkA=∬ΩAf(x,y)ΨkA(x,y)dxdy,k=1,…,4.
If x=(d-c)/(b-a), the values of the integral (8) for the problem (1) are given in the following 4 × 4 matrix:(10)k=(kij)i,j=1:4¯=m6[2x+2x-2x+1x-x-1xx-2x-2x+1x2x+2xx-2x-x-1x-x-1xx-2x2x+2x-2x+1xx-2x-x-1x-2x+1x2x+2x]+n1d-c12[-221-1-221-1-112-2-112-2]+n2b-a12[-2-112-1-221-1-221-2-112]+α(b-a)(d-c)36[4212242112422124].Thus the differential equation (1) is approximated in a grid point (xi,yj),i,j=1,…,nl by the following system of linear equations (see [9]):(11)[[k24Dk23D+k14Ck13Ck34A+k21Dk33A+k44B+k11C+k22Dk43B+k12Ck31Ak32A+k41Bk42B]]ci,j=f3A+f4B+f1C+f2D,
where(12)[[abcdefghk]]ui,j=aui-1,j+1+bui,j+1+cui+1,j+1+dui-1,j+eui,j+fui+1,j+gui-1,j-1+hui,j-1+kui+1,j-1.


## 3.1. 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 function v∈U(Ω)={v∣v∈H1(Ω),v=g on Γ}, then is integrated on the domain Ω; H1(Ω)⊂C0 is the space of functions with square integrable derivatives on Ω,(3)-m∬ΩΔcvdxdy+∬Ωn∇cvdxdy+α∬Ωcvdxdy=∬Ωfvdxdy.
Using Green’s formula, the equation above becomes as follows:(4)m∬Ω∇c∇vdxdy-m∫δΩ∂c∂nvds+∬Ωn∇cvdxdy+α∬Ωcvdxdy=∬Ωfvdxdy,u,v∈U(Ω).
The functions c and v are approximated using some continuous functions, Φi (Φi(xj,yj)=δij,i,j=1,…,N, N=nl2 being the number of interior points of the grid on level l), through the following relations:(5)c≈∑i=1NuiΦi,v≈∑j=1NvjΦj,
where ci=c(xi,yi),i=1,…,N. Replacing these approximations in (4), the system obtained is:(6)∑j=1NKijcj=Fi,i=1,…,N,
where(7)Kij=∬Ω[m(∂Φi∂x∂Φj∂x+∂Φi∂y∂Φj∂y)+n1Φi∂Φj∂x+n2Φi∂Φj∂y+αΦiΦj]dxdy,Fi=∬ΩfΦidxdy.The restrictions ofK and F on a domain ΩA=[a,b]×[c,d] are as follows:(8)kijA=m∬ΩA(∂ΨiA∂x∂ΨjA∂x+∂ΨiA∂y∂ΨjA∂y)dxdy+n1∬ΩAΨiA∂ΨjA∂xdxdy+n2∬ΩAΨiA∂ΨjA∂ydxdy+α∬ΩAΨiAΨjAdxdy,i,j=1,…,4,(9)fkA=∬ΩAf(x,y)ΨkA(x,y)dxdy,k=1,…,4.
If x=(d-c)/(b-a), the values of the integral (8) for the problem (1) are given in the following 4 × 4 matrix:(10)k=(kij)i,j=1:4¯=m6[2x+2x-2x+1x-x-1xx-2x-2x+1x2x+2xx-2x-x-1x-x-1xx-2x2x+2x-2x+1xx-2x-x-1x-2x+1x2x+2x]+n1d-c12[-221-1-221-1-112-2-112-2]+n2b-a12[-2-112-1-221-1-221-2-112]+α(b-a)(d-c)36[4212242112422124].Thus the differential equation (1) is approximated in a grid point (xi,yj),i,j=1,…,nl by the following system of linear equations (see [9]):(11)[[k24Dk23D+k14Ck13Ck34A+k21Dk33A+k44B+k11C+k22Dk43B+k12Ck31Ak32A+k41Bk42B]]ci,j=f3A+f4B+f1C+f2D,
where(12)[[abcdefghk]]ui,j=aui-1,j+1+bui,j+1+cui+1,j+1+dui-1,j+eui,j+fui+1,j+gui-1,j-1+hui,j-1+kui+1,j-1.


## 4. Solving Method

The systems generated in Section3 can be written on any level l. Each system contains nl2 unknowns. The solution is exactly computed on a level l0, for example, on l0=2 or l0=3 using Gauss elimination method.The exact solution on the levell0, for the problem is approximated by ci, i∈{1,2,…,nl02} (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 ni subdomains, where ni=2li+1-1 and li=l-l0-1.Figure 2On each subdomainΩk, the discretization of the differential equation leads to a system whose matrix has the same form as the one on l0 level. But on the level l0 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 ni values on each of the sides of the domain Ωk. Two possible ways to accomplish this are given in the following subsections.


### 4.1. Pondered Arithmetic Mean Prolongation

As in [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:(13)cjN+1,iN+1(l)=c(i-1)n+j(l0),cjN+1,iN+1+k(l)=1N(kcin+j(l0)+(N-k)c(i-1)n+j(l0)),i=0,…,n,j=1,…,n,N=ni+1;cjN+1+k,iN+1(l)=1N(kc(i-1)n+j+1(l0)+(N-k)c(i-1)n+j(l0)),i=1,…,n,j=0,…,n,k=1,…,ni.


### 4.2. Stellar Prolongation

In [9], a new type of prolongation which we called “stellar prolongation” because the nodes 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 ai and bi, i=1,2,…,n2+n from Figures 3 and 4.Figure 3Figure 4The valuesak, k=1,2,…,n(n+1) are computed by solving a system with the following matrix:(14)A=[CDΘ⋯ΘΘSCD⋯ΘΘΘSC⋯ΘΘ⋮⋱ΘΘΘ⋯SC].If the discretization is made by the finite element method, using the notation:(15)[[l1l2l3l4l5l6l7l8l9]]=[[k24Dk23D+k14Ck13Ck34A+k21Dk33A+k44B+k11C+k22Dk43B+k12Ck31Ak32A+k41Bk42B]],
where kij is given by (10), the matrix A has:(16)C=[l5l60⋯0l4l5l6⋯00l4l5⋯0⋮⋱000⋯l5],D=[l2l30⋯0l1l2l3⋯00l1l2⋯0⋮⋱000⋯l2],S=[l8l90⋯0l7l8l9⋯00l7l8⋯0⋮⋱000⋯l8].The components of the constant terms vector are in this case:(17)zin+j=f3A+f4B+f1C+fD-[[l1l2l3l4l5l6l7l8l9]]cfr(jh,(i+x0)h),i=0,…,n,j=1,…,n.cfr 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 l0 level.The valuesbk, k=1,2,…,n(n+1) are obtained from a system whose matrix is also of the form (14), but in which:(18)C=[l5l20⋯0l8l5l2⋯00l8l5⋯0⋮⋱000⋯l5],D=[l6l30⋯0l9l6l3⋯00l9l6⋯0⋮⋱000⋯l6],S=[l4l10⋯0l7l4l1⋯00l7l4⋯0⋮⋱000⋯l4].The components of the constant terms vector are now(19)zin+j=f3A+f4B+f1C+f2D-[[l1l2l3l4l5l6l7l8l9]]cfr((i+x0)h,jh),i=0,…,n,j=1,…,n.In the matrixA: x=x0 and y=x0 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-x0, y=1-x0 on ΩC and ΩD. For the remainder of the lines: x=1, y=1 (see Figure 5).Figure 5For the matrixB: x=1/x0 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-x0), y=1. For the other lines: x=1, y=1 (see Figure 6).Figure 6The system obtained by discretizing the problem on every subdomainΩiN+j, i=0,…,n, j=1,…,n has ni2 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).Figure 7Reuniting the solutions computed on the grid corresponding to the levell0 and the ones from every subdomain, one gets the final solution on the work level l.


## 4.1. Pondered Arithmetic Mean Prolongation

As in [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:(13)cjN+1,iN+1(l)=c(i-1)n+j(l0),cjN+1,iN+1+k(l)=1N(kcin+j(l0)+(N-k)c(i-1)n+j(l0)),i=0,…,n,j=1,…,n,N=ni+1;cjN+1+k,iN+1(l)=1N(kc(i-1)n+j+1(l0)+(N-k)c(i-1)n+j(l0)),i=1,…,n,j=0,…,n,k=1,…,ni.


## 4.2. Stellar Prolongation

In [9], a new type of prolongation which we called “stellar prolongation” because the nodes 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 ai and bi, i=1,2,…,n2+n from Figures 3 and 4.Figure 3Figure 4The valuesak, k=1,2,…,n(n+1) are computed by solving a system with the following matrix:(14)A=[CDΘ⋯ΘΘSCD⋯ΘΘΘSC⋯ΘΘ⋮⋱ΘΘΘ⋯SC].If the discretization is made by the finite element method, using the notation:(15)[[l1l2l3l4l5l6l7l8l9]]=[[k24Dk23D+k14Ck13Ck34A+k21Dk33A+k44B+k11C+k22Dk43B+k12Ck31Ak32A+k41Bk42B]],
where kij is given by (10), the matrix A has:(16)C=[l5l60⋯0l4l5l6⋯00l4l5⋯0⋮⋱000⋯l5],D=[l2l30⋯0l1l2l3⋯00l1l2⋯0⋮⋱000⋯l2],S=[l8l90⋯0l7l8l9⋯00l7l8⋯0⋮⋱000⋯l8].The components of the constant terms vector are in this case:(17)zin+j=f3A+f4B+f1C+fD-[[l1l2l3l4l5l6l7l8l9]]cfr(jh,(i+x0)h),i=0,…,n,j=1,…,n.cfr 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 l0 level.The valuesbk, k=1,2,…,n(n+1) are obtained from a system whose matrix is also of the form (14), but in which:(18)C=[l5l20⋯0l8l5l2⋯00l8l5⋯0⋮⋱000⋯l5],D=[l6l30⋯0l9l6l3⋯00l9l6⋯0⋮⋱000⋯l6],S=[l4l10⋯0l7l4l1⋯00l7l4⋯0⋮⋱000⋯l4].The components of the constant terms vector are now(19)zin+j=f3A+f4B+f1C+f2D-[[l1l2l3l4l5l6l7l8l9]]cfr((i+x0)h,jh),i=0,…,n,j=1,…,n.In the matrixA: x=x0 and y=x0 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-x0, y=1-x0 on ΩC and ΩD. For the remainder of the lines: x=1, y=1 (see Figure 5).Figure 5For the matrixB: x=1/x0 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-x0), y=1. For the other lines: x=1, y=1 (see Figure 6).Figure 6The system obtained by discretizing the problem on every subdomainΩiN+j, i=0,…,n, j=1,…,n has ni2 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).Figure 7Reuniting the solutions computed on the grid corresponding to the levell0 and the ones from every subdomain, one gets the final solution on the work level l.


## 5. 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.The values of the parameters used here are as follows:(i)
the diffusion coefficient:m=10-3cm2s-1;(ii)
the velocity vector:n=(10,10)cms-1;(iii)
the wall permeability:ξ=1cms-1;(iv)
the concentration onΓin: g1=2·10-3;(v)
the concentration onΓout: g1=1.98·10-3;(vi)
the forcing term for the solute concentrationf=0;(vii)
the vessel dimensions have been chosen to be1cm×0.2cm.For the convection-diffusion problem (1), the results are presented in Figures 8 and 9.Figure 8
The concentration computed on the domain[0,1]×[0,0.2].Figure 9
The concentration map on the domain[0,1]×[0,0.2].



---

*Source: 204926-2012-01-26.xml*