CC BY
12
УДК 519.63
Energy Method for Mathematical Modeling of Heat Transfer in 2-D Flow
Valery V. Denisenko*
Institute of Computational Modelling RAS SB Academgorodok, Krasnoyarsk, 660036 Institute of Mathematics and Computer Science Siberian Federal University Svobodny, 79, Krasnoyarsk, 660041
Russia
Received 10.07.2014, received in revised form 10.08.2014, accepted 28.09.2014 A multigrid finite element method for 2-D convection-heat transfer problem is proposed. The method is based on minimization of the energy functional. Effectiveness of the method is demonstrated by calculation of temperature distribution in a microchip that separates living cells from fluid stream by electrophoresis.
Keywords: mathematical modeling, energy method, elliptical equation, nonsymmetric operator, elec-trophoresis
The problem of heat transfer in a moving medium and similar convection-diffusion problem arise in various fields of research, for example, in mathematical modeling of the transport of pollutants and energy in the atmosphere or in mathematical modeling of devices created for biological research.
Heat transfer is described by second order elliptic equation with proportional to velocity first order term. Then operators of boundary value problems are asymmetric operators. Because of this, the principle of minimum of a quadratic energy functional is not formulated. For problems with symmetric operators this principle allows one to use the most effective methods for the approximate and numerical solutions [10]. For problems with asymmetric operator the principles of nonequilibrium thermodynamics of minimal entropy production are not formulated, etc. [11]. Some problems with asymmetric operators are reformulated as problems with symmetric positive definite operators, and for them the energy principles are proved in [2,4]. Three-dimensional problems are discussed in [3].
The aim of this work is to create a multigrid finite element method for solving two-dimensional heat transfer problem in a moving medium and to demonstrate effectiveness of the method.
As an example temperature distribution in the microchip, that separates living cells from a fluid stream by means of electrophoresis is calculated.
1. The original problem
According to [9] stationary heat transfer equation in inviscid incompressible fluid has the form
pCu ■ grad T - div (A grad T)= Q, (1)
*[email protected] © Siberian Federal University. All rights reserved
where T is the temperature, p is the density of the fluid, C is the heat capacity per unit mass at constant pressure, A is the coefficient of thermal conductivity, u = (u, v) is the velocity vector of the fluid and Q is given density of the energy source. The same equation is included in the system of equations describing the heat transfer in a viscous fluid. Then Q includes the production of thermal energy due to friction.
We solve the equation in the two-dimensional domain Q with boundary r. To use the embedding theorem by S.L. Sobolev, we assume that the domain Q is a connected union of a finite number of domains. Each domain is starlike with respect to some circle. It is assumed that boundary r is piecewise smooth. We use the arc length I as the coordinate on the boundary. Indices n, I denote the normal and tangential components of the vectors. The positive normal is directed outward and positive sense of rotation is counterclockwise.
We consider the Dirichlet problem when the temperature is specified on the whole boundary:
T lr = To(l). (2)
We assume that heat capacity C is constant along each stream line. This assumption is true, for example, for perfect gases
u • grad C = 0. (3)
Coefficient of thermal conductivity is always positive. We assume that it is uniformly bounded from above and below. It means that there are positive numbers £1 and £2 such that the following inequality
0 < £1 < A < £2 < ^ (4)
is true for all points of the domain.
2. Description of heat transfer in gyrotropic medium
In the stationary case, the mass conservation law has the form:
div(pu) = 0. (5)
In view of (3, 5) div(Cpu) = Cdiv(pu) + pu • gradC = 0. For any divergence-free vector function one can construct a stream function p(x, y) such that
pC u = -rot(^ez), (6)
where ez is the unit vector along the z axis. Here we calculate only x and y components of the rot operator.
We suppose that p(x, y) is a limited function. Thereby we exclude from consideration only flows with infinite parameters.
Then equation (1) can be written as
-d^(( P -A) gradt) = Q. (7)
Tensor coefficient is denoted by <r. It is symmetrical only if p = 0 when fluid is at rest. Symmetric part of <r is positive definite due to (4). The vector
j = -a grad T (8)
is the density of some current but it does not coincide with the true density of the heat flux which divergence is given by equation (1). Density of the heat flux is j — rot(^Tez). The difference has zero divergence and that is why the law of energy conservation (1) can be written as conservation law for some auxiliary transfer process occurring from the law (8).
Let us consider tensor a. As is evident from (7) the tensor is invariant with respect to rotation about the z —axis. It defines rotation around this axis and isotropic extension. Therefore, (8) is the law of heat transfer in a gyrotropic medium.
By analogy with the problems of the electrical conductivity we denote E = —grad T. Since rot of grad is identically zero the vector E satisfies the equation:
rotz E = 0. (9)
The problem to determine vector E is equivalent to problem of finding T(x, y) if we also set the average value of T on Q or r. Therefore, the original problem (1), (2) can be replaced by the following problem:
divj = Q, rotzE = G, j = aE, E | r = g(l). (10)
Boundary condition is obtained by differentiating (2) along the boundary. As this takes place, the average value of T0(l) disappears. It is needed to determine T(x,y) after the solution of problem (10).
It is usually assumed that G = 0. The necessary condition for solvability of (10) is obtained from the equality of two expressions for circulation of vector E. The first expression is obtained by integrating the second equation over the domain and the second expression is obtained by integrating boundary condition over the boundary. Then we have
jj Gdxdy = j) gdl. (11)
We assume that the integrals of functions Q, G, g squared are finite. This condition can be replaced by the weaker condition. It is referred to as the divergence condition in [10]. It allows the existence of certain singularities. For example, such singularities occur near the edges of the electrodes as in example solved below.
Heat transfer problem (10) or problems with other boundary conditions can occur regardless of problem (1), (2). Specific for a direction of rotation around the z—axis can be associated with presence of magnetic field strength in z direction when the heat-carrying particles have a charge or if the medium rotates.
3. Energy principle
The operator of boundary value problem (10) is asymmetric as operator of the original problem (1), (2). The way of its symmetrization is presented and justified in [1]. Here we present only the main formulae and statements.
Let us introduce two new unknown functions F and P. The old unknowns j, E can be expressed in terms of F and P as follows
j = a E = — — aS<7*grad F + aSrot P, (12)
ao
where P is z-component of the vector (0,0, P); S is arbitrary symmetric and uniformly positive definite tensor with bounded coefficients; a0 is a constant which makes equations dimensionless and symbol * denotes the transposition.
The best estimations are obtained for
S-1 = (a + a* )/(2ao), ao = (13)
Here the constant £1 is constant from condition (4) and £3 is chosen so that the following inequality in the domain Q is satisfied
A + p2/A < £3 < to.
This condition is satisfied because we assume that p is bounded.
Let us consider the set of pairs of functions F, P, satisfying the following conditions:
F |r =0, ^Pdxdy = 0. (14)
n
We introduce the energy scalar product
P MP :
1 [f grad F \* -L âSâ * -<r S \ ( grad F ' , , , ,1C, rot P J -W* a0S К rot P' (15)
Energy space of the problem is the set of pairs of functions, satisfying (14) and having bounded energy norm.
Let us define the energy functional
W (F,P ) = 2
PMP
^(FQ/Oq + PG) dxdy + j> Pg(l) dl. (16)
It is proved that the energy functional has unique minimum with respect to functions of the energy space. The minimum value of the energy functional is attained at a pair of functions. Such pair is a generalized solution for the problem which is written here in a simpler form based on tensor <r given in (7) and relations (13):
/A2 + p \
div (--^ a grad F + a grad PJ = Q/oo,
divf A grad F - 00 grad Pj = G, (17)
в dF a0 dP A dn A dn
= g(i). (18)
г
The boundary condition (18) is equivalent to the condition (10). The main boundary conditions (14) originally define the set of admissible functions F, P. The existence of generalized solution of (17), (18) is proved. The problem of diffusion in gyrotropic medium (10) has the same generalized solution. If the functions F, P are smooth we have classical solution. Uniqueness of the classical solution for (10) is also proved.
Function T is the solution of problem (10) plus an arbitrary additive constant. The constant is determined by the average value (2) given on the boundary.
The proposed symmetrization of the problem can be regarded as right multiplication of the initial operator by adjoint operator. It is usual practice to do when the potential for unknown (usually vector) functions are introduced, although usually it is explaind as satisfaction of one of the equations of the original system. For example, we can observe this in the transformation of the Cauchy-Riemann equations to the Laplace equation.
i
4. Implementation of finite element method
Since the principle of minimum of the quadratic energy functional (16) for the elliptic boundary value problem (14), (17), (18) is justified it is easy to apply the finite element method. One need to triangulate the region, choose the approximating functions and write down the conditions of minimum energy in relation to the parameters of these functions. We use block-structured grids and piecewise linear approximating functions. Nodal values satisfy the system of linear algebraic equations. The system follows from laws (10), and hence from original equation (1), integrated over the half-integer grid cells. Linear relation j = aE (10) is satisfied with tensor function a averaged over each triangle of the mesh.
The matrix of this system is symmetric and positive definite as a consequence of similar properties of the boundary value problem. We solve the system of equations with the use of the multigrid method. From the point of view of the theory of finite element method the problem has all the properties which the Dirichlet problem for the Poisson equation has because the energy norm is equivalent to the norm of W2(1)(Q). The presented implementation of finite element method for the numerical solution of such problems is described and justified in [2].
5. The design of the microchip
To extract a small number of cancer cells from the blood for diagnosis one can use elec-trophoresis. It is movement of the dispersed phase in a fluid medium under the influence of the electric field. In flow type devices [5] cells approach the electrodes with the fluid flow and strong field keeps them in the vicinity of the electrodes without disturbing the fluid flow. The electric field in a conducting fluid generates an electric current and Joule dissipation Q = aE2, where a is the conductivity and E is the electric field strength. To be effective electrophoresis demands high field strength but it produces too much heat that can kill cells. Therefore, an essential part of mathematical modeling of such devices is the calculation of heat transfer. To solve this problem, we use the method described above.
The channel of the microchip under consideration has a rectangular cross section. We use Cartesian coordinates with the axis x along the channel and y—axis directed vertically. The fluid flows along the axis x in a flat layer 0 < y < y1. A sufficiently wide channel can be considered in the context of 2-D model so all parameters do not depend on the coordinate z. Fig. 1 shows cross section of the microchip in the plane x, y.
The bottom wall is made from copper, which is a good conductor of heat. It is kept at a constant temperature, which is due to the linearity of the problem can be equal to zero. The upper wall is a glass plate at y1 < y < y2. The outer surface of the glass plate has the same temperature due to air flow. The channel is long and fluid that enters the channel has the same temperature. This can be approximately described as supporting the same temperature for long enough distance x1. Since both channel walls have the same temperature and the source of energy is concentrated in a small area, the temperature of the fluid for long distances downstream tends to wall temperature. It can be approximately described as supporting the same temperature for sufficiently long distance x2. We choose the values of the parameters x1 and x2 in the test calculations so as not to disturb the solutions in the domain of interest but the computational domain should not be too large. It permits to improve the efficiency of calculations. Thus we obtain a Dirichlet condition on the whole boundary
TI =0. (19)
y = y2 = 0.3 mm
glass
y = yi =0.1 mm
dV = 20 V dl = 30 ^A/mill
copper
y=0
Fig. 1. The design of the microchip and electric field distribution: thick lines — equipotentials, thin lines — current lines
6. Calculation of the electric field
To determine the spatial distribution of the Joule dissipation which plays the role of the source of energy density in equation (1), it is necessary to solve the problem of the electrical conductivity and compute Q = aE2. The problem is formulated in this section.
In the frequency range of interest a vortex electric field can be ignored. This allows us to introduce the electric potential V such that the electric field E is represented as E = —grad V.
Potential V satisfies the electrical conductivity equation. In a homogeneous conductor the equation has the form of the Laplace equation
AV = 0.
(20)
The electric field and the current in this device are generated by a pair of electrodes. The first electrode takes the entire bottom wall ri. The second electrode is the thin strip deposited on glass and it takes the region r2: y = y1, |x| < x0. The electrodes have potential values equal to 0 and V0, respectively. Glass is a good insulator so the normal component of the current density is equal to zero at its surface. Because of the isotropy of the conductor it means that the potential derivative in the normal direction is negligibly small. Electric current flows from the second electrode to a small portion of the bottom wall which is comparable with the height of the channel. Asymptotic condition of potential decreasing with increasing distance from the second electrode assumes that it is vanishingly small at sufficiently large distance x3. In such a way we obtain a mixed boundary value problem for the Laplace equation (20):
V In
0,
V Ir
Vo.
dV dy
0,
y=yi,|x|>xo
dV
dx
0.
(21)
x=±x3
Let us note that the Laplace equation (20) is a special case of system (17) when p = 0. Therefore, boundary value problem (20), (21) can be solved by the same multigrid finite element method [2]. The result is shown in Fig. 1. The equipotential lines and current lines are presented. The following values of the input parameters are used: the channel height y1 =0.1 mm, the conductivity a = 0.00176 / m corresponds to sugar solution which provides the same osmotic pressure as blood [8], the voltage V0 = 220 V. It is enough to take x3 = 0.4 mm, so as not to
2
make a disturbance in the solution that it is seen in Fig. 1, 2. Total current I = 0.52 mA/mm and total Joule dissipation W = 115 mW/mm are obtained.
Fig. 2 shows the corresponding distribution of the density of the Joule dissipation Q in the fluid. We use the logarithmic scale. The values on adjacent lines differ by the factor 101/4 ~ 1.8. Exponential decay at |x| >> y1 is seen.
Fig. 2. The density of the Joule dissipation. Values at adjacent lines differ by the factor 101/4 ~ 1.8. Bold segment is the second electrode. The solid curves show the lines with values Q > 1 W/mm3 and dashed lines are smaller values
The boundary condition at the second electrode (21) gives the value of potential. The normal derivative of potential is given at the rest of the flat surface y = y1. In such cases, the solution of the Laplace equation (20) and hence the function Q have the singularity of the type 1/r, where r is the distance to the point where the type of boundary condition is changed. It can be shown that singularities in the temperature distribution do not arise since these functions satisfy the divergence condition [10]. It will be seen in the resulting temperature distribution.
7. Poiseuille flow
For the considered flow the Reynolds number is of order unity. Therefore, the flow is laminar. Because of the pressure drop in the direction of the axis x for sufficiently large length and width of the channel the Poiseuille flow is formed. The velocity in Poiseuille flow is directed along the axis x and has a parabolic profile:
u = (u(y), 0), u(y) = uo^ 1 21 - 1) )> £ =3PCuoy1^(y0 - y0 - ^ . (22)
We use the following parameters in our calculations: water density p = 1000 kg/m3, C = 4200 J/(kg-K), the thermal conductivity A = 0.7,0.58, 384 J/(m • s • K) for glass, water and copper, respectively, the thickness of the glass y2 - y1 = 0.2 mm, the velocity at the center of the layer u0 = 15 mm/s. The current function is obtained such that 0 < £ < 4.2 J/(m-s-K).
The exact solution of the one-dimensional heat transfer towards uniform flow is T = To exp (x/xb). An important parameter characterizing this solution is the thickness of the boundary layer xb = A/(u0pC) = 0.009 mm. It is about one-tenth of the height of the channel y1 =0.1 mm that is small but discernible in the scale of interest.
8. Heat transfer in the microchip
Numerical solution for boundary value problem (14), (17), (18) with density of the energy source Q presented in Fig. 2 and G = 0, g = 0 gives the temperature distribution in the structure
described above. Strictly speaking, it is necessary to include the dissipation of mechanical energy due to viscosity into the energy source but in this microchip the dissipation is negligibly small in comparison with the Joule dissipation.
Fig. 3 shows the computational domain and the coarsest mesh for which the results are presented below. We should note that we mean the main grids and the multigrid method coarsens the mesh until there is only one line of nodes within the domain.
y2 = 0.3 mm
-0.35 y = 0 x, mm 0.9
Fig. 3. Example of the computational grid, 16 * 64 cells. Bold segment is the second electrode
Fig. 4 and 5 show the obtained distribution of temperature when the fluid is at rest and when the fluid moves. The mesh has 8 times more cells in each direction in comparison with that shown in Fig. 3, that is 128 * 512 cells. Fig. 5 also demonstrates the error introduced by the transfer of the asymptotic boundary conditions to finite distances. It is enough to take x1 = —0.35 mm, x2 = 0.9 mm so as not to make a disturbance of the solution that would be evident in the following figures. With further expansion of the domain the shift of the line in the figures is less than their thickness. No difference is also seen from the results obtained with the use of 64 * 256 cells, i.e. the approximation error is small in this scale.
y2 = 0.3 mm
Fig. 4. Temperature distribution when the fluid is at rest
Fig. 6a shows the corresponding temperature distribution at the vertical line passing through the right edge of the second electrode. Sections are given for a set of grids. The first grid is twice finer than one shown in Fig. 2, i.e., 8*32 cells. The last grid is 8 times finer than one shown in Fig. 2, i.e., 64*256 cells. Finer grids provide solutions that are indistinguishable in the scale of the figure. We see the convergence of the solutions when grid size is decreased. Let us note that the function T is naturally calculated on the lines of half-integer grid in terms of the potentials F, P. In contrast to the integer cells half-integer lines are not in place when grid is refined as
y2 = 0.3 mm
Fig. 5. Temperature distribution. The bold lines is the solution in the domain -0.35 < x < 0.9 mm. Thin lines is the solution in the domain -0.25 < x < 0.7 mm
shown in Fig. 6a. We use linear interpolation between half-integer points. These graphs show the error of approximation. For smooth solutions it is quadratic with respect to the step size of the grid. In our case, the error decrease is slower due to the presence of singularities in the solution.
Fig. 6b shows the convergence of multigrid iterations. Iteration parameter in the method of successive over-relaxation(SOR) t = 1.75. Thin lines correspond to the results when the fluid is at rest. In this case, the optimal t = 1.4. As you can see, the rate of convergence of multigrid iterations is reduced by 2.5 times when the motion of the fluid is taken into account.
Fig. 6. a — Convergence with grid refinement. Cross section T(x, y) at x = 0.025 mm corresponds to the edge of the electrode. Thin lines correspond to the fluid at rest; b — The convergence of multigrid iterations in the norm || H ||= £|Hj|, divided by the norm of the right-hand side. I is the number of multigrid iterations. Thin lines correspond to the fluid at rest. Curves for meshes from 8 * 32 to 64 * 256 — upwards
A twofold increase of the fluid velocity provides a twofold decrease of the rate of convergence. Let us note that with this increase of the velocity the boundary layer thickness xb decreases to 0.005 mm and it becomes equal to the smallest grid size in a remote part of the downstream region. In such cases, the accuracy of the solution also decreases. At high speeds, numerical methods based on the separation of the flow direction and boundary layers should be used [12,13]. Our method is designed to solve problems with a relatively slow movements of the complex structure
when the implementation of the scheme with the separation of the flow direction and boundary layers is difficult. In particular, the problem solved here contains the regions with and without flow.
Let us note that the solution reaches the accuracy shown in Fig. 6a after 4 multigrid iterations. The accuracy shown in Fig. 5 in the downstream region is achieved after 8 iterations. This means that accuracy becomes equal to the approximation error and further iterations do not make sense. Each multigrid iteration is equivalent in respect of number of arithmetic operations to 10 iterations of SOR method.
9. The condition numbers of the matrices
An important characteristic of a system of linear algebraic equations is its condition number. For symmetric positive definite matrices, which are under consideration here, the condition number is defined as the ratio of the maximum Amax and minimum Amin eigenvalues [7]. This number, in particular, is the ratio of the relative error of the solution to the relative error of the right-hand side of the system of equations. Since the condition number for our problem is about one million, we carry all calculations with double precision. The accuracy of these 8-byte numbers ~ 10-15.
As a test, first we consider the special case where the fluid is at rest and thermal conductivity A =1 J/(m • s • K) in the entire domain. Then boundary value problem (14), (17), (18) turns in two independent Dirichlet and Neumann problems for the Poisson equation for the functions F and P, respectively. We use a uniform grid and slightly reduce the domain in the direction of x to x2 — x1 = 1.2 mm. Then the grid sizes h in both directions are equal. As already mentioned, the equations of our finite element method mean the laws (10) integrated over the half-integer cells. We multiply equations corresponding to the corner grid points by 4 and equations corresponding to the rest of the boundary nodes by 2. Then all the equations differ from the equations used in the theory of difference schemes by only common factor h2. With these caveats, in this particular case, the equations of our finite elements are transformed into the equations of classical five-point difference scheme for the Poisson equation. Eigenfunctions and the corresponding eigenvalues of the difference schemes of the Dirichlet and Neumann problems are known. For a united matrix, combining the equations for the nodal values of F and P, we obtain Amax = 8 and Amin = 2(1 — cos (n/(x2 — x1))) ~ (n/Nx)2, where Nx is the number of steps in the direction of x in which the region is more extended. Test calculations give these values with high accuracy. We use simple methods to find the maximum and minimum eigenvalues: power method and the method of inverse iteration. There are much more effective methods but these are the simplest ones. The condition number increases as about 1/h2 when h ^ 0 and for the finest grids 128 * 512 cells we have Amax/Amin ~ 2.1 • 105.
If we use not uniform grids similar to that shown in Fig. 3 the condition number is increased by the factor of less than 1.5 for the same number of nodes.
For the problem with the moving fluid for the grid 128 * 512 cells we obtain Amax ~ 232, that is, the condition number is increased by 30 times in comparison with the problem when fluid is at rest. This scale roughly corresponds to the maximum eigenvalue of the matrix of coefficients in system (17)
1/(A2 + ^2)/ac — £
A —£ ^Q
that reaches 30 for the predetermined range of parameters values and aQ = 1.
For this problem Лт1п ~ 6.7• 10-5 and the condition number ~ 3.5 • 106, that is, it is increased by 17 times in comparison with the problem when fluid is at rest.
Conditionality can be improved almost threefold by reducing the |в| twofold since the function is defined up to an additive constant. A half of the maximum value of в can be deducted from the constructed function в. Additional small improvement can be achieved by choosing the constants ст0 instead of ст0 = 1 used here. This optimization in в and ст0 has no effect on the rate of convergence of the iterations shown in Fig. 7. This is due to the fact that we solve the equations for the pair of F, P at each node during iterations of SOR method at each grid in our multigrid method. This is equivalent to matrix (23) inversion when the coefficients are constant. This is close to the matrix inversion for sufficiently fine grids when coefficients differ little in neighboring grid cells.
Let us note that for the solution of original equation (1) it is sometimes recommended to use the least squares method, This is actually equivalent to multiplying the original equation by the adjoint operator from the left. In the particular case of a stationary homogeneous fluid the Poisson equation in this approach is replaced by the biharmonic equation. The condition number of matrices arising in the approximation of the biharmonic equation is of order 1/h4 [6]. This is much worse than the 1/h2 in our approach and in the general case the estimate of 1/h4 cannot be improved.
Thus, the proposed method for the problem of heat transfer in a two-dimensional fluid flow without thin boundary layers is almost as effective as the classical multigrid methods for the Poisson equation.
The obtained results of heat transfer in the microchip can be used in designing such devices. This work was supported by the Russian Science Foundation under project 14-15-00805.
References
[1] V.V.Denisenko, A boundary value problem for an elliptic equation in two variables with asymmetric tensor coefficients, Sib. Math. J., 35(1994), no. 3, 495-505.
[2] V.V.Denisenko, Energy Methods for Elliptic Equations with Asymmetric Coefficients, Publ. house of the Russian Academy of Sciences, Siberian Branch, Novosibirsk, 1995 (in Russian).
[3] V.V.Denisenko, Energy method for three-dimensional problems of transfer in moving media, Russian Journal of Numerical Analysis and Mathematical Modelling, 14(1999), no. 1, 37-58.
[4] V.V.Denisenko, The energy method for convection-diffusion problems, Zh. Prikladnoi Matematiki i Tehnicheskoi Fiziki, 38(1997), no. 2, 197-203 (in Russian).
[5] J.Gao, X.-F.Yin, Z.-L.Fang, Integration of single cell injection, cell lysis, separation and detection of intracellular constituents on a microfluidic chip, Miniaturisation for chemistry, biology and bioengineering, 4(2004), 47-52.
[6] L.V.Gileva, V.V.Shaidurov, Two multigrid iterative algorithms for a discrete analogue of the biharmonic equation, Sib. Zh. Chislennoi Matematiki, 7(2004), no. 3, 213-228 (in Russian).
[7] S.K.Godunov, V.S.Ryabenky, Finite Difference Schemes: Introduction to Theory, Nauka, Moscow, 1973 (in Russian).
[8] C.-T.Huang, T.G.Amstislavskaya, G.-H.Chen, H.-H.Chang, Y.-H.Chen and C.-P.Jen., Selectively concentrating cervical carcinoma cells from red blood cells utilizing dielectrophoresis with circular ITO electrodes in stepping electric fields, Journal of medical and biological engineering, 33(2012), no. 1, 51-58.
[9] L.D.Landau, E.M.Liphshits, Theoretical Physics, V. 6, Hydrodynamics, Moscow, Nauka, 1986 (in Russian).
[10] S.G.Mikhlin, Variational Methods in Mathematical Physics, New York, Pergamon Press, 1964.
[11] Yu.B.Rumer, M.Sh.Ryvkin, Thermodynamics, Statistical Physics and Kinetics, Mir, Moscow, 1980.
[12] V.V.Shaidurov, G.I.Shchepanovskaya, M.V.Yakubovich, Numerical simulation of viscous heat conducting gas in a chanal, Vychislitelnye Tehnologii, 18(2013), no. 4, 77-90 (in Russian).
[13] E.O'Riordan, M.L.Pickett, G.I.Shishkin, Parameter-uniform finite difference schemes for singularly perturbed parabolic diffusion-convection-reaction problems, Mathematics of Computation, 75(2006), 1135-1154.
Применение энергетического метода при математическом моделировании переноса тепла в двумерном течении
Валерий В. Денисенко
В статье предложен многосеточный вариационно-разностный метод решения двумерной задачи конвекции,—теплопроводности, основанный на минимизации функционала энергии. Эффективность метода продемонстрирована при расчете распределения температуры в микрочипе, выделяющем живые клетки из потока жидкости с помощью электрофореза.
Ключевые слова: математическое моделирование, энергетический метод, эллиптическое уравнение, несимметричный оператор, электрофорез.
