key: cord-0045664-92k8bwio
authors: Chistyakov, Viktor F.; Chistyakova, Elena V.; Levin, Anatoliy A.
title: Application of Underdetermined Differential Algebraic Equations to Solving One Problem from Heat Mass Transfer
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50426-7_7
sha: d8616a6279d28d5796ba89ee47dab273a1e63f6d
doc_id: 45664
cord_uid: 92k8bwio

This paper addresses a mathematical model of the boiling of subcooled liquid in an annular channel. The model is presented by a mixed system of ordinary differential equations, algebraic relations and a single partial differential equation, which, written together, can be viewed as an underdetermined differential algebraic equation with a partial differential equation attached. Using the tools of the differential algebraic equation theory, we reveal some important qualitative properties of this system, such as its existence domain, and propose a numerical method for its solution. The numerical experiments demonstrated that within the found existence domain the mathematical model adequately represents real-life boiling processes that occur in the experimental setup.

Many technical systems and processes can be described by mathematical models in the form of differential algebraic equations (DAEs). Generally speaking, a DAE is a generalization of ordinary differential equations and can be viewed as an incomplete system of differential equations closed with a set of algebraic relations. DAEs have received much attention due to a wide range applications in electric circuits simulation [1] [2] [3] , mechanics of multibody systems [4] [5] [6] , flow networks [7, 8] , etc. The complete coverage of the topic can be found, for example, in [9] [10] [11] [12] and the references therein. In this paper, we focus on a particular case of DAEs that arises as a mathematical model of boiling of subcooled liquid in a steam generating channel of an experimental unit. This model appears to be an underdetermined DAE with a partial differential equation attached to it. We present a qualitative analysis of this particular system and reveal some important properties of underdetermined DAEs. We propose a numerical method that takes into consideration all important features of the model and discuss the compliance of the results of numerical experiments with the real-life process.

We consider a mathematical model of the coolant boiling-up under the conditions of unsteady heat generation on the wall of the flow-through duct, which detailed description can be found in [13] , where the model was addressed by solving all equations individually. However, a closer look at a set of equations that comprise the model shows that it can be written as a system of the following form 

where u(t) is the desired vector-function with the following components: u 1 (t) is the temperature of the surface on which the unsteady heat transfer crisis develops; u 2 (t) is the thickness of the liquid microlayer, which evaporation dynamics determines heat transfer; u 3 (t) and u 4 (t) are the interfacial velocities in the upper and lower part of the channel, respectively; u 5 (t) is the vapor layer thickness which is determined by the interfacial velocity; u 6 (t) is the temperature inside the vapor layer; u 7 (t) is the vapor pressure inside the vapor layer; u 8 (t) is the liquid temperature gradient at the boundary with a cooling metal surface:

∂T ∂z |z=0 , where z is a coordinate axis perpendicular to the heating surface; k j , j = 1, 7, A kof , T s , p s , R v , L, p 1 , p 2 are some given parameters (their detailed description and values can be found in [13] ). All units of measurement are given in the SI metric system.

The initial conditions are also given

(2) One of the issues in the problem statement is that the domain [0, ϑ] is unknown to us and should be determined on the basis of numerical calculations.

The DAE (1) is underdetermined and in [13] the authors of the model used a solution to a partial differential equation to obtain the necessary closure for the system

However, from a mathematical viewpoint, such underdetermined DAEs require special research. If we could not express the function ∂T ∂z |z=0 via the other seven components of the vector-function u(t), we would have to rely only on experimental observations. The practical meaning of studying the model (1) is ensued by the recently discovered heat transfer enhancement under conditions of the rapidly oscillating interface. It was shown in [13] [14] [15] [16] that this enhancement is mainly associated with the short-term presence of a metastable liquid on a metal surface. The direct experimental observation and analysis of the dynamics of such processes are limited in view of the smallness of time periods and the geometric dimensions of the objects. Therefore, one of the most important tools in this area is the numerical integration of the equations that comprise the mathematical model.

The DAE (1) is underdetermined in the sense that the number of equations exceeds the dimension of the desired vector-function u(t). Such systems can be generalized as

where A(u, t) is an (ν × n)-matrix, B(u, t) is an ν-dimensional vector-function, u = du(t)/dt, u ∈ U ⊆ R n , and the following condition holds

Usually, we also have a set of initial data

We consider that U = {u : a − u ≤ ∈ R 1 , > 0}. Below we assume that all entries are sufficiently smooth. One of the most important characteristics of DAEs is index. An index is a way of measuring the distance from a DAE to its related ordinary differential equation (ODE). The concept of index for DAEs has several definitions, most of which are extensively covered in [12] and [18] . However, each new type of DAEs often requires a refinement of one of the popular notations.

where

where rank A l (a, α) = min (n, ν). Then, (6) it is said to be the left regularizing operator for the DAE (3) and the smallest possible i is said to be the index of the DAE (3).

Now we will show one of the ways of building the LRO. Let there exist a (ν × ν)-matrix P(u, t) ∈ C 1 (U × T ) with the properties: det P(u, t) = 0∀(u, t) ∈ C 1 (U × T ),

where the zero block has the dimension ([ν−r]×n), r = max{rankA(u, t), (u, t) ∈ U × T }. Differentiate with respect to t the algebraic relations in (8) . We obtain

Let rank A 1 (a, α) = min (n, ν). Then the system

can be rewritten as (10) where the block A 1 11 (u 1 , u 2 , t) has the dimension (ν × ν), det A 1 11 (a, α) = 0, u = (u 1 , u 2 ) . Set u 2 = ψ(t), where ψ(t) is an arbitrary smooth function, ψ(α) = a 2 , (u 1 (α), u 2 (α)) ) = (a 1 , a 2 ) = a . Then, taking into consideration that the matrix A 1 11 (u 1 , u 2 , t) is continuous, we can obtain a system of regular ODEs in the neighborhood of the point (a, α):

The reasoning presented above can be formulated as the following theorem.

If rank A(a, α) = rank(A(a, α)|B(a, α)), then (8) satisfies the condition B 2 (a, α) = 0, and there exists a segment [α, α + ε] ⊆ T , on which the solution (u 1 (t), ψ (t)) to the initial problem (3), (5) can be found.

,

If in (9) rank A 1 (u, t) < min (n, ν) ∀(u, t) ∈ U × T, then, similarly, we can build an operator

Suppose that after l steps of such a process we arrive at

where rank A l (a, α) = min (n, ν) and

If the conditions

are satisfied, we can prove solvability of the initial problem (3), (5) . The number l is said to be a rough index of the DAE (3).

If the matrices that define operators Ω i in (13) are sufficiently smooth, then the product l−1 i=0 Ω i can be represented in the form of the operatorΛ l (u).

If a DAE is closed, then the notion of index for underdetermined DAEs introduced in this paragraph fully coincides with the LRO based index from [17] , which was further developed in [20] and [19] . (13) can be built differently (see, for example, [12] ). Example 1. Consider the following DAE:

The rough index for this DAE is not determined, because the matrix A 1 (u, t) in (9) is singular on the solutions of the system. However, the LRO exists and can be written as

and the index of the system is 1 in terms of Definition 1. Indeed,

Let the DAE have the form:

In this case, the LRO is defined and

so the system has index 2 in terms of Definition 1. It can be readily seen that

It can be easily verified that this system has a rough index and it is equal to 2.

If in Example 1 we set u 3 as a given smooth function ψ(t), we obtain closed systems with unique solutions. Technically, these systems have index 1: the LRO is of the first order, but if we take into account the derivatives that participate in the LRO's coefficients, the rough index coincides with the index in terms of Definition 1.

If we differentiate (1) with respect to t, we can see that its index is equal to 1 in terms of Definition 1. Therefore, (1) has a unique local solution in the neighborhood of the initial data. To find out if the system has a solution on a given segment, investigate its Lyapunov stability. Find the stationary solution of (1) by solving the following system of non-linear equations

From the equations B 1 (u) = 0, B 3 (u) = 0, B 7 (u) = 0 and the expression for u 8 , we obtain

Now multiply the first equation of the system by k 2 /k 1 and subtract it from the second equation. We arrive aṫ

where c is an arbitrary constant. For calculate c, we can use the stationary solution of the system. Taking into consideration the values for the initial data and coefficients from [13] , we have

The component u 2 (t) is present in the denominator of the right-hand part of the system: it participates in B 1 (u(t)), B 2 (u(t)), B 7 (u(t)). At some point t = ϑ, u 2 (t) turns into zero, since u 2 (0) = 0.0001 > 0. Therefore, we can find the segment [0, ϑ], on which the solution to the initial problem (1), (2) exists. However, the system is quite stiff. We have the following values for the derivatives of the solution's components in the starting point: , and the eigenvalues of the Jacobi matrix are of the same order. This fact preconditions the choice of the numerical method.

Due to the fact that the right-hand part of (1) has only the first continuous derivative, we have to choose among the first order methods. Therefore, to solve (1) numerically, we applied the implicit Euler method:

where N is a number of integration steps, is a square matrix that was obtained from (1) by excluding a zero column from the matrix A, t * is derived from the empirical data. We solve the nonlinear system (14) for each i by the Newton method taking only the first iteration. As the result, the calculation algorithm has the form:

Since (1) has index 1, the matrix [Â+τ J(w i )] is non-singular for a sufficiently small τ [17] and the method has the first convergence order. The integration step was selected taking into account the eigenvalues of the Jacobi matrix J(w 0 ). The solution to the initial value problem (1), (2) in a mathematical sense exists only on the semi-interval [0, ϑ). However, numerical experiments showed that we can apply the difference scheme (15) beyond ϑ. Numerical calculations also revealed that when u 2 (t) reaches a zero-crossing, the solution surges up, which is fully consistent with the physical experiment ( Fig. 1) : at some moments, the liquid microlayer described by the component u 2 (t) completely evaporates and gets refilled.

In particular, the performed numerical calculations allowed us to define the semi-interval [0, ϑ) = [0, 0.02] millisecond. This result coincides with the empirical findings from [13] .

We presented a qualitative and numerical study of the model for the boiling of subcooled liquid in an annular channel. We obtained the solvability conditions and proposed a numerical method of solution. The analysis of the mathematical model (1) was performed under the assumption that the elements of the matrix A and a number of other given parameters are constant. In fact, their values depend on the solution to the system and are located within a certain range determined by the physical meaning of the problem. Further research implies refinement of the model and investigation of its properties to determine the possible ranges for variable parameters.

CAD-based electric-circuit modeling in industry I: mathematical structure and index of network equations

Structural analysis of electric circuits and consequences for MNA

Differential-Algebraic Systems: Analytical Aspects and Circuit Applications

Nonholonomic Motion of Rigid Mechanical Systems from a DAE Viewpoint

Numerical Methods in Multibody Systems

Numerical methods in vehicle system dynamics: state of the art and current developments

A unified (P)DAE modeling approach for flow networks

Reduced-Order Modeling (ROM) for Simulation and Optimization

Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations

Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems

Differential-Algebraic Equations: Analysis and Numerical Solution

Differential-Algebraic Equations: A Projector Based Analysis

Self-oscillatory regime of boiling of a highly subcooled liquid in a flow-passage annular duct

Self-excited pressure pulsations in ethanol under heater subcooling

Heat transfer at cooling high-temperature bodies in subcooled liquids

Subcooled water quenching on a super-hydrophilic surface under atmospheric pressure

Algebro-differentsial'nye operatory s konechnomernym yadrom (Algebraic Differential Operators with a Finite-Dimensional Kernel)

Index concepts for differential-algebraic equations

Solution of differential algebraic equations with the Fredholm operator by the least squares method

Evaluation of the index and singular points of linear differential-algebraic equations of higher order