Thermal Accelerometer Simulation by the R‑Functions Method

As well as many modern devices, thermal accelerometers (TAs) need a sophisticated mathematical simulation to find the ways for their performance optimization. In the paper, a novel approach for solving computational fluid dynamics (CFD) problems in the TA’s cavity is proposed (MQ-RFM), which is based on the combined use of Rvachev’s R-functions method (RFM) and the Galerkin technique with multiquadric (MQ) radial basis functions (RBFs). The semi-analytical RFM takes an intermediate position between traditional analytical approaches and numerical methods, such as the finite-element method (FEM), belonging to the family of the so-called meshless techniques which became popular in the last decades in solving various CFD problems in complex-shaped cavities. Mathematical simulation of TA by using the MQ-RFM was carried out with the purpose to simulate the temperature response of the device and to study and improve its performance. The results of numerical experiments were compared with well-known analytical and numerical benchmark solutions for the circular annulus geometry and it demonstrated the effectiveness of the MQ-RFM for solving the convective heat-transfer problem in the TA’s cavity. The use of solution structures allows one to take a relatively small number of expansion terms to achieve an appropriate accuracy of the approximate solution satisfying at the same time the given boundary conditions exactly. The application of the MQ-RFM gives the possibility to obtain semi-analytical solutions to the diffusion-convection problems and to identify the main thermal characteristics of the TA, that allows one to improve the device performance.


Introduction
Many modern devices, in particular, sensors based on microelectromechanical systems (MEMS), need a sophisticated mathematical simulation to find the ways of their optimization. Among them are thermal accelerometers (TAs) [1] which principle of operation is based on the effect of fluid or gas convection in closed cavities under the influence of external forces of inertia. It is very important to estimate the scale factor and bias stabilities of TAs under an external thermal slope, and to test different types of cavity geometry (cylindrical, rectangular, hexagonal, etc.) to achieve the best performance of the device.
To solve the aforementioned problems, it is necessary to simulate diffusive and convective heat and mass transfer in arbitrarily shaped enclosures, which in turn is a problem of great importance due to a large number of its practical applications. Since analytical solutions exist for a very narrow class of simple-shaped domains only, a majority of modern computational fluid dynamics (CFD) techniques

Thermal Accelerometer and Its Principle of Operation
The principle of operation of a TA is based on the effect of fluid or gas convection in a closed cavity under the influence of external forces of inertia. The device ( Figure 1) includes a heating element H that creates around itself a symmetrical thermal field. The thermal sensors are located on opposite sides of this element. In the absence of an external acceleration along the sensitive axis X, the system is balanced and the heat detectors generate the same signal ( Figure 2a). In presence of an external acceleration, the temperature profile changes, which results in the temperature difference dT between the sensors S1 and S2, depending on the amplitude of the acceleration (Figure 2b). Nonzero temperature difference between the heat sensors converts the input impedance to the output electrical signal. This type of accelerometer has several important advantages over conventional sensors for acceleration based on microelectromechanical systems (MEMS). In particular, due to the absence of moving mechanical parts, convective accelerometers have high reliability, low cost of production, as well as the ability to withstand and measure high loads caused, for example, by shock action.
To find the optimal parameters of the device: the type of gas and its pressure, the size and geometry of the cavity, the material and the size of the heating element and heat sensors, etc., it is necessary to solve the boundary value problem described by a system of the Navier-Stokes differential equations. This problem can be solved analytically only for a limited class of simple regions (cylinder, sphere), and in the general case, a numerical simulation is required. In the absence of an external acceleration along the sensitive axis X, the system is balanced and the heat detectors generate the same signal ( Figure 2a). In presence of an external acceleration, the temperature profile changes, which results in the temperature difference dT between the sensors S1 and S2, depending on the amplitude of the acceleration (Figure 2b). Nonzero temperature difference between the heat sensors converts the input impedance to the output electrical signal.

Statement of the Problem and Governing Equations
Appl. Sci. 2020, 10, x FOR PEER REVIEW  3 of 16 and identify main thermal characteristics of the TA of more complicated geometry [37], that allows one to improve the device performance.

Thermal Accelerometer and Its Principle of Operation
The principle of operation of a TA is based on the effect of fluid or gas convection in a closed cavity under the influence of external forces of inertia. The device ( Figure 1) includes a heating element H that creates around itself a symmetrical thermal field. The thermal sensors are located on opposite sides of this element. In the absence of an external acceleration along the sensitive axis X, the system is balanced and the heat detectors generate the same signal ( Figure 2a). In presence of an external acceleration, the temperature profile changes, which results in the temperature difference dT between the sensors S1 and S2, depending on the amplitude of the acceleration (Figure 2b). Nonzero temperature difference between the heat sensors converts the input impedance to the output electrical signal. This type of accelerometer has several important advantages over conventional sensors for acceleration based on microelectromechanical systems (MEMS). In particular, due to the absence of moving mechanical parts, convective accelerometers have high reliability, low cost of production, as well as the ability to withstand and measure high loads caused, for example, by shock action.
To find the optimal parameters of the device: the type of gas and its pressure, the size and geometry of the cavity, the material and the size of the heating element and heat sensors, etc., it is necessary to solve the boundary value problem described by a system of the Navier-Stokes differential equations. This problem can be solved analytically only for a limited class of simple regions (cylinder, sphere), and in the general case, a numerical simulation is required. This type of accelerometer has several important advantages over conventional sensors for acceleration based on microelectromechanical systems (MEMS). In particular, due to the absence of moving mechanical parts, convective accelerometers have high reliability, low cost of production, as well as the ability to withstand and measure high loads caused, for example, by shock action.

Statement of the Problem and Governing Equations
To find the optimal parameters of the device: the type of gas and its pressure, the size and geometry of the cavity, the material and the size of the heating element and heat sensors, etc., it is necessary to solve the boundary value problem described by a system of the Navier-Stokes differential equations. This problem can be solved analytically only for a limited class of simple regions (cylinder, sphere), and in the general case, a numerical simulation is required.
Usually, the temperature distribution inside the 3D cavity of the device is almost uniform in the direction along the heater and sensors. Thus we need to solve the 2D heat-transfer problem in the middle cross-section of the cavity. One of the most important models of the TA is the circular annulus model (Figure 3a) [34,36] because for such a domain the analytical solution to the heat diffusion problem is known as well as the asymptotic solution of the diffusion-convection problem [24,25].
We consider the Boussinesq approximation with the assumption of the steady state incompressible flow inside an arbitrary 2D closed cavity  The dimensionless governing equations for the two-dimensional convective heat and mass transfer inside Ω have the following form: ; , , where ( )/( )   where n is the outward normal to boundary ∂Ω . We consider the Boussinesq approximation with the assumption of the steady state incompressible flow inside an arbitrary 2D closed cavity Ω ⊂ R 2 (Figure 3b), when all constants, except density, do not depend on temperature. Without loss of generality, we take the domain with boundary composed of two parts: ∂Ω C , ∂Ω H , with given constant temperatures T C and T H on each of them respectively. Let external acceleration Γ be applied towards the positive direction of sensitivity axis OX.
The dimensionless governing equations for the two-dimensional convective heat and mass transfer inside Ω have the following form: where θ = (T − T C )/(T H − T C ) is the dimensionless temperature; T is the temperature; ψ is the stream function; ζ is the vorticity; X = x/L, Y = y/L are dimensionless coordinates; L is the specific dimension; U = uL/ν; V = vL/ν are dimensionless components of velocity; u, v are horizontal and vertical components of velocity; ν is the kinematic viscosity; Gr = Γβ(T H − T C )L 3 /ν 2 is the Grashof number; Γ is the acceleration; β is the temperature volume expansion coefficient; Pr = c p µ/λ is the Prandtl number; µ is the dynamic viscosity; λ is the heat conductivity; c p is the specific heat at constant pressure.
Therefore, on ∂Ω the following dimensionless boundary conditions are given: where n is the outward normal to boundary ∂Ω. Boundary conditions for the vorticity function ζ are not defined explicitly and are usually approximated by its Taylor series expansion in the neighborhood of the boundary.

Generalized Numerical Procedure
First of all, before applying common variational techniques considered in [38,39], we should pass from inhomogeneous boundary conditions (2) to homogeneous ones by introducing the new function: where Φ is a continuous function which is equal to unity on the part ∂Ω H and vanishes on ∂Ω C . Function Φ can be constructed in the way similar to that for the traditional Lagrange interpolation polynomial but here instead of interpolation points we take corresponding segments of the boundary. This technique also is called the transfinite interpolation [40].
The iterative process (6) stops when some convergence conditions are fulfilled, for example, In a number of variational and projection techniques [6,38,39], approximations to θ (k) at each iterative step are found in the form of truncated generalized Fourier series with respect to functions of some basis f n N n=0 : where N ) are undefined coefficients. In addition to their differentiability, all functions f n must vanish on the boundary, i.e., From Equation (9) it follows that approximation (8) at any iterative step will also satisfy the homogeneous boundary condition (5) for θ exactly.
Leonid Kantorovich [38] proposed the following technique for constructing a functional basis satisfying Equation (9). First we should take an arbitrary basis system of functions {χ n } N n=0 . There may be algebraic or trigonometric polynomials, splines, etc. Then a function ω(X, Y) is constructed such that ω > 0 inside Ω; ω < 0 outside Ω ∪ ∂Ω; ω = 0 and |∇ω| 2 0 on ∂Ω. For simple domains, for example circular, these functions are trivial, but the common approach to obtaining analogous expressions for arbitrary domains was developed by V. Rvachev [10]. This approach (the R-functions method-RFM) will be considered in the next section.
Then we take the new system of functions which, due to the properties of the function ω, satisfy the conditions (9). Taking into account Equation (9), Equation (8) can be written as: To strictly obey conditions (5) for the stream function, we must take another basis ω 2 χ n N n=0 with undefined coefficients e (k) = (e (k) Then the velocity components will be expressed as Substituting Equation (12) into the left-hand side of the third equation of system (4), we obtain a vorticity function expansion: The latter expression may be considered as a representation of the vorticity function with respect to the set of basic functions ∇ 2 (ω f n ) with a set of unknown coefficients Thus, it is necessary to solve recurrently only two equations of System (6) to find coefficients c n . It should also be noted here that Equation (14) provides us an expression for the vorticity functions without the necessity to approximate additionally its boundary values.
At each iteration, the realization of various numerical techniques requires different procedures of reducing an original problem to a corresponding system of linear algebraic equations (SLAE) with respect to unknown vectors c (k) and d (k) . Here we propose to use an approach on the base of the Petrov-Galerkin procedure [6].
At the first step, we find the components c  (11) into the first equation of (6), we obtain the residual: Pr Appl. Sci. 2020, 10, 8373 7 of 17 Here the right-hand part is One must choose the set of coefficients c (k) n minimizing residual (16). In the Galerkin scheme, the orthogonality of δ θ to all functions of the system f n where dσ ≡ dXdY. Substituting Equation (16) in Equation (17), we get the SLAE with respect to elements of the vector c (k) : where components of matrices A, A are determined respectively as: Pr Pr Some optimization of Equation (19) can be done to decrease the time for computations. Then we pass to finding coefficients of expansion (14). Analogously to the previous case, we write the corresponding residual for the second equation of the System (6) as: where To find coefficients d (k) n , we make an orthogonal projection of the residual (20) to all functions of the basis f n N n=0 that is the main idea of the Petrov-Galerkin method [6]: Finally we obtain the following SLAE for d (k) : Appl. Sci. 2020, 10, 8373 Integrals in Equations (19) and (23) can be evaluated numerically, and differential operators are approximated with their finite-difference analogs.

R-Functions Method and Transfinite Interpolation
The RFM [10] allows all prescribed boundary conditions to be satisfied exactly at all boundary points. The R-functions are real-valued functions that behave as continuous analogs of logical Boolean functions. With R-functions, it became possible to construct functions with prescribed values and derivatives at specified locations. Furthermore, the constructed functions possess desired differential properties and may be assembled into a solution structure that is guaranteed to contain solutions to the posed boundary value problems. For example, the homogeneous Dirichlet conditions may be satisfied exactly by representing the solution as the product of two functions: (i) a real-valued R-function that takes on zero values at the boundary points; and (ii) an unknown function that allows one to satisfy (exactly or approximately) the differential equation of the problem.
A function f (x, y) is called an R-function if its sign is completely determined by the signs (but not magnitudes) of its arguments. The most popular system of R-functions is the system with R-operations (R-conjunction, R-disjunction, and R-negation) defined in the following way: The above R-functions correspond to the Boolean logic functions ∧, ∨, ¬ in a piecewise sense and allow constructing normalized implicit functions for complex-shaped geometric objects. Let the geometric domain Ω = B(Ω C , Ω H ) be constructed as a Boolean (union and intersection) combination of primitive regions Ω C , Ω H , defined by real-valued functional inequalities ω C (x, y) > 0, ω H (x, y) > 0 respectively. If f is an R-function corresponding to the Boolean function B, then the implicit function of the resulting geometric domain is immediately given by It is known that the equation of the boundary ∂Ω ( f = 0) is called normal if the value of f (x, y) is equal to the Euclidean distance from the point (x, y) to the boundary ∂Ω. Similarly a function f that coincides with the normal function only on the boundary ∂Ω is called normalized and has a property that: If both implicit functions, ω C , ω H , are normalized on the boundaries Ω C , Ω H respectively then all of the R-functions above preserve this property, and the function f (ω C , ω H ) is normalized on the whole boundary ∂Ω.
With the help of R-functions, it is possible to make the transfinite interpolation, i.e., to construct a continuous expression satisfying different conditions on boundary parts ∂Ω C , ∂Ω H . Here, the "bonding" operation for boundary conditions must be used. For example, the function: satisfies the conditions The latter expression is a generalized analog of the Lagrange interpolation formula and allows us to pass from inhomogeneous boundary conditions to homogeneous ones.
The detailed justification of the RFM with Galerkin technique in connection with solving the natural convection problem in enclosure regions is presented in [12,41,42]. In particular, the natural convection in presence of local heat is investigated in [12]. This justification is based on variational principles [43] and is appropriate for the wide class of bases composed of both spectral and compactly supported functions.
It is convenient to take normalized functions due to the fact that with them it is easier to evaluate some important differential characteristics of the solution, in particular, the Nusselt number, Nu. The local Nusselt number at a point of boundary ∂Ω is expressed as: Taking into account Equations (3) and (11) and properties of the normalized function ω, we can write that: Substituting Equation (26) into Equation (29) and taking into account properties of the normalized functions ω C , ω H , we get: Therefore on each part of the boundary ∂Ω we have the simple expressions: and Average Nusselt numbers on boundaries ∂Ω C , ∂Ω H are found by integration of local Nusselt numbers along corresponding contours:

Numerical Experiment
Consider the model boundary value problem in circular annulus (Figure 3a) with boundary conditions: The common way to obtain a function describing the boundary of the domain shown in Figure 3b is the use of R-conjunction operation: For the circular annulus, expressions for "hot" and "cold" parts of the boundary can be written as follows: Here, the normalizing terms 1/(2R i/o ) are taken to provide the normalized functions ω C , ω H whose values, as well as values of ω, are close to the distances to corresponding boundaries at points located in their vicinity.
At first glance, for such a simple domain as shown in Figure 3a, by analogy with the approach proposed in [38,39], one would take instead of Equation (34) another expression: However, this way of constructing the function ω is restricted by the narrow class of domains (eccentric circular annuli, elliptic annuli, etc.) and it is inappropriate, for example, for annular domains with more complicated inner boundaries (triangular, rectangular, et al. [44,45]). Additionally, due to the fact that the function (36) is not normalized, namely, we shall obtain more complicated formulae for the local Nusselt numbers than Equation (31) while the application of functions (35) yields the following expressions: With the help of the transfinite interpolation (26), we pass to the problem with homogeneous boundary conditions: In this work, we take the basis of multiquadric RBF (MQ-RBF) which is often used in various meshless techniques [18,19]: where X i(n) Y j(n) are the centers of functions χ n and α is their shape parameter. For simplicity we take X i(n) Y j(n) as knots of a regular two-dimensional square mesh with widths h X = h Y = α.
Instead of the latter representation we may construct 2D MQ-RBF by means of the tensor product: We used both formulae, (38) and (39), and numerical experiments demonstrated that from the point of view of computation costs and accuracy they are almost equivalent.
To define the appropriate parameters of our technique, we tested it on the problem of the stationary heat diffusion problem in the annulus cavity, whose stationary axisymmetric analytical solution is: (40) Figure 4 illustrates temperature isolines for stationary heat transfer problem in circular annulus with different ratios between inner and outer radii, which are in good accordance with the solution (40). Only 36 basic functions were taken with centers in knots of the square regular mesh. Numerical integration was realized by the two-dimensional method of trapezoids on the regular rectangular mesh of 32 × 32 nodes.   Table 1 presents the relative error as a function of the number of mesh nodes. The same resuls are shown in Figure 5 in logarithmic scale. Then we studied the stationary convective-diffusive heat transfer. For this set of parameters, we have Pr ≈ 0.67. In the case of a circular annulus, the characteristic length L in the expression for the Grashof number is equal to i R , i.e., To increase the accuracy, instead of an expression (37) for "bonding" function Φ, we can take immediately the stationary heat diffusion distribution (40): ( )  The geometrical and thermophysical parameters were as follows: Relative errors between approximate and exact solutions in L 2 -norm were ε ≈ 4 · 10 −4 for R i = R o /2.6 and ε ≈ 7 · 10 −3 for R i = R o /10. Table 1 presents the relative error as a function of the number of mesh nodes. The same resuls are shown in Figure 5 in logarithmic scale. Table 1. Relative errors between approximate and analytical solutions as a function of the number of mesh nodes along each direction (R i = R o /10).   Then we studied the stationary convective-diffusive heat transfer. For this set of parameters, we have Pr ≈ 0.67. In the case of a circular annulus, the characteristic length L in the expression for the Grashof number is equal to R i , i.e., Gr = Γβ(T H − T C )R 3 i /ν 2 . To increase the accuracy, instead of an expression (37) for "bonding" function Φ, we can take immediately the stationary heat diffusion distribution (40):

Number of Nodes
. Figures 6 and 7 illustrate temperature and streamfunction isolines for both configurations respectively. Temperature distributions along the symmetry axis OX are shown in Figure 8.   Our results were compared with a well-known asymptotic solution obtained by Hodnett [24] and Mack and Bishop [25]: Our results were compared with a well-known asymptotic solution obtained by Hodnett [24] and Mack and Bishop [25]: Here the first, "diffusive", term is defined by Equation (40) and the second, "convective", term is expressed by a rather cumbersome formula for function f with a set of coefficients presented, for example, in [24].
It observed a good correspondence between the results obtained by means of our semi-analytical approach on the base of MQ-RFM and representation (41) (within the restrictions imposed on the asymptotic solution: the inner cylinder radius R i and the temperature difference (T H − T C )/T C should be small enough). Our method demonstrated a very high rate of convergence. For example, in the case R i = R o /2.6, the relative error between approximate and analytical solutions became less than 1 · 10 −2 after only two iterations. Here, Nu ∂Ω C ≈ 242.5, Nu ∂Ω H ≈ 238.6, with relative errors between 1 · 10 −2 and 3 · 10 −2 respectively as compared with those evaluated with using expression (41) which gives us Nu ∂Ω C = Nu ∂Ω H = 245.9. The obtained temperature profiles are also in good accordance with both experimental and numerical results obtained in the work of Garraud [36]. However, in her work, to obtain an appropriate accuracy, either a uniform finite-element grid must be used with 11,000 nodes or an adaptive non-uniform grid with 2000 nodes. At the same time, our approach provides a semi-analytical solution with the same accuracy in the form of a series of 36 terms only with respect to very simple 2D MQ-RBFs (38).
Beside temperature profiles, the temperature difference between two opposite points along the sensitivity axis OX was evaluated ( Figure 9) to find an optimum relative position X * = X(∆T max ) of sensors of a TA with respect to the inner radius. The relative position is as follows: The sensitivity S was evaluated as a function of temperature difference ∆T max from external acceleration Γ = Γ(g), where g is acceleration of gravity (g ≈ 9.8 m/s 2 ) ( Figure 10). A linear part of these dependence between Γ = 0 and Γ = 1000 g corresponds to the range of "small" and "medium" accelerations and shortens as R i decreases.
sensors of a TA with respect to the inner radius. The relative position is as follows: It does not depend on Ra number and almost coincides with the theoretical estimation according to Hodnett's model (41) as well as to numerical and experimental data [36]. The sensitivity S was evaluated as a function of temperature difference max T Δ from external acceleration Γ = Γ(g), where g is acceleration of gravity ( 9.8 g ≈ m/s 2 ) ( Figure 10). A linear part of these dependence between Γ = 0 and Γ = 1000 g corresponds to the range of "small" and "medium" accelerations and shortens as i R decreases. It does not depend on Ra number and almost coincides with the theoretical estimation according to Hodnett's model (41) as well as to numerical and experimental data [36]. The same shape of the plots, with the steep linear slope in the range of "small" accelerations, maximum achieved at "medium" accelerations, and smooth deceleration at "high" accelerations, was observed experimentally [36]. In that work, the sensitivity S = 0.157 °C/g was measured in the interval between Γ = 0 and Γ = 500 g for the sensor with o 750 μm R = , 50 μm It should be noted that increase of acceleration to the values greater than some hundred g yields the turbulence phenomenon [46] which can explain partially the behavior of S(Γ) plot. However, the simplified RFM-Galerkin scheme proposed in this article has low accuracy for such case and requires further modification to simulate turbulent flows. The results shown in Figure 10 and corresponding to values Γ > 1000 g can only be considered qualitatively.

Conclusions
The new modification of the Galerkin method for solving stationary convection-diffusion problems in arbitrarily-shaped domains was proposed. It is based on the combined use of the RFM with Boolean representation of the domain boundary and the Petrov-Galerkin iteration procedure with multiquadric RBFs. Numerical experiments showed the high accuracy and high rate of convergence of the novel approach. The semi-analytical solution was obtained in a closed form of the series with respect to MQ-RBFs and it satisfied the boundary conditions exactly. The technique was applied to the well-studied benchmark problem of convection in the circular annulus, which is the simplest model of the thermal accelerometer. The obtained results were in good accordance with experimental data, numerical and asymptotic solutions. MQ-RFM can be applied directly for evaluation of thermal fields in more complicated domains and easily generalized to the case of other The same shape of the plots, with the steep linear slope in the range of "small" accelerations, maximum achieved at "medium" accelerations, and smooth deceleration at "high" accelerations, was observed experimentally [36]. In that work, the sensitivity S = 0.157 • C/g was measured in the interval between Γ = 0 and Γ = 500 g for the sensor with R o = 750 µm, R i = 50 µm. It should be noted that increase of acceleration to the values greater than some hundred g yields the turbulence phenomenon [46] which can explain partially the behavior of S(Γ) plot. However, the simplified RFM-Galerkin scheme proposed in this article has low accuracy for such case and requires further modification to simulate turbulent flows. The results shown in Figure 10 and corresponding to values Γ > 1000 g can only be considered qualitatively.

Conclusions
The new modification of the Galerkin method for solving stationary convection-diffusion problems in arbitrarily-shaped domains was proposed. It is based on the combined use of the RFM with Boolean representation of the domain boundary and the Petrov-Galerkin iteration procedure with multiquadric RBFs. Numerical experiments showed the high accuracy and high rate of convergence of the novel approach. The semi-analytical solution was obtained in a closed form of the series with respect to MQ-RBFs and it satisfied the boundary conditions exactly. The technique was applied to the well-studied benchmark problem of convection in the circular annulus, which is the simplest model of the thermal accelerometer. The obtained results were in good accordance with experimental data, numerical and asymptotic solutions. MQ-RFM can be applied directly for evaluation of thermal fields in more complicated domains and easily generalized to the case of other types of boundary conditions including mixed ones in different parts of the boundary. For analytical investigation of bandwidth, a transient study based on a combination of MQ-RFM and the Rothe method can be realized by analogy with [11]. As for solving 3D problems, it should be noted that, with the help of RFM it is possible to describe any 3D object geometry. However, we cannot directly transfer the technique described in the article to 3D problems due to difficulties in extending the vorticity-stream function formulation to the multidimensional case [47].