# Dual reciprocity boundary element formulation applied to non-linear Darcian diffusive-adptive problems

Carlos Loeffler, Felipe Carlos Neves, Patrícia Oliveira

<figure><img src="image_1.png" /></figure>

## To cite this version:

Carlos Loeffler, Felipe Carlos Neves, Patrícia Oliveira. Dual reciprocity boundary element formulation applied to non-linear Darcian diffusive-adptive problems. 11th International Conference on Advances in Boundary Element Techniques, Jul 2010, Berlin, Germany. hal-01580915

HAL Id: hal-01580915

https://hal.science/hal-01580915

Submitted on 3 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dual Reciprocity Boundary Element Formulation Applied to Non-Linear Darcian Diffusive-Advective Problems

C.F. Loeffler¹, F. P. Neves² and P.C. Oliveira³

¹ DEM/PPGEM, Universidade Federal do Espírito Santo, Vitória, ES – Brasil
carlosloeffler@bol.com.br

² Corpo de Bombeiros Militar do Espírito Santo, Vitória, ES – Brasil
patriciofire@gmail.com

³ Centro Agropecuário, Universidade Federal do Espírito Santo, Alegre, ES - Brasil
pacololiveira@yahoo.com.br

**Keywords:** Diffusive-Advective Problems, Dual Reciprocity Boundary Element Formulation.

**Abstract.** This paper presents a Dual Reciprocity boundary element model to simulate computationally the distribution of pressures, velocities, temperatures and heat fluxes in two-dimensional control volumes for steady state condition. Darcy's law is assumed to represent the Momentum Equation, making the link between pressure and velocity fields, mathematically given by Laplace's Equation for the linear case. The nonlinear analysis is characterized by the dependence of the velocity field with temperature, resulting in a model represented by Poisson's Equation. Velocity results are then implemented in the Diffusive-Advective Equation, generating temperatures and heat fluxes. These values are introduced in Poisson's Equation to generate new pressure and velocity fields. The process is repeated until the results converge to stable values.

## Introduction

Engineering problems involving mass or energy transport through diffusion associated with advection are very common. The most well known cases concerning to heat transfer on a fluid boundary layer in laminar flows, typical of heat exchanger fins and wings, and incompressible fluid transport with low viscosity in pipes subject to temperature gradients. There are other situations of no lesser importance, such as the surface dispersion of pollutants or mixtures in a homogeneous aqueous medium, drying of agricultural products, absorption of liquids in unsaturated region, oil extraction in porous rock and aeration problems in closed environments. In the general case, the complexity of physical phenomena means a very difficult mathematical model, consisting of nonlinear partial differential equations, which necessarily requires the use of numerical methods for its approximate solution. The general treatment of these cases is usually substituted for the solution of simplified models, which allow achieving a more accessible solution while at same time attaining the practical requirements.

## Basic Equations

The physical domain consists of a two-dimensional control volume Ω (x, y), where exists a field of velocities and pressures. Along the boundaries Γ (x, y) are prescribed temperatures (essential conditions) or fluxes (natural conditions). Fig. 1 illustrates these features. Therefore, considering homogeneous, isotropic and incompressible fluid flow, the continuity equation [5] is expressed by:

$$ \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} = 0 \qquad (1) $$

In the previous equation, $v_x$ and $v_y$ are the components of the flow velocity vector, while x and y are the spatial coordinates. The Energy equation [5] corresponds to the Diffusive-Advective equation, given by:

$$ \rho C_v \left[ \frac{\partial T}{\partial x} v_x + \frac{\partial T}{\partial y} v_y \right] = K \left[ \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right] \qquad (2) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Two-dimensional control volume with a velocity field given by ($v_x, v_y$)</figcaption></figure>

Where $C_v$ is the specific heat at constant volume and K is the thermal conductivity. The Momentum Equation [2], assuming the presence of body forces $b_x$ and $b_y$, it is written as:

$$ \rho\left[\frac{\partial\vec{V}}{\partial x}v_x + \frac{\partial\vec{V}}{\partial y}v_y\right] = -\vec{V}\rho + \mu\nabla^2\vec{V} + \rho\vec{b} \quad (3) $$

Where $\rho$ is the density, $\rho$ is the pressure field in the fluid and $\mu$ is the viscosity. The total derivative present on the left hand side of eq (3) is proportional to the force required to overcome the inertia of the fluid element. This term can be omitted from the equation when viscous flow is assumed, i.e., Reynolds number is lower [2]. Under these conditions, considering incompressible fluid flow, the former equation is simplified [3]:

$$ \vec{V}\rho = \rho\vec{b} + \mu\nabla^2\vec{V} \quad (4) $$

This hypothesis applies to the description of the particles with negligible mass flowing through the fluid medium, such as settlement and sedimentation. This assumption is also relevant in porous media problems.

### Darcian Model for Momentum Equation

Since the porosity is relatively low in the continuous middle, Darcy's law can be applied to simplify the Momentum Equation [6, 7]. It is valid only for slow viscous flow. Thus, the velocity components are given by:

$$ \frac{\mu\Phi}{\chi}\vec{V} = \rho\vec{b} - \vec{V}\rho \quad (5) $$

Where $\chi$ is the permeability of the medium and $\Phi$ is the porosity. Therefore, for laminar flow conditions, comparing the value of the velocity field given by Darcy's law and the expression correlated given by the Momentum Equation (eq 4) for incompressible viscous flow, it follows that:

$$ \nabla^2 \vec{V} = -\frac{\Phi}{\chi} \vec{V} \quad (6) $$

### Boussinesq's Approximation

For saturated porous rock, the Boussinesq approximation [5] can be used to describe density changes due to the temperature effect. This approach provides that where the gravitational force acts, changes in the density are not small enough to be neglected. Expressing the density as a function of temperature, it follows that:

$$ \rho_F = \rho_O[1 - \beta(T_F - T_O)] \quad (7) $$

Where $\rho_F$ is the density at a general temperature $T_F$, while $\rho_O$ is the density at the reference temperature $T_O$ and $\beta$ is the thermal expansion coefficient of the porous medium. Considering now the density variation:
$$ \frac{\mu\Phi}{\chi} \vec{V} = \rho_o[1 - \beta(T_F - T_o)]\vec{b} - \vec{V}p \qquad (8) $$

Putting the previous equation on the dimensionless form, its components can be rewritten as:

$$ v_X^* = R_A^X T^* - \frac{\partial p^*}{\partial x^*}, \quad v_Y^* = R_A^Y T^* - \frac{\partial p^*}{\partial y^*} \qquad (9) $$

Where $R_A$ is the Rayleigh Number, defined in each direction x and y as a function of the body force $b_x$ or $b_y$. Usually, only the gravitational force is computed, so that the other component is null. In these conditions, substituting eq (9) into the continuity equation, the pressure field is given by the well known Poisson's Equation:

$$ \frac{\partial^2 p^*}{\partial x^2} + \frac{\partial^2 p^*}{\partial y^2} = -R_A^x \frac{\partial T^*}{\partial x^*} - R_A^y \frac{\partial T^*}{\partial y^*} = q(x, y) \qquad (10) $$

## Solution of the Coupled Problem

Due to the Boussinesq approximation the Diffusive-Advective Equations is coupled to the Poisson equation, since both depend on the temperature. This means that the whole problem is not linear. Here, an iterative process is used to solve this engagement, as well as to start the process. It takes the right hand side of eq (10) as a source term, in which the initial temperatures and fluxes are identified according to a given reference. Solving this equation, the field of pressures and velocities necessary to comprise the Diffusive-Advective Equation is determinate. Achieved the new temperatures and fluxes, the source term previously arbitrated is readjusted and the Poisson Equation is solved again, giving new values of pressures and velocities. The process carries on until the magnitudes no longer significantly change between two successive iteration steps.

## BEM Formulation for the Momentum Equation

The traditional mathematical treatment with the BEM to the left hand side of the Momentum Equation is well known in the literature [1], using an auxiliary function $u^*(\xi;x,y)$, called the fundamental solution. Being a self-adjoint operator, the Laplacian allows the application of the integration by parts scheme twice into the kernel of integral equation, so that it can rewrite as:

$$ C(\xi)p(\xi) + \int_{\Gamma} \{pq^*\}d\Gamma - \int_{\Gamma} \{vu^*\}d\Gamma = \int_{\Omega} qu^*d\Omega \qquad (11) $$

The coefficient $C(\xi)$ is related to the position of the source point $\xi$ on the physical domain $\Omega(x,y)$. If the source point is located along the smooth boundaries, $C(\xi)=1/2$. In this last equation appear:

$$ q^* = \frac{\partial u^*}{\partial x} \frac{\partial x}{\partial n} + \frac{\partial u^*}{\partial y} \frac{\partial y}{\partial n} = \frac{\partial u^*}{\partial n}, \quad v = \frac{\partial p}{\partial x} \frac{\partial x}{\partial n} + \frac{\partial p}{\partial y} \frac{\partial y}{\partial n} = \frac{\partial p}{\partial n} \qquad (12) $$

The domain source term is approached by the Dual Reciprocity formulation (DRF) [3], which uses suitable interpolation radial basis functions to describe the domain term and transform the associated domain integrals onto boundary integrals. This is done expressing the domain source term $q(x,y)$ through a finite linear combination of new auxiliary functions $\alpha^j$ and $\beta^j$ where $\alpha^j$ and $\beta^j$ are the unknown coefficients and the $F^j$ are the mentioned interpolation functions. The approach is constructed considering basis points $X^j$ selected in both the boundary and inside the domain. Thus, according to the DRF approach, the right hand side of eq (11) is rewritten as:

$$ \int_{\Omega} q u^* d\Omega = \int_{\Omega} \alpha^j \Psi_{ii}^j u^* d\Omega = \alpha^j \left[ \int_{\Gamma} \eta^j u^* d\Gamma - \int_{\Gamma} \Psi^j q^* d\Gamma - C(\xi) \Psi^j(\xi) \right] \qquad (13) $$

The new function $\eta^j$ represents the normal spatial derivative of the potential $\psi^j$. The choice of the interpolation functions is object of many researches, but the radial basis functions are the most interesting class of functions due to its flexibility, invariance, among other properties.
## BEM Formulation for the Energy Equation

The Energy Equation is also solved through to DRF, which does not impose constrains on the velocity field and uses the diffusive fundamental solution $u*(\xi;x,y)$ to comprise the integral boundary element equation. According to similar procedure foregoing presented, the advective term is approached as a field source $g(x, y)$, that is:

$$K \iint_{\Omega} \left( v_x \frac{\partial T}{\partial x} + v_y \frac{\partial T}{\partial y} \right) u^* d\Omega = \rho C_v \iint_{\Omega} g(x, y) u^* d\Omega = \rho C_v \gamma^j \left\{ \iint_{\Gamma} \left[ \eta^j(X^j; X) u^*(\xi; X) - \Psi^j(X^j; X) q^*(\xi; X) \right] d\Gamma - C(\xi) \Psi^j(\xi) \right\} \quad (14)$$

## Discretization and Matricial BEM Equation

The BEM discretization procedure is well known and very simple: traditional shape functions are used to approximate the basic variable and its normal derivative. Using the DRF, the auxiliary functions are also approximated by shape functions along the boundary elements. In the case of Poisson's equation, the left side generates the traditional matrices H and G for Laplace's problem [1]; on the right side, due to the reciprocity of the integral components, it follows that:

$$H\rho - G\nu = \{H\Psi - G\eta\}\alpha = \{H\Psi - G\eta\}F^{-1}q = Q \quad (15)$$

For the Diffusive-Advective equation the procedure is similar, but the replacement of the vector of coefficients $\alpha$ is more elaborated:

$$HT - Gq = \frac{\rho C_v}{K} \{H\Psi - G\eta\}\gamma = \frac{\rho C_v}{K} \{H\Psi - G\eta\}F^{-1}g = \frac{\rho C_v}{K} [H\Psi - G\eta]F^{-1}V_i T_{,i} \quad (16)$$

The temperature spatial derivatives can be eliminated easily by interpolation, such that the definitive matrix equation becomes:

$$HT - Gq = \frac{\rho C_v}{K} [H\Psi - GN] F^{-1}V_i F_{,i} F^{-1}T = RT \rightarrow (HT - RT) = Gq \quad (17)$$

It must be pointed out that the velocity field must be completely known at the boundaries of the control volume, whereas the solution of the Poisson Equation gives only the values of normal velocities. The tangential velocities need to be calculated to comprise the complete velocity field in the Diffusion-Advection Equation. These velocities can be obtained using a simple radial basis function interpolation, that is:

$$p_{,i} = v_{,i} = F_{,i} \beta = F_{,i} F^{-1} p \quad (18)$$

The same requirement concerns to the tangential derivatives of temperature, which comprising the source term in Poisson's Equation for the iterative process. Also, values of fluxes in internal points need to be calculated, especially in order to improve the accuracy of the FDR procedure in the solution of Poisson's Equation.

## Example

<figure><img src="image_4.png" /><figcaption>Figure 2: (a) geometrical characteristics; (b) pressure boundary conditions; (c) thermal boundary conditions; (d) Specification to the different faces of the silo.</figcaption></figure>
The example chosen is a vertical silo, where fields of pressure, velocities and temperature are applied as shown in Fig. 2. The mesh consists of 40 constant boundary elements, with the same size. The numeric order of the elements obeys the trigonometric sense along the boundary. For a better understanding of the results, the control volume that makes up the silo is divided into four sections, also pointed out at Fig. 2.

As previously discussed, the problem is non-linear. Two distinct parts comprise the complete solution, in which the velocity field is obtained initially and then determines the distribution of temperatures and fluxes. New values of fluxes comprise the source term of Poisson's Equation. For the proposed problem the prescribed values of temperature are dimensionless, being the reference value taken equal to 100°C. The Rayleigh number considered is equal to 6.

The plots in Fig. 3 show the values of normal component of velocities calculated on face 3 for the three iterations carried out by solving the momentum equation. Comparing the first with the second iteration, it is significant the difference in velocity results calculated, due to the introduction of gravitational effect. However, a greater number of iterations were not necessary for the value of Ra adopted, since no meaningful difference appears for the third iteration results when compared to the second. Also on Fig. 3, it can be noted that the presence of the gravitational effect reduces the intensity of the normal velocities along face 3, such as expected.

<figure><img src="image_5.png" /><figcaption>Figure 3: Normal velocities at the nodal points located along face 3.</figcaption></figure>

The gravitational effect also reduces the pressure intensity on the bottom, i.e., on face 1, as shown in Fig. 4. The pressure drop along the horizontal direction is due to the remoteness of the higher values of pressure, located on face 4.

<figure><img src="image_6.png" /><figcaption>Figure 4: Pressure results at the nodal points located on face 1.</figcaption></figure>
Unlike the calculation of pressures and normal velocities, the solution for temperatures and fluxes practically was not affected by the gravitational action, as shown in Fig. 5, where the results for the three iterations are practically the same.

<figure><img src="image_7.png" /><figcaption>Figure 5: Diffusive flux results at the nodal points located on face 3.</figcaption></figure>

## Conclusions

Many preliminary numerical experiments made with linear models showed that the BEM applies with necessary accuracy in problems governed separately by Laplace, Poisson and Diffusive-Advective Equations [4], this latter considering cases in which the Peclet number is not high, such as when the diffusive phenomenon is not negligible. Thus, it was expected that the method could be applied successfully in cases where nonlinearities were present. One of these forms of nonlinearity is the inclusion of the gravitational action through the coupling between the diffusive-advective model and a pressure field. This is a problem that typically appears in cases of porous media.

In this more complex problem, the application of the Dual Reciprocity BEM was satisfactory. Thus achieving quick and easy results for the various state variables comprise the proposed physical problem, such as the values of normal and tangential velocity, pressures, temperatures and heat flows, normal and tangential, or in the walls of the silo and inside. For the sake of space, only a small fraction of these can be presented here.
For the numerical model proposed, the solution presented very low computational cost. It must be highlighted the easiness of data entry and interactivity between the Poisson and Diffusive-Advective problems.

## References

[1] C.A. Brebbia, *The Boundary Element Method for Engineers*, Pentech Press, London (1978).

[2] M.N. Özisik, *Boundary Value Problems of Heat Conduction*, Dover Publ., New York, (1968).

[3] P.W. Partridge, C.A. Brebbia, L.C. Wrobel, *The Dual Reciprocity Boundary Element Method*, Computational Mechanics Publications and Elsevier, London, 1992.

[4] P.A. Ramachandran, *Boundary Element Methods in Transport Phenomena*, Computational Mechanics Publication and Elsevier Applied Science, London (1994).

[5] F.M. White, *Fluid Mechanics*, McGraw-Hill Int., Singapore (1986).

[6] Zhao, C., Hobbs, B.E., Muhlhaus, H.B., *Finite Element Modelling of Temperature Gradient Driven Rock Alteration and Mineralization in Porous Rock Masses*, Int. J. Numer. Anal. Meth. Geomech., **21**, 863-881 (1997).

[7] Zhao, C., Hobbs, B.E., Hornby, P., Ord. A., Peng, S., *Numerical Modelling of Fluids Mixing, Heat Transfer and Non-equilibrium Redox Chemical Reactions in Fluid-saturated Porous Rocks*, Int. Journal for Num. Methods in Engineering, **66**, 1061-1076 (2006).