Padé approximants of canards and critical regimes of Darrieus wind turbine model

1.   Introduction

When the first publications by French mathematicians appeared several decades ago describing and investigating a new mathematical object called canards, many mathematicians, including one of the authors of this paper, were convinced that canards (or French ducks) were hardly interesting from an applied point of view because of their exceptionally high sensitivity to parameter changes. However, after a relatively short period, the canard apparatus has become an effective means of investigating critical phenomena of various natures. The theoretical aspects of the canards technique are reflected, for example, in the ^{[14,27,28,42,47]}. In applications, the canards are used in two roles. The use of periodic canards is based on the phenomenon of the so-called canard explosion ^[26,27,42], the main feature of which is a sharp increase in the amplitude of oscillations at a slight change in one of the parameters. Non-periodic canards serve as a model watershed between processes with fundamentally different characteristics. For example, in combustion processes, they allow us to find the critical conditions of thermal explosion that separate the conditions of slow safe combustion from the explosive combustion reaction ^{[2,18,19,20,21,22,23,34,35,37,38,39,42,44,46]}.

Restricting ourselves to applications in the field of engineering mathematics, we can single out works related to the study of the dynamics of chemical reactors ^{[10,15,16,17,31,36]} and combustion processes ^{[2,18,19,20,21,22,23,34,35,37,38,39,42,44,46]}, rotating machinery ^[6], two-wheel vehicle ^[45], metal structures at ultra-low temperatures ^[8], aircraft ground dynamics ^[32,33], laser dynamics and optical amplifiers ^{[24,30,40,41]}.

One of the central problems arising in using the canards theory in solving specific applied problems is the problem of computing the so-called canard parameter values. The main difficulty lies in the fact that it is necessary to determine these values with a very high degree of accuracy. As it is shown in ^[4], for the classical van der Pol system at $\varepsilon = 0.01$ it is necessary to carry out calculations with accuracy to the tenth decimal place, inclusive. This means that when using the asymptotic expansions for the canard values in the Maclaurin form, it is desirable to find as many terms of the expansion as possible. At the same time, however, the complexity of analytical calculations increases rapidly, rendering the practical implementation of theoretical concepts challenging. A well-established solution to this issue relies on the use of Padé approximations, which enable a significant improvement in the accuracy of asymptotic formulas without requiring additional cumbersome analytical computations ^[1,3].

Although Padé approximations ^[3] were discovered much earlier than canards, they are not at all more common in the applied mathematics literature. Therefore, we will give the necessary information without delving into theoretical problems.

Let some function $f$ of $\varepsilon$ can be represented of form

where

and

The rational function $[k/m] = \frac{A(\varepsilon)}{B(\varepsilon)}$ is the Padé approximant if

Here $(k+m = n)$.

Equating powers of $\varepsilon$ in the equality

we get the formulas for $a_0, a_1, \dots, a_k$ and equations for $b_1, \dots, b_m$.

If we restrict ourselves to the case $n\leq 4$, we obtain the following equations and the corresponding diagonal, i.e., $k = m$, Padé approximants. Thus, we have

where

can be found from the equalities

Analogously, we have

where

can be found from the equalities

2.   Padé approximants for van der Pol system

Consider the well-known van der Pol system

with $F(x) = x^3/3-x$. As it was shown in ^[47] for some special value $\alpha^* = \alpha(\varepsilon)$ there exists the periodic canard and the following asymptotic expansion takes place

i.e., $f_0 = 1, f_1 = - 1/8, f_2 = -3/32, f_3 = -173/1024, f_4 = -7593/16384$. Note that this formula contains the term $-7593 \varepsilon^4/16384$, derived by the first author, who refined the canard value.

Thus, the following Padé approximants for $\alpha^*$ take place

and

In particular, for $\varepsilon = 0.01$

here (2/3) means that 2 corresponds to the canard with a head while 3 corresponds to the canard without the head. Let $\check\alpha = 0.998 740 4512$, and $\hat\alpha = 0.998 740 451 3$, then any value of $\alpha\in [\check \alpha, \hat\alpha]$ corresponds to the canard. Further,

and

This means that for the van der Pol equation, the Padé approximation $[2/2]$ gives precisely the canard value, which within the framework of the method can be considered as precisely as it belongs to the interval of canard $[\check \alpha; \hat \alpha]$ described above. Note that the values of $\alpha_{M4}$ do not belong to this interval.

Note that in the theory of canards, there is a statement that is formulated as follows: "The life of ducks is short". The meaning of this statement is that the transition from a small cycle to a relaxation cycle corresponds to a change in the parameter $\alpha$ over an interval of the order of $exp(-1/c \varepsilon), c > 0$. Thus, within an exponentially small range of the parameter, the very fast transition occurs from a small amplitude limit cycle via canard cycles to a large amplitude relaxation cycle. This very fast transition is called a canard explosion ^[26,27,42].

3.   Darrieus wind turbine model

In this section, a mathematical model of a small vertical-axis wind turbine known as a Darrieus wind turbine is presented. The mathematical model of the wind turbine is a three-scale singularly perturbed differential system. In this model, the aerodynamic moment is approximated by polynomials based on experimental data. The conditions for the occurrence of dangerous large amplitude oscillations, which are modeled by canard trajectories, were found, which made it possible to find the critical values of the model parameters. A mathematical model of a Darrieus wind turbine installation is considered, which is a three-speed differential system ^[11,12,25]:

Here $\Omega$ is the turbine rotational speed, $I$ is the current in the armature winding, $R$ is the external resistance, $M$ is the relative moment of aerodynamic forces, $J$ is the moment of inertia of the turbine, $L$ is the armature inductance. The remaining constants are: $k$ is the electromechanical interaction coefficient, and $r$ is the small internal resistance of the armature. In the system under consideration, the small parameters are the quantities $L$ and $\varepsilon$, which in turn means that the external resistance $R$ is the slowest variable, and the current variable $I$ is the fastest of the three variables.

The quantities in the Figure 1 are dimensionless ($M = \frac{M_a}{0.5 \rho SbV^2}$, $\Omega = \frac{b\omega}{V}$, where $\omega$ is the angular velocity, $b$ is the distance from the effective pressure of the blades to the axis of rotation, $V$ is the air velocity, $M_a$ is the moment of aerodynamic forces, $\rho$ is the air density, $S$ is the blade area) ^[11,12]. The coefficients are chosen to achieve agreement with the experimental data ^[5]. The experimental data were obtained by blowing a Darrieus wind turbine with a NASA 0012 profile ^[5], see also ^[11,12,25].

Based on the experimental data for the interval $[0, 8.4]$, it is convenient to choose a fifth degree polynomial as an approximation of the aerodynamic moment

This polynomial can be smoothly extended to the interval $[0, 8.89]$. The role of such continuation in the Figure 1 is played by the following fourth degree polynomial

The graph of the resulting function $M(\Omega)$ is presented in Figure 1.

Note that the moment function presented in the article approximates the available experimental data over a known interval $[0, 8.89]$. Therefore, it makes sense to consider her behavior only in this interval. The behavior beyond this interval depends on the choice of a specific approximating function. Different approximating functions can be arbitrarily close in the selected interval but differ significantly beyond it.

Since the small parameter $L$ is multiplied by the derivative of the current strength, system (3.1) is singularly perturbed. System (3.1) has a two-dimensional invariant manifold ^[9,43]

The function $h$ can be found from the invariance equation:

Neglecting terms of order $O(L)$ since $L \ll \varepsilon$, we obtain:

Movement along an invariant manifold is described by a differential system:

The resulting system is slow/fast since the right side of the equation for the slow variable $R$ is multiplied by a small parameter $\varepsilon$. Based on this, to analyze it, one can apply the apparatus of the theory of relaxation oscillations ^[29] and the theory of canards (see, for example, ^[26,27,42]). The slow curve of system (3.2) is given by the equation:

From (3.3), we can express $R$ and get the equation of the slow curve in the explicit form:

The expression obtained here from the invariance equation in (3.3) equals the critical manifold that is found for $L = 0$ in (3.1), i.e., if the flow on the fastest scale is assumed to have reached a quasi-steady state.

The plot of the slow curve is shown in Figure 2, at $k = 0.5, r = 0.1$. Let us find the derivative of the function $R(\Omega)$ to study the stability of the slow curve (3.4). Direct analysis shows that the stable parts of the slow curve are:

and the unstable part corresponds

Points $A (1.01265035468611, 5.66632535089612)$ and $B (3.65182989752896, 0.593912765877734)$ are the fold points of the slow curve specified by Eq (3.4) and correspond to the points of its extrema. In addition, they are points of change in stability.

Let $F(\Omega, I, R) = \varepsilon(\alpha -\Omega)$. Then, setting $\Omega = x,$ $R+r = y$ and $J = 1$ we obtain the differential system which was analyzed above with

Depending on changes in the value of the parameter $\alpha$ the equilibrium position moves along the slow curve. If the equilibrium position lies on the stable part of the slow curve, then it is asymptotically stable. From a physical point of view, this situation seems preferable, since it corresponds to the stationary operating mode of the system. The case when the equilibrium position is in an unstable area and a relaxation cycle occurs seems dangerous from a physical point of view since oscillations with a sufficiently large amplitude arise in the system (Figure 3). Between these modes (stable equilibrium position and relaxation cycle) there are transition modes. When the value of the parameter $\alpha$ changes, the singular point passes the extremum point (fold points) and merges with it, while losing stability and becoming unstable, which corresponds to the Andronov-Hopf bifurcation. With a further change in the value of the parameter $\alpha$, the equilibrium position moves beyond the stall point to the unstable part of the slow curve, while remaining in a small neighborhood of the stall point, on the order of $O(\varepsilon)$ for $\varepsilon > 0$, the size of the cycle begins to increase (Figure 4a,b). The critical value of the parameter $\alpha$ corresponds to the fold points. At a critical value of the parameter $\alpha$, the cycle becomes a canard trajectory, and then a canard explosion occurs ^[26,27,42]. It should be noted that this mathematical model (8.1) is extremely sensitive to changes in the parameter $\alpha$ in the vicinity of critical values. Thus, with a slight change in the parameter $\alpha$, a very rapid transition occurs from the limit cycle of small amplitude to the relaxation cycle (Figures 3 and 4).

Substituting the selected numerical values of the parameters into the corresponding formulas of the Section Appendix gives the following results

Taking into account the Maclaurin series

we get the coefficients of Padé approximants $[1/1]$ of form

and consequently,

Comparing the values $\alpha_{M2}$ and $[1/1]$ with the results of numerical experiments we can see the advantage of the Padé approximants

since the Padé approximant gives an order of magnitude better accuracy.

4.   Conclusions

At least two natural questions have been left unanswered in this paper. The first one concerns the necessity of a more accurate construction of a two-dimensional invariant manifold when reducing a three-dimensional wind turbine model to a two-dimensional one. In this paper, it is implicitly assumed that $L = o(\varepsilon^2)$. An increase in accuracy when computing the critical values of the control parameter requires an increase in accuracy when performing the model reduction. Fortunately, the methods that allow us to do this are well-developed and quite effective ^[9,43]. Since $\frac{\partial}{\partial I}(k\Omega-(r+R)I) = -(r+R)$, and $R\ge0, r > 0$, the critical manifold is normally attractive.

The second issue is related to the dependence of the results on the choice of the function approximating the experimental data on the aerodynamic moment. The results of the preliminary analysis performed by the first author of the paper show that there are no fundamental differences even if a piecewise linear approximation is used. It should be noted that, despite significant differences in the approximating functions, the recommended critical values of the parameter $\alpha$ differ only by the second decimal place. Nevertheless, the authors believe it is necessary to revisit these issues in future studies.

The problem of the application of Padé approximations to calculate asymptotic expressions of the values of the control parameter at which the well-known phenomenon of canard explosion is observed is considered. Examples demonstrating the advantage of Padé asymptotics over asymptotic Maclaurin expansions are proposed. Asymptotic formulas of the sufficiently general form have been derived for certain classes of differential systems, distinct from the van der Pol systems. The obtained mathematical results are applied to the analysis of the mathematical model of the wind power plant with the vertical axis of rotation, known as the Darrieus wind generator. The polynomial approximation of the aerodynamic moment made it possible to find the conditions for the occurrence of dangerous large-amplitude oscillations, which are modeled by canards, which made it possible to find the critical values of the model parameters. An analysis of the influence of the choice of approximation of the relative moment of aerodynamic forces on the dynamics of the system is carried out.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank the anonymous referees and the handling editor for their valuable comments and suggestions which have led to a significant improvement to the paper.

Conflict of interest

The authors declare no conflicts of interest.

Appendix

Most of the works devoted to duck trajectories consider van der Pol systems, i.e., the system ^[47]:

Consider an auxiliary differential system on the plane which are not van der Pol type systems

Suppose that for the system under consideration, the conditions under which the canard explosion occurs are satisfied (see, for example, ^[13,14,28] and references therein). In particular, function $y = \varphi_0 (x) = \frac{k^2x}{p(x)}$ on some segment $[x_1, x_2]$ containing the point $x_0$ is sufficiently smooth and $\varphi'_0 (x_0) = 0, \, \, \varphi''_0 (x_0)\neq 0.$ Note that

Following to ^[47] the canard can be found as an asymptotic expansion

Our goal is to obtain the asymptotic expansion for the corresponding canard value of $\alpha$.

The invariance equation is

or

Setting

we obtain the invariance equation of form

Equating the powers of $\varepsilon$ we obtain

$\varepsilon^0:$

which implies

$\varepsilon^1:$

which implies

This means that

$\varepsilon^2:$

which implies

and, therefore,

$\varepsilon^k:$ Let

for $k > 1$, then

which implies

and

Note that it is possible to formalize the derivation of the asymptotic expansion for the parameter $\alpha$ by formulating a proposition and a corresponding proof. The given reasoning can be made sufficiently rigorous by the mathematical induction method following to ^[47] but that is beyond the scope of this article.