IMPLEMENTATION OF DRBEM FOR THE DETERMINATION OF THE HEAT FLUX IN AN INVERSE PROBLEM

A numerical investigation of inverse unsteady natural convection ow in a square cavity lled with Cu water nanouid is performed. In the direct problem, the enclosure is bounded by one isothermally heated vertical wall at temperature Tm and by three adiabatic walls. In the inverse problem, the enclosure is bounded by right hostile wall on which no boundary condition can be prescribed or measured and by left accessible wall on which both the boundary temperature and heat ux data are overspecied. The dual reciprocity boundary element method (DRBEM) with the fundamental solutions of Laplace and modied Helmholtz equations is used for the solutions of direct and inverse problems. Inhomogeneities are approximated with radial basis functions. Computations are performed for several values of Rayleigh number (Ra), solid volume fraction ( ) and percentage of noise ( ), and accurate and stable results are given for three forms of heat ux namely, steady heat ux (q = q(y)), time dependent uniform heat ux (q = q(t)) and non-uniform time dependent heat ux (q = q(y; t)).


Introduction
In the thermo- ‡uid areas, ‡ow heat transfer analysis is one of the most interesting subject among the researchers. The researchers who study in this area focused to the idea of using nano ‡uid instead of a using conventional heat transfer base ‡uids with the low thermal conductivity such as water, engine oil etc. Nano ‡uid is the combination of the base ‡uid and added nanoparticles, and this term is made use of for the …rst time by Choi in [1]. The aim of using nano ‡uids is to increase the thermal performance of conventional ‡uids. Some studies published as review on nano ‡uids can be found in [2]- [9].
In some applications of heat transfer, because of the high temperature it may not be possible to measure of boundary conditions on some part of the surface. The nonmeasurable temperature or heat ‡ux boundary conditions can be determined by the temperature measurements from the other parts of the surface or inside of the body. The estimation of temparature or heat ‡ux boundary conditions is called inverse heat conduction problem for temperature and the excellent study about this subject has been done by Beck et al. [10]. After that, a number of studies for the inverse problem have been published in [11]- [16]. An inverse forced-convection problem in a channel has been investigated in [11] while convection-di¤usion equation has been solved in [12]- [14]. In addition, advection-di¤usion and reaction-di¤usion equations for an inverse problem case have been considered in [15] and [16], respectively. Later on, one-dimensional inverse heat conduction problem has solved in [17].
The inverse natural convection ‡ow has been solved by using a solution algorithm based on the sequential function speci…cation method in [18], and by using an implicit control volume discretization procedure in [19]. Also, DRBEM applications have been presented for the inverse natural convection ‡ow under a magnetic …eld problem in [20] and for the Magnetohydrodynamics (MHD) ‡ow problem in [21]. In all these studies mentioned above, the governing equations are coupled nonlinearly which leads to extra di¢ culty for the inverse problem.
The purpose of this study is to solve the inverse natural convection of copper based nano ‡uid problem using di¤erent type heat ‡ux boundary conditions. To do this aim, computations are carried using for three types of adiabatic boundary conditions, namely, steady heat ‡ux (q(y) = cos( y)), time dependent uniform heat ‡ux (q(t) = sin( t=t f )) and non-uniform time dependent heat ‡ux (q(y; t) = sin( t) cos( y)). Since the analytical solution of corresponding direct problem is not available, it is necessary to solve the corresponding direct problem using any numerical method. In this study both direct and inverse problems are solved using DRBEM. In the DRBEM procedure, it is only necessary to discretize the boundary of the domain to be analyzed, so it requires less storage and less computing time.
The fundamental concepts about this method can be found in [22]. Since DRBEM uses the fundamental solution of di¤erential equation for converting the domain integral into a boundary integral, the original di¤erential equation should be written as a di¤erential equation whose fundamental solution can be obtained. In this current work, vorticity transport and temperature equations are written in the forms of modi…ed Helmholtz equations by considering all the nonlinear terms as the nonhomogeneous terms.
A well-posed problem in generally can be de…ned as a mathematical problem where its solution(s) has/have three properties namely existence, uniqueness and stability. If any of these three properties mentioned above fail then the problem becomes an ill-posed problem and this kind of a problem is called an inverse problem. Due to the ill-posedness property of the problem, the small errors occurred in the measurement causes the big deviations in the inverse solution. Therefore, Tikhonov regularization method [23] have been utilized to regularize ill-posedness of the problem.
The inverse natural convection in a nano ‡uid …lled cavity presented in this study has not been considered and investigated yet. This article is planned in the following way. In Section 2 de…nition of the direct and inverse problems are given. Section 3 contains the mathematical model of the governing equations. DRBEM is described for both direct and inverse problems in Section 4. Results based on numerical simulation are discussed and presented for three types of heat ‡ux boundary conditions in Section 5. Finally, Section 6 presents the conclusions obtained from this study.

Definition of the Direct and Inverse Problems
Consider the two-dimensional unsteady equations of motion and energy for nano ‡uid (Cu-water) in the non-dimensional velocity (u; v), pressure (p) and temperature (T ) form as @u @x + @v @y = 0; (1) @v @t + u @v @x + v @v @y = @p @y + nf nf f 400 N. ALSOY AKGÜN where ; ; ; and are dynamic viscosity, e¤ective density, thermal di¤usivity and thermal expansion coe¢ cient, respectively. The above non-dimensional equations are obtained by using the following non-dimensional variables where L 0 ; g; and T = T hot T cold are the characteristic length, gravitional acceleration, kinematic viscosity and temperature di¤erence, respectively. Also, Ra is the Rayleigh number and P r is the Prandtl number. The prime in Equation (5) refers that the variables are dimensional. The thermo-physical properties of the nano ‡uid, water and copper can be given as in [24] by (6) where is nanoparticle volume fraction, k is thermal conductivity, ( C p ) nf is heat capacitance, ( ) nf is thermal expansion coe¢ cient, is dynamic viscosity and the subscripts 'nf ', 'f 'and 's'refer to nano ‡uid, ‡uid and solid, respectively. Two-dimensional unsteady equations of motion and energy for nano ‡uid in the non-dimensional stream function ( ), vorticity (w) and temperature (T ) form can be de…ned as nf nf f where u = @ @y ; v = @ @x ; w = @v @x @u @y : 2.1. Direct Problem: The schematic of the physical situations of both direct and inverse problems, which have the same boundary conditions as those previously investigated in [19], are shown in Figure (1). The boundary conditions for the direct problem are @T @x (1; y; t) = q(y; t); @T @x (0; y; t) = @T @y (x; 0; t) = @T @y (x; 1; t) = 0: where q(y,t) is the most general form of nondimensional heat ‡ux boundary condition at x = 1: Vorticity boundary conditions are obtained using the Equation (10). Direct problem is subject to the initial conditions w(x; y; 0) = T (x; y; 0) = 0:

Inverse Problem:
In the inverse problem one of the terms (boundary conditions at the right wall) is not known explicitly for the temperature and it should be speci…ed from a direct problem. In this study, the time varying temeperature measurement (T m ) at the left wall, which is given in Figure (1) and will be obtained by solving the direct problem, is considered as an overspeci…cation for the compensation of the missing boundary condition at the right wall. Thus the temperature boundary conditions for the inverse problem are Stream function and vorticity boundary conditions are the same as the direct problem. Initial conditions for vorticity transport and temperature equations are also the same as the direct problem.
In the inverse problem, the small measurement errors appeared in the computation of the T m can be caused the large errors. Thus, ill-posedness occurs in the inverse solution procedure. This means that the solutions T (1; y; t) and @T @x (1; y; t) can be unstable. Therefore, in order to obtain a stable numerical solution for the inverse problem, regularization methods, for instance Tikhonovs regularization should be used in the solution procedure.

Mathematical Model
For the solution of both direct and inverse problems DRBEM will be used as a numerical technique and so we need to use fundamental solutions of the Laplace and modi…ed Helmholtz equations. Stream function equation is in the form of Poisson equation and vorticity transport and energy equations are transformed to the nonhomogeneous modi…ed Helmholtz equations. To do this aim …rst the time derivatives in the vorticity transport and energy equations are approximated by using the forward …nite di¤erence approximations where w (s) = w(x; y; t s ); T (s) = T (x; y; t s ); t s = s t and t is the time step. Then vorticity and temperature variables in the Laplace terms are also approximated at the two time levels by using relaxation parameters w and T as Inserting the approximations in Equations (15) and (16) into the related Equations (8) and (9), two nonhomogeneous modi…ed Helmholtz equations are obtained for vorticity transport and energy equations

DRBEM Formulation of The Problem
The governing equations given in (17)- (19) for the direct and inverse formulations are solved by using DRBEM. DRBEM implementation of Poisson and modi…ed Helmholtz equations for the inverse problem can be done by using the same way as the direct problem. The DRBEM discretization of the inverse problem leads to an ill-conditioned linear system of equations for temperature due to overspeci…ed boundary conditions on the left wall and underspeci…ed conditions on the right wall. In order to obtain the stable results, the ill-conditioned linear system of equations for temperature is regularized with Tikhonov regularization method. The governing equations in the form of Poisson and modi…ed Helmholtz equations can be written in the compact forms as where b 1 ; b 2 ; b 3 are the right hand sides of corresponding equations (17) (21) and (22) are weighted in domain by the fun- equations as in [25]. By applying the Green's second identity to the resulting weighted residual statements, left hand side of the governing equations are transformed into boundary integral equations as in [25] c i where c i are the constants which depend on whether the source point i lies in the interior (c i = 1; i = 1; :::; K I ) or on the smooth boundary (c i = 1 2 ; i = 1; :::; K B ). The domain integrals caused by the inhomogeneous terms are transformed to the boundary integral equations using the DRBEM idea as in [22], [25]. To this aim the inhomogeneous terms are expanded as Here K B and K I are the number of the boundary and interior points, respectively.
1j ; 2j (t) and 3j (t) are initially unknown coe¢ cients. f j = 1+r and f j = r 2 log r are the coordinate functions. f j and f j are related to particular solutions c j ; and c w j ; c T j as After inserting f j and f j into the corresponding inhomogeneous terms in Equations (26)-(28), we have the Laplace and modi…ed Helmholtz operators inside the domain integrals (34) Now the DRBEM idea can be used for the domain integrals appearing on the right hand sides of Equations (32)-(34) in order to obtain the boundary integrals. (36) After discretizing the boundary using K B constant elements and taking K I interior nodes, the ultimate form of Equation (35)-(37) can be written as one global scheme in a matrix form for both boundary and interior nodes as where b Y represents c W and b T , represents w and T : Here the de…nition of r is r 2 = r 2 x + r 2 y where r x and r y are the components of r in the direction of x and y axes.
The entries of the matrices H; G; H 0 and G 0 are de…ned on the boundary elements j as in [22,26,27] where ij Kronecker delta function, j is the j-th boundary element and refers to w and T : The collocation of the inhomogeneities at K B + K I points results in the systems for k , (k = 1; 2; 3) where F and F are the coordinate matrices constructed by taking f j = 1 + r j and f j = r 2 j ln r j as columns. Substituting the 1 , 2 , 3 into equation (38), the …nal system of equations for the unknown vectors w; and T can be given as In this system all the square matrices have the size (K B + K I ) (K B + K I ) and all the vectors have size (K B + K I ) 1: First and second derivatives of the variables located in b i , (i = 1; 2; 3) can be approximated as where S represents w (s) ; (s+1) and T (s) : Also, we need r 2 w (s) and r 2 T (s) ; and they can be obtained by using the coordinate matrix F such that The boundary conditions for vorticity are determined by taking the curl of the velocity vector and with the help of the DRBEM coordinate matrix F as w = @v @x @u @y A system of linear equations for the variables ( -w-T ) can be constructed by using the equations (43) and by inserting the corresponding known boundary conditions as A 1 x 1 = y 1 ; Here, x 1 ; x 2 and x 3 are constructed by using the unknown information of variables and their normal derivatives. These system of linear equations are obtained for both direct and inverse problems. Any numerical method can be used to obtain the unknown x i ; (i = 1; 2; 3) in the solution of the direct problem. Thus, some part of the known temperature boundary conditions are prepared in order to use as a sensor in the solution of the inverse problem. In the solution procedure of the inverse problem, similar to the solution of direct problem, the classical DRBEM can be used to obtain x 1 and x 2 . However, the system A 3 x 3 = y 3 for the inverse temperature equation is ill-conditioned and it needs to use any regularization method to obtain a regularized solution. In this study, Tikhonov regularization method [23] is prefered since it is simple compared to the other regularization methods [16]. In the regularization procedure, the system A 3 x 3 = y 3 is written in the form Here to obtain stable regularized solutions positive regularization parameter should be chosen carefully.
In the measurement of the temperature T m at the right wall small random errors inherently occur and these errors cause big instability in the solution of the inverse problem. These measurement errors are taken into account by adding random noise data to the temperature T (0; y; t) = T m obtained at the right wall from the direct problem as T noise m (0; y; t) = T m (0; y; t)(1 + ); y 2 (0; 1); where is the percentage of the noise and is a random variable in the interval [ 1; 1] generated by IMSL library subroutines RNUN and SSCAL.

Numerical Results
In this study, DRBEM application is derived for the unsteady natural convection problem in a square cavity …lled with Cu water based nano ‡uid with three di¤erent Neumann type boundary conditions, which are imposed at the right wall and are steady heat ‡ux (q(y) = cos( y)), time dependent uniform heat ‡ux (q(t) = sin( t=t f )) and non-uniform time dependent heat ‡ux (q(y; t) = sin( t) cos( y)), depicted in Figure (1). The numerical results are discussed by using graphs to emphasize the e¤ect of Ra, ; ; and on the ‡uid behaviour. For the computation of both direct and inverse problem, the stopping criteria is taken as 10 5 : In order to achieve steady-state results by using a small number of iterations, the relaxation parameters are taken w = T = 0:9 [27].
At the beginning of the procedure, the forward …nite di¤erence approximation is used for both the time derivatives in the vorticity transport and in the temperature equations. Since forward …nite di¤erence scheme is an explicit method, the choice of t is an important point for the stability of the solutions. Also, in the modi…ed Helmholtz equations, t is located in the parameters w and T . So there is a close relation between the behavior of function K 0 (x) and the choice of t. Since K 0 (x) goes to in…nity as x goes to in…nity, too small t can lead to unstable solutions. In this present study, t values used in the computations changes between 0.01 to 0.1. In Figure (3) the e¤ect of the regulation parameter is presented within the range = 10 4 to = 10 7 : From this …gure it can be seen that when the regularization parameter takes the higher values than 10 5 and takes the smaller values than 10 7 ; the numerical results becomes unstable which is expected behavior for the ill-posed system.
The e¤ect of Rayleigh number on the steady heat ‡ux is presented in Figure (4). From the …gure it can be observed that the behavior of the inverse solution for heat ‡ux are similar to each other for all values of Rayleigh number. The best result is obtained when Ra = 10 3 . As the Rayleigh number decreases up to Ra = 0, the di¤erence between the exact and inverse solutions starts to grow especially near the top and bottom right corners of the domain. When Ra = 10 4 an improvement occurs in the results. But when the Rayleigh number takes the higher values, for instance Ra = 10 5 ; the inverse solutions lose the accuracy and stability. Ra = 0 For the same analysis, the isotherms are given for t = 0:5 in Figure (6). As can be noticed from the …gure, since the heat transfer mechanism is dominated by the conduction up to Ra = 10 3 , there is no considerable change over the isotherms. This means that the convection in the system is weak and it can not be e¤ective on the heat transfer in the system. Also, for all values of Rayleigh number, the isotherms obtained both direct and inverse problem are well matched with the each other.   The global heat ‡ux for the unsteady case is presented in Figure (9) which is depicted for each iteration up to t = 1:0: By comparing Figures (5) and (9), it can be seen that in the unsteady case the global heat ‡ux obtained from the direct and inverse problems are in good agreement with each other for all values of Rayleigh number. Figure (10) shows the e¤ect of the variation of Rayleigh numbers on the isotherms at t = 0:5. From the …gure, the e¤ect of convective forces on heat transfer mechanism can be seen clearly. When the system is dominated by conduction, the isothermes are parallel to the vertical walls. But as the Rayleigh number increases, Figures (11) shows the e¤ect of volume fraction for di¤erent Rayleigh numbers on the unsteady heat ‡ux. Similar to the steady case, as Rayleigh number increases from 0 to 10 2 ; there is no remarkable e¤ect of changes of volume fraction on the heat ‡ux. However after the Rayleigh number reaches 10 3 and 10 4 ; as the volume fraction increases from 0 to 0:05; the inverse solutions approach to the exact values. 5.3. Case III: Nonuniform time dependent heat ‡ux. As the last case, heat ‡ux is considered as q(y; t) = sin( t) cos( y) (52) which is nonuniform time-dependent boundary condition. The Rayleigh number e¤ect is reported by giving the global heat ‡ux simulation and isotherm variation for both direct and inverse problems. In Figure (12) global heat ‡ux is presented for Case III by drawing the heat ‡ux at each iteration up to t = 1:0: From the …gure it can be expressed that the level of agreement between the exact and inverse solutions is better than the solutions in the Case I and Case II.
The corresponding isotherms to the same analysis are presented in Figure (13). From the …gures it can be concluded that similar to the steady heat ‡ux case, there is no noticeable di¤erence on the isotherms as the Rayleigh number increases. As in Case I, this situation can be explained by the weak convection in the system.  From the Figure (14) it can be understood that up to t = 0:2; the noisy data does not a¤ect the heat ‡ux. But after t = 0:2; there is a deviation occurs for all values of and this deviation does not grow further and exhibits unchanging behavior up to t = 1: Particularly for = 0:05, the deviation in the graphs is much more pronounced than the lower values of . Similar result can also be inferred from Figure (15). In the case of no noise ( = 0), the numerical results agree with the exact solution. But when the noisy data is used in the solution of inverse problem, a de ‡ection is occurred in the numerical results, particularly = 0:05, as t approaches to the ultimate time 1: In Figure (  can be explained by the fact that T m ; which is obtained from the direct problem and considered as a sensor for the inverse problem, is located at the left wall which is the farthest away from right wall. Thus, it can be concluded that the location of the sensor is an important factor for the solution of the inverse problem.

Conclusion
In this study, the solution of two-dimensional unsteady inverse natural convection problem in a square enclosure …lled with Cu-water nano ‡uid has been obtained by using DRBEM. The method is used for the solutions of both direct and inverse problems. Computations are carried for three types of adiabatic boundary conditions, namely, steady heat ‡ux (q(y) = cos( y)), time dependent uniform heat ‡ux (q(t) = sin( t=t f )) and nonuniform time dependent heat ‡ux (q(y; t) = sin( t) cos( y)). DRBEM is an extension of BEM and is used to transform the domain integrals caused by non homogeneous terms of the partial di¤erential equations occuring in the BEM procedure into the boundary integrals. So DRBEM does not need to discretize the domain integrals and this is the main advantage of the method over the domain discretization methods. In the DRBEM Figure 15. E¤ect of the percentage of the noise on the global heat ‡ux for Ra = 10 3 ; = 0:05 and = 10 6 for Case II.
solution procedure, all the terms except the Laplace or modi…ed Helmholtz operators depending on the fundamental solution to be used are considered as nonhomogeneous terms. The nonhomogeneous terms are approximated using the radial basis functions.
The results are given for a range of Rayleigh number (0 Ra 10 4 ), particle volume fraction (0 0:05) and percentage of the noise (0 0:05). The e¤ects of the parameters on the inverse heat ‡ux and isotherm patterns are presented graphically. The results of this study can be summarized as follows: Figure 16. E¤ect of the percentage of the noise on theTemperature T for Ra = 10 3 ; = 0:05 and = 10 6 for Case II.
The analyses done for di¤erent values of Rayleigh number to show that in all three cases the Rayleigh number has a crucial e¤ect in obtaining heat ‡ux from the inverse problem. Reasonable accurate and stable numerical results are obtained at various values of Ra using both noiseless and noisy temperature data. Especially, it is observed that when the heat transfer mechanism is dominated by weak conduction or convection, the heat ‡ux results obtained from the inverse problem are much better. On the other hand, when the Rayleigh number takes the higher values, it not possible to obtain stable solutions from the inverse problem. Accuracy and stability of the inverse solutions also depend on type of the boundary conditions. When the graphics of global heat ‡ux in all three cases are compared, it can be concluded that the closest solutions to the exact solution is obtained by using the non uniform time-dependent heat ‡ux. Also, when steady and non-uniform time-dependent heat ‡uxes are used, the heat transfer in the system is governed by weak conduction, while it is governed by convection when time-dependent heat ‡ux is used for the same Rayleigh numbers.
The sensor position becomes key factor to obtain a stable results when the noisy temperature data is used in the solution of inverse problem. As used in this study, if the sensor is located at the left wall which is the farthest away from right wall, the heat ‡ux at the right wall is obtained with larger deformation. In addition, the percentage of the noise is another important factor for the inverse solutions. When the increased value of percentage of the noise is used, the deformation becomes more pronounced.
These results mentioned above are good agreement with the results for the inverse natural convection problem and inverse heat conduction problem in [19] and [17], respectively.
The temperature values at the left wall of the cavity obtained from the direct problem are considered as a sensor for the inverse problem. These values are used as overprescribed boundary condition for the solution of inverse problem. In this case, the accuraccy of the solutions of the direct problem at the boundary become quite important. DRBEM enables us to get unknowns and their normal derivatives at the boundary without discretizing the whole domain. Thus, the stable and reasonable accurate results can be obtained at a small computational expense. DRBEM is the most appropriate numerical solution technique for the solution of unsteady inverse natural convection of nano ‡uid problem. It can be concluded that the DRBEM solutions are well regularized for the inverse problem.

Declaration of Competing Interest
The author has no competing interests to declare.