Reinforcement learning-based optimization of locomotion controller using multiple coupled CPG oscillators for elongated undulating fin propulsion

1.   Introduction

Aeroengine is an important equipment to propel the aircraft forward and its total thrust is the sum of the thrust generated by the core engine and the turbine fan. The air flow state between adjacent blades of the turbine fan determines the thrust value of aeroengine. Since the experiments cost a lot, numerical simulations become an essential part of the blade geometry design and optimization [5,10,2,12,7]. Nevertheless, many difficulties arise in numerical simulation such as nonlinearity, high Reynolds number, complex three-dimensional (3D) geometrical domains, and boundary layer effect (see, e.g., [9]). In order to alleviate these difficulties, we will use the dimension splitting method to simulate the airflow state of aeroengine blade fan.

Based on the article [6], we establish a semi-geodesic coordinate system (called R-coordinate system), whose two basis vectors are on the manifold, and the other is along the hub circle. Thus, Navier-Stokes equations (NSEs) in the R-coordinate system can be rewritten as a set of membrane operator equations on the blade surface, and the bending operator equations along the hub circle [8]. By using Euler central difference scheme to approximate the third variable, the 3D NSEs become a series of two-dimensional (2D) equations with three variables. After successively iterations, the approximate solution to the NSEs can be obtained. Obviously, the significant feature of this new method is this method only solves the 2D problem in each sub-domain. In addition, it can alleviate the boundary layer effect by approaching adjacent surfaces, and the parameterized surface provides convenience for blade design and optimization.

The purpose of this work is to introduce our proposed method to simulate the flow state of the channel of aeroengine turbine fan. A lot of work has been carried out in the field of viscous flow and its applications in aeroengine turbine [16,15,4,11,14,13,3]. In this paper, a toy model is designed to give a comparison between our novel method and traditional methods. It turns out the new method shows a good performance for the toy model. Then we apply the proposed method to simulate the flow state of aeroengine turbine fan.

The present paper is built up as follows. In Section 2, some essential differential geometry knowledge is briefly introduced, then the R-coordinate system is established. Meanwhile, the NSEs' new form, splitting method and the variational form are formulated in Sections 3. Furthermore, we derive the finite element form in Section 4. Section 5 presents the numerical results, which contain the comparison of the new method and traditional methods, and the simulation results of aeroengine turbine fan.

2.   The model problem and the new coordinate system

In this section, we establish a new coordinate system according to the geometric shape of the blade and give the relationship between the new coordinate system and the rectangular coordinate system and the cylindrical coordinate system. To express concisely and clearly, we let Greek letters $\alpha, \beta, \cdots$ and Latin letters $i, j, \cdots$ range over the values $\{1, 2 \}$ and $\{1, 2, 3 \}$, respectively. Einstein summation convention is adopted in tensor analysis in the sequel.

Compared with the size of the aircraft engine fan, the thickness of the fan blade can be neglected. Thus the blade surface $\Im$ is considered as a 2D surface in $\mathbb{R}^3$ in this article, which is a connected subset $D$ defined in $\mathbb{R}^2$ and mapped to the range $\boldsymbol\Re(D)$ by the injective map $\boldsymbol\Re$ into $\mathbb{R}^3$. Suppose $\boldsymbol\Re$ is smooth enough, any point ${\boldsymbol x} = (r, z) \in\overline{D}$ in the Gaussian coordinate system on the surface $\boldsymbol\Re ({\boldsymbol x})$ can be expressed as in the cylindrical coordinate system

Here $\Theta \in C^2(D)$ is a smooth function of radian, $(\boldsymbol e_r, \boldsymbol e_\theta, \boldsymbol k)$ are basis vectors of cylindrical coordinate system on the fan.

The channel $\Omega_\varepsilon$ is determined by the boundary $\partial \Omega_\varepsilon = \Gamma_{in}\cup\Gamma_{out}\cup\Gamma_t\cup\Gamma_b\cup \Im_+\cup \Im_{-}$, where $\Gamma_{in}, \, \Gamma_{out}$ are inlet and outlet, $\Gamma_t, \, \Gamma_b$ are top surface (shroud) and bottom surface (hub), and $\Im_+, \, \Im_{-}$ are positive pressure surface and negative pressure surface (see Figure 1 and Figure 2). Denote $N_b$ as the number of blades and $\varepsilon = \pi/N_b$, then rotating one blade $2\varepsilon$ can get another which means there exists a family of single-parameter surfaces $\Im_\xi$ covering the flow channel by mapping $\boldsymbol\Re({\boldsymbol x};\xi):\ D\rightarrow \Im_\xi$, that is

where $\theta = \varepsilon\xi+ \Theta(r, z)$ is the rotation angle.

Let

it is clear that the Jacobian matrix $J(\frac{\partial (r, \theta, z)}{\partial(x^1, x^2, \xi)}) = \varepsilon$ is nonsingular. Thus we establish a new curvilinear coordinate system (called R-coordinate system) $(x^1, x^2, \xi)$

In the R-coordinate system, the fixed region $\Omega = \{(x^1, x^2, \xi)|\, (x^1, x^2)\in D, -1\leq\xi\leq 1\}$ is mapped to a channel $\Omega_\varepsilon = \{\boldsymbol \Re(x^1, x^2, \xi) = x^2\boldsymbol e_r+x^2(\varepsilon\xi+\Theta(x^1, x^2))\boldsymbol e_\theta+x^1\boldsymbol k, \ \forall (x^1, x^2, \xi)\in \Omega \}$.

Since the new coordinate is established, these basis vectors become a bridge to communicate different coordinate systems. Let $(x, y, z)$ and $(r, \theta, z)$ denote Cartesian coordinate system and cylindrical coordinate system, respectively, which are fixed on impeller with angular velocity ${\boldsymbol \omega}$ in a 3D Euclidean space. Let $(\textbf{i}, \, \textbf{j}, \, \textbf{k})$ and $(\textbf{e}_r, \, \textbf{e}_\theta, \, \textbf{k})$ be the basis vectors of $(x, y, z)$ and $(r, \theta, z)$, respectively, which satisfy

The basis vectors of the R-coordinate system are denoted as $(\textbf{e}_1, \; \textbf{e}_2, \; \textbf{e}_3)$, which are defined as

where

By straightforward calculation, the relationship between base vectors is as follows

where $\Theta_\alpha = \frac{\partial\Theta}{\partial x^\alpha}.$

The metric tensor $a_{\alpha\beta}$ of the surface $\Im_\xi$ is defined as

and it is easy to see that $a_{\alpha\beta}$ is nonsingular and independent of $\xi$ as follows, that is

where $|\nabla\Theta|^2 = \Theta_1^2+\Theta^2_2$.

Similarly, the covariant and contravariant components of metric tensor are $g_{ij}$ and $g^{ij}$, which are defined as

By calculation, they can be expressed as (cf.[5])

As described in our previous article [6], for given $\xi$, both surfaces $\Im_\xi$ and $\Im$ have the same geometric characteristics and interested readers can refer to the article [6].

3.   The NSEs' new form and its variational formulation

3.1.   Navier-Stokes equations in the R-coordinate system and its dimension split method

In this section, we derive the flows state governed by the incompressible rotational NSEs through employing differential operators in the R-coordinate system. Because of the lighter mass of air, its body force is neglected and the rotational NSEs can be expressed as

For a time-dependent problem, we usually discretize the problem in time and solve the static problem at each time step, then the time-dependent problem is reduced to a static problem at each time step. If we consider an implicit time discretization (e.g., the backward Euler method) to equation (3.1) with time-step size $\kappa$, and multiply $\kappa$ to it, then we have a static problem

where ${\boldsymbol f}_u$ is the velocity value of previous time step.

Assume that $\Gamma^i_{jk}$ and $\nabla_ju^k$ are the Christoffel symbol and the covariant derivative in space, while $\stackrel{\ast}{\Gamma^\alpha}_{\beta\sigma}$ and $\stackrel{\ast}{\nabla}_\beta u^\alpha$ are the Christoffel symbol and the covariant derivative on surface, their expression can be obtained by

where $\textbf{e}_{ij} = \partial_i\textbf{e}_j$ is the first-order partial derivative of basis vectors.

Lemma 3.1. [6] By denoting $(\boldsymbol e_\alpha, \boldsymbol e_3)$ as the basis vectors in the coordinate system $(x^1, x^2, \xi)$, the parts of NSEs in the new coordinate system have the following conclusions.

1). The Laplace operator $\Delta {\boldsymbol u}: = g^{ij}\nabla_i\nabla_j {\boldsymbol u}$ can be written as

where

and $\widetilde{\Delta}\Theta = \delta^{\alpha\beta}\Theta_{\alpha\beta}, \; \widetilde{\Delta} u^\alpha = \delta^{\lambda\sigma}\stackrel{\ast}{\nabla}_\lambda\stackrel{\ast}{\nabla}_\sigma u^\alpha$.

2). Its related items coriolis force and centrifugal force satisfy, respectively,

where $C^1_j = 0, \; C^2_\alpha = -2\omega r\Theta_\alpha, \; C^2_3 = -2\varepsilon\omega r, \; C^3_\alpha = 2\omega(r\varepsilon)^{-1}a_{2\alpha}, \; C^3_3 = 2\omega r\Theta_2.$

3). The pressure gradient is $\nabla p = g^{ij}\frac{\partial p}{\partial x^j}$, which can be represented as

4). The nonlinear term in equation (3.1) is

where

5). The mass conservation formula can be rewritten as

By Lemma 3.1, the incompressible rotational NSEs in the R-coordinate system can be rewritten as

where

And its channel region and boundaries are

The initial and boundary value conditions are

3.2.   Variational-difference formulation of Navier-Stokes function in the R-coordinate system

In this section, the finite-element-difference method will be presented. The first step is to divide the interval $[-1, +1]$ into N-subintervals with step size $\tau = \frac{2}{N}$, that is

Then the domain $\Omega$ is split into N-layers $\Omega = \cup_{k = 0}^{N-1}\{D\times\{\xi_k, \xi_{k+1}\}\}$ and the central difference is used to replace the derivative with respect to the variable $\xi$, i.e.,

where $w_k: = w({\boldsymbol x}, \xi_k)$. In order to simplify the equation, we denote that

Plugging (3.13) into (3.9), we obtain

where

And its boundary conditions are

where $\gamma_0 = (\Gamma_t\cup\Gamma_b)\cap\{\xi = \xi_k\}, \; \gamma_{in} = \Gamma_{in}\cap\{\xi = \xi_k\}, \; \gamma_{out} = \Gamma_{out}\cap\{\xi = \xi_k\}.$

Introduce the Hilbert space $V(D)$ by

then its inner product and norms are given, respectively, by

Without the ambiguity, the index `D' is often omitted.

For clarity and simplicity, we denote

Meanwhile, we denote ${\hat{F}^i_k} = \widetilde{F}^i_k-{\mathcal{S}}^i-{\mathcal{P}}^i$, then the 2D-3C NSEs on manifold $\Im_\xi$ can be written as

Thus, the variational problem corresponding to the boundary value problem of the NSEs is as follows:

where the linear, bilinear and trilinear forms are given, respectively,

4.   Finite element solution of 2D-3C equation

In this section, we apply the Taylor-Hood elements, i.e., $(P_2, P_1)$ Lagrange finite elements for the pair $({\boldsymbol u}, p)$. Let $V_h$ and $M_h$ be the finite element subspaces corresponding to space $V(D)$ and $L^2(D)$, respectively, which are

The product space $Y_h = V_h\times M_h$ is subspace of $Y = V(D)\times L^2(D)$ obviously. Then the variational problem (3.19) approximated by the standard Galerkin finite element method is

Suppose the finite element basis functions are denoted as

where $NG1$ and $NG2$ are the total number of nodes, respectively. Finite element expansion of $\boldsymbol w_h, \; p_h$ are

Assume solution vector $W^m$ and test vector $V^k$ are

Substitute (4.4) into (4.2), we obtain the result on 2D membrane operator:

where

Proof. In this proof, we only prove $(L(\boldsymbol w_h), \boldsymbol v_h)$. Noting that

we obtain

The remainder of the argument is analogous and is left to the reader.

Next, we will give the discrete scheme on the right side of the formula (4.2). According to the definition of inner product in space $V(D)$, this discrete scheme can be written as

Thus, we obtain

where ${\mathcal{F}}^\beta_j = \widehat{F}_j^\beta-\widehat{H}_j^\beta.$ Similarly, the bending operator can be discreted as

where

And its right terms are

It is easy to verify it, so its proof will not be given here. Then the bending operators becomes

where ${\mathcal{F}}^3_j = \widehat{F}_j^3-\widehat{H}_j^3.$

Finally, we obtain the finite element algebraic equations of equation (3.19), i.e.,

where

5.   Numerical simulations and discussions

In this section, for certain examples, the results of the new algorithm program and the traditional algorithm program will be compared to verify the accuracy of the new algorithm program. And the new algorithm program shows good performance. Then, based on the new algorithm, the flow state of the gas in the fan channel of an aeroengine is given.

5.1.   The comparison of new method and traditional method

In this part, a simple model is provided to give the comparison between results obtained by the dimension splitting method (DS method) and traditional 3D method (T3D method). In this example, we adopt stationary model and assume $\Theta = 0$ and $\theta = [-7.5^\circ, 7.5^\circ]$, see Figure 3. In this case, the central surface of both methods is completely coincident. Figure 4 presents one blade of this model and it is also the shape of the central surface. DS method is concerned with the solution of 2D surface. Figure 5(a) and Figure 5(b) show the mesh of central surface, respectively. By comparison, we can see that the mesh of Figure 5(a) is more regular than Figure 5(b). That is because Figure 5(a) is one of the solving planes while Figure 5(b) is the projection of a 3D mesh on the central plane.

51 2D-manifolds and 4987 elements per surface are used to partition the channel and for DS method while 254867 elements are used for T3D method. In this model, we assume $\boldsymbol v_{in} = 1$m/s and $\omega = 10 \ $rad/s. Figure 6 – Figure 9 show the comparison of results of the velocity distribution on the central surface. And Figure 10 shows the comparison of results of the pressure distribution on the central surface. By comparison, we know that the results obtained by 3D method and DS method are almost identical.

5.2.   Simulation results of aeroengine turbine fan

As mentioned above, the flow state determines the maximum thrust that aeroengine turbine fan can provide. Our main purpose is using the DS method to simulate the flow state of the aeroengine turbine fan. In addition to the alleviation of boundary layer effects and make parallelism easier, DS method can give one clear design objective $\Theta$ for blade shape. Because the flow state mainly depends on the shape of blades, $\Theta$ can be determined by the inverse problem method according to the required thrust value.

In this article, we adopt the blades provided by the partner airlines which shape can be shown in Figure 15. Its meshes in the R-coordinate system are shown in Figure 1, and 51 2D-manifolds are used in this case. The velocity at inlet is 80 m/s and the rotating angular velocity is 80 rad/s. Figure 11 – Figure 12 show the pressure distribution on the blade surface. Because of the rotation, positive pressure surface and negative pressure surface have different pressure distribution. The pressure distribution of the positive pressure surface is higher than that of the negative pressure surface. Because the blade is cocked up at the bottom (see Figure 15), the pressure distribution near the bottom of the blade is obviously different from other places. Through formula transformation, we present the 3D model of blade pressure and velocity in Figure 13 – Figure 14.

Meanwhile, we show the velocity distribution at outlet in Figure 16. The velocity on the outside is significantly higher than that on the inside, and the velocity in the middle of the channel is also higher than that near the blade, which is consistent with common sense. From Figure 16, we can find that the maximum velocity is located at the edge of the blade, which can reach 130 m/s.

5.3.   Discussion

In this work, we conduct numerical simulations for 3D flow in the flow channel of the aeroengine turbine fan based on the dimension splitting method. However, the simulation of the flow states between the blades is not the ultimate goal. Our ultimate goal is to use the simulation results to design or optimize the blades. Meanwhile, the parameterized blade surface provides convenience for blade design and optimization. Thus, it will be the focus of follow-up research to optimize the blade by combining the inverse problem. In addition, compressible flow and boundary layer phenomenons are also the focus of the follow-up research.

Acknowledgments.

The work of G. Ju was partially supported by the NSF of China (No. 11731006), the Shenzhen Sci-Tech Fund (No. JCYJ20170818153840322), and Guangdong Provincial Key Laboratory of Computational Science and Material Design No. 2019B030301001. The work of R. Chen was supported by the NSF of China (No. 61531166003) and Shenzhen Sci-Tech fund (No. JSGG20170824154458183 and ZDSYS201703031711426). The work of J. Li was partially supported by the NSF of China (No. 11971221) and Shenzhen Sci-Tech fund (No. JCYJ20190809150413261). The work of K. Li was supported by the NSF of China under the grant No. 10971165 and 10771167.