The study of thermodynamic properties of complex fluid involving electrical solvent is of interest of many biological and physiological applications [21, 16, 7]. For example, electrokinetic microfluidic devices have been widely used in biomedical and biotechnological applications and chemical synthesis, which can be used for pumping and controlling the liquid flow in microfluidic systems or separating constituents suspended in a liquid [16, 12]. While low thermal conductivity devices are widely used, temperature effect should be includeded as it affects the electrical conductivity and can produce substantial heat driven flow. In addition, the interaction between heat flux and electric field could result in Joule heating effect as a byproduct of presense of electric field [25, 13]. This effect is generated by the ohmic resistance of the electrolyte subjected to electric current and, clearly, the traditional isothermal models aer not adequate in such situation. More recent works account for thermal effects and show that these interactions can be modeled in a thermodynamic consistent way (see ).
Another important pool of applications stems from biology, namely, modeling ion channels. Such channels are formed by charged walls and play a fundamental role in controlling and regulating the nervous system. In addition to the three main types of ion channels: voltage gated, ligand gated and chemically gated ion channels, temperature controlled ion channels are also been found and studied in the literatude (see [4, 24]). For example, six members of the mammalian Transient receptor potential ion channels (TRP) respond to varied temperature thresholds. TRPV1, TRPV2, TRPV3, and TRPV4 are heat activated, whereas TRPM8 and TRPA1 are activated by cold .
Despite there is vast amount of works devoted on related experimental studies[13, 16], mathematical modeling[27, 22] and simulation, , thermodynamically consistent models are much less studied and developed.
To begin the discussion of the non-isothermal model we first consider the isothermal case and the Poisson-Nernst-Planck (PNP) system. This model is widely used to describe the transport of charged particles [23, 1, 11, 28]with low concentration and couples the well-known drift-diffusion equations for the concentration with the Poisson equation for the electrostatic potential with appropriate initial and boundary conditions :
|(Drift Diffusion equation)|
Here 0 is the diffusivity, is the valence number, and is the mobility of the -th ion species. The PNP system models the interaction of N ionic species through an electrostatic field, usually . is the electrostatic potential, stands for the charge density of the ith species.
Both ion channels and electroosmosis can be modeled by taking the background fluid into consideration and coupling the PNP system with the Navier-Stokes equation. The energetic variational formulation which deals with the isotropic case for the resulting system from a mechanical prospective is well elaborated in[31, 14]. The energy law for this system, as stated in  is:
Here, denotes the part of boundary where the Robin boundary condition: is imposed. The total energy consists of the sum of the kinetic energy of the background fluid, the electro potential energy and the free energy of each of the species. The dissipation of the total free energy is driven by the viscosity and diffusion. This coupled system can be derived through energetic variational approach (EVA, introduced in ) and is as follows:
|(Drift diffusion equation)|
As is evident from the equations above, this model does not take temperature effect into consideration.
A number of numerical solution approaches for the PNP system and relevant mathematical models models have been considered in the literature. For example, A. Flavell et al. (2014)  applied a conservative finite difference scheme which achieves second-order accuracy in both space and time and also conserves total concentration for each ion species. C. Liu et al. (2015)  used an energetically stable finite element method for the PNP-NS(Poisson-Nernst-Planck-Navier-Stokes) system. D. Xie et al. (2016)  used a non-local finite element method for the PB (Poisson-Boltzmann) equation to tackle the problem of solution singularity caused by point charge term.
In order to capture the temperature change ,  considered a thermodynamic approach for incompressible Newtonian fluid coupled with temperature model. The resulting model is based on the first and second law of thermodynamics and in accordance with the Fourier’s Law assumes that the internal energy is proportional to the temperature and the heat flux is proportional to the temperature gradient. The equation for the temperature is derived through energy conservation and the whole system satisfies the inequality constraint for entropy production and this is in agreement with the second law of thermodynamics.
Some works in biology ,, take both temperature and electrokinetics into consideration by a simple coupling which adds electrokinetic force into the Navier-Stokes equation and use PNP or Poisson-Boltzmann to model the electric field. The temperature effects are incorporated by including the Joule heat force term into the heat equation. One issue with this model is that it is unclear what the energy associated with this system, and, moreover, it is not evidennt whether the second thermodynamics law is violated by the model.
In our study, we are using a more consistent model to capture the interactions mentioned above. According to first law of thermodynamics, heat and work can convert to each other, so it is natural to consider this energy interchange in order to derive the equation for the temperature. By the second law of thermodynamics, the conversion between heat transfer and work follows the rule that keeps the entropy increasing. We follow the ideas in  and add this as an inequality constraint for the system. As it turns out, maintaining such inequality constraints is also crucial for the stability of the discretized model. This motivated the choice of discretization and we have applied the Edge Average Finite Element (EAFE, see ) discretization, which satisfies the discrete maximum principle for the temperature. The discrete maximum principle turns out to be the key to numerical stability and one of the key results that we present is the energy estimate satisfied by the numerical model.
The paper is organized as follows. We introduce the Heat-PNP equations in section 2 and the corresponding energy law in section 3. In section 4 we propose our discretization and prove an energy estimate for the discretized Heat-PNP system. In section 4 some numerical experiments are provided to validate the stability of our numerical scheme. We also provide numerical tests which show consistency with qualitative phenomena observed experimentally .
2 The Heat-PNP equations and their discretization
As it is well known  the PNP system is as follows:
We treat the evolution equations for ions’ concentration as a system of continuity equations. After introducing the flux variable for each ion species, this sytem can be written as:
We use the same fashion to rewrite the PNP system with temperature throughout this paper.
We consider the the temperature PNP system in without background fluid. According to , this system has the form:
Here, is the electrostatic potential , stands for the charge density of the ith species, T is the temperature function and is the macroscopic velocity of the ith species. All these variables are space time functions defined on . 0 , , and are correspondingly the diffusivity constant, the valence number and the constant mobility of the ith ion species.
2.1 Initial and Boundary Conditions
For the system (2), the initial conditions at are
The boundary conditions are of different types: We have homogeneous no-flux boundary conditions for each of the concentrations of ion species,
For the Poisson equation, we partition the boundary into disjoint parts:
For the temperature equation we have Dirichlet boundary conditions which, as we show later help us in deriving the discrete energy law.
3 Energy of the Heat-PNP system
According to , the Heat-PNP system satisfies simultaneously first and second energy laws of thermodynamics, which are:
These laws involve the internal energy functional , the kinetic energy functional of the system, and the entropy functional of the whole system. In addition, is the applied electric field, which, for simplicity, we assume to be independent of time in our system here.
More precisely, if we substitute the predefined entropy functional and entropy production functional for Heat PNP system we obtain:
Next, we use the divergence theorem to derive the second law for our system in terms of the unknon variables:
In fact, to derive the second energy law of thermodynamics from the equation set is straightforward. We multiply by both sides of the heat equation (2) and integrate over the domain. For the left hand side we then we have
Next, for the right hand side we have:
Notice that in the derivation we used the continuity equations. Also, since this energy law is for closed systems, zero flux boundary condition is used when integrating by parts and finally we used that the heat source vanishes, namely, we have .
Clearly, the energy law (3.3) makes sense when ion concentrations and temperature are positive and the solution satisfies the following regularity assumptions:
3.1 Log-density formulation and its energy
We introduce a change of variables (also known as a transformation for the ion concentrations) which is as follows
We note that such change of variables requires that the concentrations are positive and we have
The embedding above holds for space dimensions . The Heat-PNP system then is written in a variational (weak) form: Find , with and such that:
Here we denote by the macroscopic velocity of ith species defined as , and for the energy law (3.3) we have
Furthermore, if we denote , , and for simplicity we ignore the constant , and under the appropriate zero flux boundary condition we arrive at the following simplified form of (3.3):
3.2 The discrete formulation
We consider a mesh of simplices covering our computational domain (triangles in2D or tetrahedra in 3D). As approximating space we take the space for piecewise linear, with respect to , continuous polynomials ,
The finite element solution to the Heat-PNP system then is defined using these finite element spaces. As a time marching scheme we choose the backward Euler scheme which is implicit and therefore stable. The discrete variational form then is: Find and
The approximation to the initial conditions uses the standard interpolation operatorand is as follows:
3.3 A discrete energy estimate
Next result shows a discrete energy estimate which holds in case when the mesh is quasi-uniform mesh and the stiffness matrix on each Picard iteration is an -matrix. In this case, we can show that the corresponding discrete solution satisfies maximum principle and with this condition satisfied, we have discrete energy estimate and mass conservation.
Suppose and satisfy equations (3.7) for and the assumptions for the discrete maximum principle of temperature equation are satisfied. Then the mass is conserved for each ion species
Moreover, the discrete analogue of the energy estimate (second law of thermodynamics) holds
where and h is the mesh size.
For simplicity, let , , , , and , and the corresponding discrete energy estimate is done as follows. First, we choose in the first equation of (3.7) and this gives us:
For the energy estimate, we test the last equation in (3.7) with . The latter is a valid test function because and is the interpolation operator mapping to to the piecewise linear, continuous .This gives us:
We analyze each term in the above equation sum and we begin by rewriting I on the left hand side. We have
Next, we rewrite also II using integration by parts and the zero flux boundary conditions:
In the last step we used the discretized continuity equation to estimate term III on the left hand side as follows:
Here we have used the assumption that , for any and that these constants are uniformly bounded by the constant . As shown in [5, Theorem 3.1.6] we have the following interpolation estimate
term4 on the right hand side:
Combining all the terms above and by a discrete Grnwall argument of telescopic sum from time step t=0 to m and, as a result, we obtain the following energy estimate:
only keep the term in I to the right side, we have: