Optimal control of thermal and mechanical loads in activation processes of mechanical components

Abstract

This paper develops a mathematical framework that aims to control the temperature and rotational speed in the activation process of a gas turbine in an optimal way. These controls influence the deformation and the stress in the component due to centripetal loads and transient thermal stress. An optimal control problem is formulated as the minimization of maximal von Mises stress over a given time and over the whole component. To find a solution for this, we need to solve the linear thermoelasticity and the heat equations using the finite element method. The results for the optimal control as functions of the rotation speed and external gas temperature over time are computed by sequential quadratic programming, where gradients are computed using finite differences. The overall outcome reveals a significant reduction of approximately 10%, from $830\frac{\mathrm{N}}{{\mathrm{mm}}^2}$ to $750\frac{\mathrm{N}}{{\mathrm{mm}}^2}$, in von Mises stress by controlling two parameters, along with the temporal separation of physical control phenomena.

Keywords:

• Activation of mechanical components
• optimal control
• finite element method
• sequential quadratic programming

Mathematics Subject Classifications:

• 49M40
• 74S05
• 80M10

1. Introduction

When technical facilities are activated or deactivated, the built-in components undergo load cycles that cause fatigue and, ultimately, failure. The safe operation of such facilities therefore requires service concepts to counteract the material fatigue by repair and replacement, but also concepts for the gentle operation of components that minimize the loads occurring. While this is general wisdom, so far little mathematical or computational research effort has been invested to minimize negative effects of activation cycles by controlling them actively and in an optimal way. In this article, we present a first study towards the optimal control of activation and deactivation processes.

A technical device, where fatigue is highly relevant, is the gas turbine, which is mostly used in jet engines for aviation, but also for energy production. In both fields, activation and deactivation happen frequently, e.g. during start and landing, or when gas turbines are used in the energy grid to balance volatile energy influx by renewable sources of energy production. For gas turbines, the loss of single turbine disks is contained in the turbine casing, but a rotor burst leads to a disintegration of the entire device with potentially catastrophic consequences. Therefore, the mechanical integrity of gas turbine rotor disks is crucial. At the same time, peak stresses in a gas turbine rotor do not necessarily occur under full load, but during the activation process itself. This is so, as the rotor usually is the heaviest part of the turbine with the largest thermal inertia. So, during activation, the exterior of the rotor is in contact with the hot working fluid and undergoes thermal expansion, whereas the interior part is still cool, which leads to transient thermal stress. At the same time, rotation speed and burner temperature are controllable for gas turbines, where constraints have to be respected for the duration of the activation process. This complex behaviour makes gas turbines predestined for our study and consequently, we work with an example motivated by the activation of a stationary gas turbine.

The coupling of heat transfer and elastic deformation is a widely used phenomenon in many engineering applications where body parts are in varying temperature fields and have to withstand to mechanical loads. In Ref. [1] an optimal control of the heating of a thermoelastic plate by internal heat sources is given. In Ref. [2] a thermoelastic sensitivity analysis has been made. From this work, the underlying equations were adapted. In Ref. [3], an online optimal control scheme of inlet steam temperature during start up of steam turbines has been made. Contrary to this work, here we apply an optimal control with two control variables, the temperature and the rotation. In Ref. [4], an optimization of a rotor with respect to the turbine load was carried out. Also turbine topology optimization was conducted in Refs [5–7] regarding to thermal and mechanical loads. In Refs [8–10], a detailed analysis for the heat transport rate was carried out. Here, we use a simplified boundary condition for the working fluid with a fixed convective heat transfer coefficient.

In this paper, we will compute a solution of an optimal control problem for an activation process of a turbine regarding to thermal and mechanical loads. The mechanical as well as the thermal load will be the control variables over a given time and the goal is to reduce the maximum von Mises stress over the whole time. In general optimal control problem can be solved in two ways: First optimize then discretize (FOTD) or first discretize then optimize (FDTO) [11]. We will use the latter one. This means we will first discretize the volume with a 3D mesh, the time with implicit Euler and then solve the resulting finite-dimensional optimization problem.

The paper is organized as follows: Section 2 provides a detailed technical background and introduces the model formulation. Section 3 delves into the optimal control problem, discussing objectives, and constraints. In Section 4, a Finite Element Method (FEM) Simulation model is presented, comprising a weak formulation and an overview of preprocessing and implementation techniques. Section 5 explores the solution of the optimal control problem using Sequential Least Squares Programming (SLSQP) algorithm. Finally, in Section 6, a summary of the paper's findings is provided, along with an outlook on potential future research.

2. Technical background and model formulation

Gas turbines play an essential role as back up capacities in energy grids with volatile influx of renewable energy from solar panels or wind turbines. This makes it necessary to start gas turbines quickly in order to meet the consumer's demands.

During the starting procedure, the gas turbine burner heats up the flowpath, while the gas turbine starts turning. While the rotation leads to centripetal mechanical loads, the non-uniform distribution of heat during the activation process results in thermal loads. Both mechanisms jointly cause mechanical fatigue which effectively limits the life of the gas turbine. For safe operation it is therefore desirable, to minimize peak stress in important components via a controlled activation procedure that increases the angular velocity and the temperature in an optimized manner.

Thermal inertia during heating leads to the aforementioned non-uniform temperature distribution and thereby, potentially, to transient peak stress. This mechanism is more severe for massive parts, like the gas turbine rotor, and is less important for parts with little mass, as gas turbine blades and vanes. At the same time, the gas turbine rotor disks are the most critical parts of the turbine, as a rotor burst is not contained by the casing. We therefore choose a simplified model of a rotor disk as the most relevant application for our optimized activation half-cycle.

Let us now present the details of our mathematical model. Throughout the paper, matrices are represented by bold uppercase letters, vectors by bold lowercase letters, and scalars by lowercase normal letters. Let $\mathrm{\Omega }\subseteq {\mathbb{R}}^3$ be the portion of the three-dimensional space that is filled by the component, cf. Figure 1 and let $\tau =[0,t_f]$ be the underlying time sequence with a final time $t_f$. The region Ω is tied to the Cartesian coordinate system and the front centre of the cylinder in Figure 1(b) is coordinate origin $(0,0,0{)}^T$. When an external load $\mathit{f}:\tau \times \mathrm{\Omega }\to {\mathbb{R}}^3$ is applied, the elastic material filling Ω undergoes a deformation that can be modelled by a vector field $\mathit{u}=(u_x,u_y,u_z{)}^T:\tau \times \mathrm{\Omega }\to {\mathbb{R}}^3$ with $u_x$ denoting the deformation in x-direction, and $u_y$ and $u_z$ representing the deformations in y- and z-direction, respectively. Under the assumption of linear elasticity and within the context of linear small deformation theory, $\mathit{u}$ fulfils the partial differential equation of thermoelasticity (1) $\mathrm{\nabla }\cdot \mathbf{\Sigma }\left(\mathit{u}\right)+\mathit{f}=\mathbf{0}\mathrm{in}\mathrm{\Omega }$(1) where (2) $\mathbf{\Sigma }\left(\mathit{u}\right)=2\mu \mathcal{E}\left(\mathit{u}\right)+[\mathrm{\lambda Tr}(\mathcal{E}\left(\mathit{u}\right))-\alpha (3\lambda +2\mu \left)\right(T-T_0\left)\right]\mathit{I}\mathrm{in}\mathrm{\Omega }$(2) is the stress tensor with $\lambda ,\mu >0$ the Lamé constants, $\mathit{I}$ the $3\times 3$ identity matrix and $T:\tau \times \mathrm{\Omega }\to \mathbb{R}$ some given field of temperatures. $T_0$ is a reference temperature that corresponds to the temperature state of the non-activated component. $\mathrm{Tr}$ denotes the trace operator. Thermal stresses are driven by α, thermal expansion coefficient, and the difference between $T\left(x\right)$, the temperature at $x\in \mathrm{\Omega }$ and $T_0$. The elastic strain tensor is given by (3) $\mathcal{E}\left(\mathit{u}\right)=\frac{1}{2}(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^{\mathrm{\top }})$(3) and the loading due to centripetal acceleration due to rotation around the z – axis of the axial turbo machine (4) $\mathit{f}=\rho {\omega }^2{\mathit{x}}_{\mathit{p}},{\mathit{x}}_{\mathit{p}}=(x,y,0{)}^{\mathrm{\top }}.$(4) The density is denoted as ρ, ω is the rotation and ${\mathit{x}}_{\mathit{p}}$ are the projected radial components of the coordinate vector $(x,y,z{)}^{\mathrm{\top }}\in \mathrm{\Omega }$ where these rotation force apply to.

The temperature distribution $T:\tau \times \mathrm{\Omega }\to \mathbb{R}$ is obtained as the solution of the time-dependent heat equation (5) $\begin{array}{rl}\rho c_p\frac{\mathrm{\partial }T}{\mathrm{\partial }t}-k\mathrm{\Delta T}& =0\mathrm{in}\mathrm{\Omega };\\ T(t=0)& =T_0\mathrm{in}\mathrm{\Omega };\end{array}$(5) where $\mathrm{\Delta }=\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}$ is the Laplace operator, $c_p$ the specific heat capacity and k is the thermal conductivity. The boundary conditions at the surface where heat exchange with the heated working fluid is possible are given by Robin or convective boundary conditions (6) $-k\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}=h(T-T_{\mathrm{e}})\mathrm{in}(\mathrm{\partial \Omega }{)}_R$(6) where $\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{n}}$ stands for the directional derivative in the direction of the outward normal vector $\mathit{n}$. The heat transfer coefficient h>0, is assumed to be fixed¹. The notation $(\mathrm{\partial \Omega }{)}_R$ denotes a Robin boundary and the surfaces, where this type of boundary conditions is applied, are depicted in Figure 2. Note that the speed of rotation $\omega =\omega \left(t\right)$ and the external temperature $T_{\mathrm{e}}\left(t\right)$ are both variables that can be controlled during the activation process. This, together with the time-dependent heat Equation (5(5) $\begin{array}{rl}\rho c_p\frac{\mathrm{\partial }T}{\mathrm{\partial }t}-k\mathrm{\Delta T}& =0\mathrm{in}\mathrm{\Omega };\\ T(t=0)& =T_0\mathrm{in}\mathrm{\Omega };\end{array}$(5) ) induces time-dependent solutions $\mathit{u}\left(t\right)$ of (1(1) $\mathrm{\nabla }\cdot \mathbf{\Sigma }\left(\mathit{u}\right)+\mathit{f}=\mathbf{0}\mathrm{in}\mathrm{\Omega }$(1) ).

Figure 1. Full (left, (a)) and reduced (right, (b)) 3D model. (a) Full 3D model of a simplified turbine. Outer diameter is about 1.4 m and represents a larger stage of a turbine. This stage is especially interesting for thermal and rotational loading. (b) Due to axis symmetry reduced model with only one turbine blade. The maximal load is expected to be in the transition from turbine blade to rotor. This is the calculation volume Ω.

Figure 2. Boundary surface $(\mathrm{\partial \Omega }{)}_R$ for external heating. The ‘R’ subscript abbreviates for Robin boundary.

To conclude this section, let us complement the boundary conditions for a computational domain that makes use of the component's rotation symmetry. In order to avoid the redundant numerical computation of temperature and stress in identical segments, we simplify the blade count to one and introduce plane symmetry boundary condition for the cylinder part on the 'internal' boundary sections of one specific segment, see Figure 1(b), i.e. we set (7) $\begin{array}{rl}& u_y(t,0,y,z)=0,u_x(t,x,0,z)=0,\end{array}$(7) (8) $\begin{array}{rl}& \frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,0,y,z)=0,\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,0,z)=0,\end{array}$(8) (9) $\begin{array}{rl}& \frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,y,0)=0,\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,y,z_l)=0.\end{array}$(9) The boundary conditions (7(7) $\begin{array}{rl}& u_y(t,0,y,z)=0,u_x(t,x,0,z)=0,\end{array}$(7) ) for $\mathit{u}$ and (8(8) $\begin{array}{rl}& \frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,0,y,z)=0,\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,0,z)=0,\end{array}$(8) ), (9(9) $\begin{array}{rl}& \frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,y,0)=0,\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,y,z_l)=0.\end{array}$(9) ) for T are symmetric and null Neumann boundary conditions, where $z_l$ is the length of the part in z-direction. This means the deformation perpendicular to the x and y axis is zero as well as the part is isolated at any boundary except for the Robin boundary (6(6) $-k\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}=h(T-T_{\mathrm{e}})\mathrm{in}(\mathrm{\partial \Omega }{)}_R$(6) ). The first boundary conditions in (7(7) $\begin{array}{rl}& u_y(t,0,y,z)=0,u_x(t,x,0,z)=0,\end{array}$(7) ) and in (8(8) $\begin{array}{rl}& \frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,0,y,z)=0,\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,0,z)=0,\end{array}$(8) ) are applied along the bottom surface in Figure 1(b). The second boundary conditions in (7(7) $\begin{array}{rl}& u_y(t,0,y,z)=0,u_x(t,x,0,z)=0,\end{array}$(7) ) and in (8(8) $\begin{array}{rl}& \frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,0,y,z)=0,\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,0,z)=0,\end{array}$(8) ) are applied along the right-hand side surface of the cut cylinder in Figure 1(b). The boundary conditions in (9(9) $\begin{array}{rl}& \frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,y,0)=0,\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}(t,x,y,z_l)=0.\end{array}$(9) ) are applied along the front and the back surface of the cylinder, respectively. We also do not apply any forces on the surface $\mathrm{\partial }{\mathrm{\Omega }}_D$ in Figure 2. In a later, refined modelling approach, pressure boundary conditions from the static pressure of the working fluid should be applied as well. However, this involves complex 3D fluid dynamics simulation and is beyond the scope of the present article.

3. Optimal control problem

The durability of turbine (as well as other engineering system) components is known to be enhanced if they are exposed to minimal stresses inside. Stress peaks can lead to unwanted deformations, porosities or even cracks in the parts. The strain of turbine components is particularly high during the starting process.

We will therefore formulate an optimal control problem for the turbine activation. Let us consider the rotation $\omega :\tau \to \mathbb{R}$ in (4(4) $\mathit{f}=\rho {\omega }^2{\mathit{x}}_{\mathit{p}},{\mathit{x}}_{\mathit{p}}=(x,y,0{)}^{\mathrm{\top }}.$(4) ) and the external temperature $T_e:\tau \to \mathbb{R}$ in (6(6) $-k\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}=h(T-T_{\mathrm{e}})\mathrm{in}(\mathrm{\partial \Omega }{)}_R$(6) ) as control functions in the fixed and given time domain $\tau =[0,t_f]$ of the starting process. The goal of the optimal control problem is to minimize maximal stresses considered over time and space, i.e. we minimize the cost functional (10) $J(T_{\mathrm{e}},\omega ):=\underset{t}{max}\left(\underset{x,y,z}{max}\right({\sigma }_v\left)\right).$(10) Here, ${\sigma }_v={\sigma }_v(t,x,y,z)$ is the von Mises stress given by (11) ${\sigma }_v=\sqrt{\frac{3}{2}{\mathbf{\Sigma }}^{\mathbf{dev}}:{\mathbf{\Sigma }}^{\mathbf{dev}}}$(11) where (12) ${\mathbf{\Sigma }}^{\mathbf{dev}}=\mathbf{\Sigma }-\frac{1}{3}\mathrm{Tr}\left(\mathbf{\Sigma }\right)\mathit{I}$(12) is the stress deviator tensor with $\mathbf{\Sigma }$ as defined in (2(2) $\mathbf{\Sigma }\left(\mathit{u}\right)=2\mu \mathcal{E}\left(\mathit{u}\right)+[\mathrm{\lambda Tr}(\mathcal{E}\left(\mathit{u}\right))-\alpha (3\lambda +2\mu \left)\right(T-T_0\left)\right]\mathit{I}\mathrm{in}\mathrm{\Omega }$(2) ) and the matrix scalar product notation $\mathit{A}:\mathit{B}=\mathrm{Tr}\left({\mathit{A}}^{\mathrm{\top }}\mathit{B}\right)=\sum _{i,j=1}^n{a}_{\mathrm{ij}}{b}_{\mathrm{ij}}$for two $n\times n$-matrices $\mathit{A}$, $\mathit{B}$ with entries $a_{\mathrm{ij}}$, $b_{\mathrm{ij}}$, $1\le i,j\le n$, has been used. From a technical point of view, we assume the control functions to be continuous. As discussed in the next paragraph, this will also be mathematically reasonable and furthermore, ω is assumed to be differentiable.

The starting process is determined by the three conditions that the external temperature $T_e$ has reached at least some desired temperature $T_{e,f}$, that the rotation ω has reached a desired rotation ${\omega }_f$ at the final time $t_f$ and that a the maximum temperature reached at least some desired condition. We hence postulate (13) $T_e\left(t_f\right)\ge T_{e,f},\omega \left(t_f\right)={\omega }_f,\underset{x,y,z}{max}T\left(t_f\right)\ge T_f.$(13) We can see that continuity of the control functions is essential in this formulation to prohibit solutions which jump from zero to the desired state in the final time $t_f$. Furthermore, we are interested in solutions with bounded rotation acceleration to prevent peaks in the rotation velocity, i.e. we set (14) $\frac{\mathrm{\partial }\omega }{\mathrm{\partial }t}\le {\omega }_{\lim },t\in \tau$(14) as pointwise constraint for the derivative of the control ω. This assumption is also natural, as the rotational inertia forces will prevent an arbitrary acceleration of the angular velocity. In addition the control variable underlies some bounds (15) $\begin{array}{rl}& {\omega }_{min}<\omega <{\omega }_{max},\end{array}$(15) (16) $\begin{array}{rl}& T_{e,min}<T_e<T_{e,max}.\end{array}$(16) These bounds will start at ${\omega }_{min}=0$ and $T_{e,min}=0$ and will go to a maximum of ${\omega }_{max}=60$ hz and $T_{e,max}={1000}^{\circ }$C. Finally, the dynamics of our model with boundary and initial conditions, as given in Section 2, are constraints of the problem.

To summarize, the optimal control problem reads: (17) $\begin{array}{rl}& \text{Minimize}\left(10\right)\\ & \mathrm{subject}\mathrm{to}\left(1\right),\left(5\right),\left(7\right),\left(8\right),\left(9\right),\left(13\right),\left(14\right)\left(15\right)\left(16\right).\end{array}$(17)

4. FEM simulation model

FEM is a method for solving partial differential equations and it is very widely used in the engineering field to solve structural mechanics problems. We are looking on a structural mechanics problem coupled with the heat equation. One of the main strengths is the flexibility to use complex computational domains, so we can calculate 3d models easily. This section captures the used methods in this work. Firstly, we will show the weak formulation of the governing equation which we need for the implementation afterwards. Then we show our 3D model and the FEM Model where we used FEniCS [12, 13], an open source FEM solver. We also provide an overview about the workflow in the appendix because it is not trivial when getting started with complex 3d domains in FEM.

4.1. Weak formulation

To use FEM, we have to formulate the weak form of the governing equations. We adapted the formulation from Refs [2, 14] and added a transient part. Applying the weak formulation on the heat Equation (5(5) $\begin{array}{rl}\rho c_p\frac{\mathrm{\partial }T}{\mathrm{\partial }t}-k\mathrm{\Delta T}& =0\mathrm{in}\mathrm{\Omega };\\ T(t=0)& =T_0\mathrm{in}\mathrm{\Omega };\end{array}$(5) ) yields (18) ${\int }_{\mathrm{\Omega }}\rho c_p\frac{\mathrm{\partial }T}{\mathrm{\partial }t}\cdot v_T\mathrm{d}x-{\int }_{\mathrm{\Omega }}k\mathrm{\Delta T}\cdot v_T\mathrm{d}x=0$(18) where $v_T\in H^1\left(\mathrm{\Omega }\right)$ is a test function for the temperature and $\mathrm{d}x$ denotes the integration over the whole domain Ω. Hereby, $H^1$ is the Sobolev space of functions whose derivative is $L^2$ integrable; see, e.g. Ref. [11] for details. Using partial integration for the second term results in (19) ${\int }_{\mathrm{\Omega }}\rho c_p\frac{\mathrm{\partial }T}{\mathrm{\partial }t}\cdot v_T\mathrm{d}x+{\int }_{\mathrm{\Omega }}k\mathrm{\nabla T}\cdot \mathrm{\nabla }{v}_T\mathrm{d}x-k{\int }_{\mathrm{\partial \Omega }}(k\mathrm{\nabla T}\cdot \mathit{n})v_T\mathrm{d}A=0.$(19) Here dA denoted the integration over the a boundary $\mathrm{\partial \Omega }$. With help of Equation (6(6) $-k\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}=h(T-T_{\mathrm{e}})\mathrm{in}(\mathrm{\partial \Omega }{)}_R$(6) ), the last term can be written as (20) $-{\int }_{\mathrm{\partial \Omega }}(k\mathrm{\nabla T}\cdot \mathit{n})v_T\mathrm{d}A=-{\int }_{\mathrm{\partial \Omega }}k\frac{\mathrm{\partial }T}{\mathrm{\partial }\mathit{n}}{v}_T\mathrm{d}A=-{\int }_{\mathrm{\partial \Omega }}h(T-T_{\mathrm{e}})v_T\mathrm{d}A$(20) The time derivatives are now replaced by a backward Euler scheme, so that the previous weak form at the time increment n + 1 is now: (21) ${\int }_{\mathrm{\Omega }}\rho c_p{\displaystyle \frac{T_{n+1}-T_n}{\mathrm{\Delta }t}}\cdot v_T\mathrm{d}x+{\int }_{\mathrm{\Omega }}k\mathrm{\nabla }{T}_{n+1}\cdot \mathrm{\nabla }{v}_T\mathrm{d}x=-{\int }_{\mathrm{\partial \Omega }}h(T_{n+1}-T_{\mathrm{e}}))v_T\mathrm{d}A$(21) The same has to be done for the thermoelasticity Equation (1(1) $\mathrm{\nabla }\cdot \mathbf{\Sigma }\left(\mathit{u}\right)+\mathit{f}=\mathbf{0}\mathrm{in}\mathrm{\Omega }$(1) ): (22) ${\int }_{\mathrm{\Omega }}\mathbf{\Sigma }\left(u\right):\mathcal{E}\left(v_M\right)\mathrm{d}x={\int }_{\mathrm{\Omega }}f\cdot v_M\mathrm{d}x+{\int }_{\mathrm{\partial \Omega }}{v}_M\cdot \mathbf{\Sigma }\left(u\right)\cdot n\mathrm{d}A$(22) where $v_M\in H^1\left(\mathrm{\Omega }\right)$ is the test function for mechanical deformation and f the linear functional corresponding to the force of external forces. The integral part with $v_M\mathbf{\Sigma }\left(u\right)\cdot n$ is often referred as a traction force which is applied on as a Neumann boundary condition on a surface. In our case, this traction force will be zero, as we neglect external gas pressure. This is also everything we need for later implementation. We can use left hand side (lhs) and right hand side (rhs) function from FEniCS to extract the linear and bilinear forms. However, for completion purpose, we will conclude the linear and bilinear forms (23) $\begin{array}{rl}{a}_T(T,v_T)& ={\int }_{\mathrm{\Omega }}\rho c_p{\displaystyle \frac{T_{n+1}-T_n}{\mathrm{\Delta }t}}\cdot v_T\mathrm{d}x+{\int }_{\mathrm{\Omega }}k\mathrm{\nabla }{T}_{n+1}\cdot \mathrm{\nabla }{v}_T\mathrm{d}x+{\int }_{\mathrm{\partial \Omega }}h{T}_{n+1}{v}_T\mathrm{d}A,\end{array}$(23) (24) $\begin{array}{rl}{a}_M(u,v_m)& ={\int }_{\mathrm{\Omega }}\mathbf{\Sigma }\left(u\right):\mathcal{E}\left(v_M\right)\mathrm{d}x,\end{array}$(24) (25) $\begin{array}{rl}{L}_T\left(T\right)& ={\int }_{\mathrm{\partial \Omega }}h{T}_{\mathrm{e}}{v}_T\mathrm{d}A,\end{array}$(25) (26) $\begin{array}{rl}{L}_M\left(u\right)& ={\int }_{\mathrm{\Omega }}\mathit{f}{v}_M\mathrm{d}x.\end{array}$(26) The last thing we need to consider for the later implementation is the type of the finite element. We will use a second-order Lagrange element for the displacement field and a first order Lagrange for the temperature. For details about finite element definition, we refer to Refs [2, 12, 15].

4.2. Preprocessing and implementation

We created our 3D data with the commercial engineering tool Autodesk Inventor. The full model can be seen in Figure 1(a). Although the part is only an academic example, it is important that we place roundings at the otherwise sharp corners to get reasonably realistic results. To save memory and time we reduce the full model to a smaller one. Here we are using symmetry condition in x and y axes as well as using only one turbine blade as seen in Figure 1(b). To use the Cartesian coordinate system will simplify these boundary conditions afterwards. To mesh the 3D geometry, we must export the data to ‘stl’ format, a common file format. Then we have to discretize the 3d model. We used ‘Gmsh’ to do this. There we marked the whole volume model (tetrahedrons) as the computation domain Ω and necessary boundary surfaces (triangles) as the boundary conditions $\mathrm{\partial \Omega }$. After the creation of the mesh we translated this to a mesh format we can use with FEniCS. We used ‘meshio’ as a python module to convert the mesh into readable data structure for FEniCS. The overview in table formatting can be found in appendix.

We implemented the thermoelasticity problem in FEniCS. This can be done by defining the problem as the weak formulation. We defined two function spaces for the two test functions from our weak form. We use two different element types in Ω: the displacement field is approximated by Lagrange (FEniCS:Continous Galerkin) elements with order 2 and the temperature field is approximated by Lagrange elements of order one. In order to minimize computation time, we opted for the smallest possible element order that still yields reasonable stress results. The stress field has to be projected afterwards onto another space to relocate the stress to the nodes. Here we take a discontinuous Galerkin (DG) space of order 1 because the stress depends on the gradient of the displacement field which is of order 2. So we get at best a piecewise discontinuous results.

The time is discretized in 20 equidistant steps. We first started with a logarithmic timestepping for better resolution in the first moments of the process. But it turns out that the time step should have some minimal value otherwise the solution of the heat equation shows some weird behaviour. The values get negative despite the initial condition are non-negative as well as the boundary conditions are non-negative with positive time stepping with a backward Euler. This behaviour is also seen in Ref. [16]. To preserve positivity, we increased the time step.

To conclude our FEM, the input for this problem is the control of rotational speed and external temperature over 20 timesteps, which results in a 40-dimensional input. The model will give transient results of the stress and temperature over a time of 1800 s. We are especially interested in the maximum von Mises stress in each timestep.

4.3. Definition of a starting guess control

The solution obtained from the finite element model corresponds to the input variables of external temperature and rotational speed. To establish a reasonable starting guess control, we adopt a linearly increasing rotation and temperature profile. This initial guess is illustrated in Figure 3. The outcome of a single finite element analysis entails determining the maximum von Mises stress (green) across the domain and the maximum heat (red) on the surface as shown in Figure 3 with the two control variables heat control (blue) and rotational speed (orange). This means we have a function (FEA) with 40 control variables, 20 external temperatures $T_e$ as indicated in blue and 20 rotational speeds ω as indicated in orange over the time domain. The response of the function is the maximum von Mises stress ${\sigma }_v$ and the maximum heat T over the whole domain. The maximum von Mises stress occurs at the final time step and reaches a value exceeding $830\frac{N}{m{m}^2}$. All imposed constraints and bounds are satisfied within this solution. This means, the maximum heat observed on the surface exceeds 400${}^{\circ }$C, with the rotation rate reaching 60 Hz ($2\pi 60$ in plot) and the heat control reaching 750${}^{\circ }$C at the conclusion of the analysis.

Figure 3. Initial guess of linear increasing speed and heat controls and solution of resulting maximum von Mises stress and maximum heat in the part over the whole simulation time. The linearly increasing control leads to the fulfilment of the constraints and seems to be a reasonable approach, since a slow heating naturally leads to a lower stress. But this initial start guess turns out to be non-effective for gradient-based optimization.

Nevertheless, in the subsequent chapter, we will demonstrate that this initial estimation is inadequate. We will present another initial guess, as indicated in Figure 4, that may not appear reasonable as a solid at first, but serves as a much better starting point for gradient-based optimization.

Figure 4. The maximum stress over the whole domain, as seen in green line, is decreasing from Iteration 1 to Iteration 17. In Iteration 1, another initial guess for the optimization procedure can be seen. A jump in heat control followed constant heat control and a speed control as late as possible will result in a better optimal control.

4.4. Statistics and performance of the FEM model

In this finite element analysis (FEA), we used a mesh consisting of 1879 nodes, 1368 triangles, and 7339 tetrahedrons. The transient time step of the simulation was set to 90 s, and the simulation covered a total of 1800 seconds or 30 min with 20 time steps. When run on a cluster with hardware specs of 2x64 Core AMD EPYC Milan, 1024GB of RAM the same simulation took only about 1 min. This demonstrates the power of utilizing a cluster for running FE simulations, as it can greatly reduce the computational time required and is even more relevant in the optimization.

5. Solving the optimal control problem

Solving optimal control problems involves finding the best set of control inputs to achieve a desired outcome within a system. One way to solve these problems is through the use of numerical optimization techniques, such as the Sequential Quadratic Programming (SQP) algorithm. In this case study, we will use an enhancement of the SQP algorithm, the Sequential Least Squares Quadratic Programming (SLSQP) implemented in the Scipy library in Python which was developed in Ref. [17]. We will take a look of the results of the optimization to understand how the control inputs affect the stress distribution within the part and understand the optimal solution that minimizes the overall stress.

In this chapter, we will explore the use of the SLSQP method for solving an optimal control problem. We will be using the Scipy library in Python, which is well suited for our problem. The implementation of Scipy is also able to deliver good results without a gradient which we will use and explain why we can not use automatic differentiation tools for the gradient.

5.1. Optimal control with SLSQP method

In this particular case, we will examine an optimal control problem that involves two control inputs: external heating and rotation. With a time step discretization of 20 and 2 input variables, we encounter a 40-dimensional optimization problem. The control variables are subject to specific limits: external heating must fall within the range of 0${}^{\circ }$C to 1000${}^{\circ }$C, while rotation must be between 0 and 60 Hz. However, it should be noted that the control temperature corresponds to the temperature of the working fluid. Furthermore, in addition to these bounds, the optimization algorithm includes additional constraints: the external heating should have a value of 750${}^{\circ }$C at the final time step, and the rotation should be 60 Hz at the final time step. Additionally, we have incorporated a constraint to indicate the maximum temperature at the surface, which must reach 400${}^{\circ }$C by the end of the simulation.

The objective of this optimization is to minimize the maximum stress over all time steps and across the entire domain. Mathematically, this can be represented as shown in Equation (17(17) $\begin{array}{rl}& \text{Minimize}\left(10\right)\\ & \mathrm{subject}\mathrm{to}\left(1\right),\left(5\right),\left(7\right),\left(8\right),\left(9\right),\left(13\right),\left(14\right)\left(15\right)\left(16\right).\end{array}$(17) ). The optimization algorithm iteratively adjusts the control inputs to minimize this objective, resulting in a decrease in the highest stress point, as illustrated by the green line in Figure 4.

Given that the SLSQP method is a gradient-based optimizer, providing a good initial guess for the control inputs is crucial. As mentioned earlier, we attempted using linearly increasing controls as the initial guess, but it turned out to be ineffective. A more effective initial solution involves maximum heat and minimal rotation, which helps in finding an optimal solution more quickly.

Various libraries offer automatic gradient calculation features. However, utilizing these automatic gradient procedures, such as Jax [18] or Pyomo [19], require the calculation to be implemented within the codebase of the respective libraries. Since we are using a C++-based code within FEniCS, compatibility is limited. Another option to obtain the gradient would be Dolfin-adjoint [20], similar to Tensorflow [21], which constructs a differentiable graph. However, this approach presents challenges as the objective function needs to be an integral over the entire domain, whereas our study specifically focuses on the maximum von Mises stress. Therefore, an integral over the entire domain is not useful in this research. So the gradient calculation is performed using finite differences in Scipy. We utilized the 2-point method, which is the default procedure when no gradient is provided to the optimization algorithm. The default step size is approximately $\sqrt{\mathrm{machine}\mathrm{precision}}\approx$ 1.49e-08, requiring one finite element method (FEM) call per optimization variable. Consequently, 40 FEM calls are made to obtain one gradient.

The SLSQP algorithm employs a BFGS approximation of the Hessian, as described by Powell [22], for constrained optimization problems. The optimal step size of the algorithm, analogous to Newton's method for nonlinear optimization problems, is $\alpha =1$ near the local optimum. However, for vectors far from the optimum, the step size needs to be modified according to the penalty function described by Han [23], which includes a non-differentiable merit function. This function is substituted by a differentiable augmented Lagrange function presented by Hock and Schittkowski [24].

In total, 179 iterations and over 7000 function evaluations (FEM simulation calls) were required to identify an optimal control with a function tolerance of ${10}^{-8}$. The optimal control is visualized in Figure 5. The SQP procedure invokes the FEM solver over 7000 times, resulting in a cumulative computation time of 7000 min or 120 h. Therefore, the need for a fast FEM simulation is clearly evident. When the SQP method necessitates multiple restarts with different initial guesses, this time requirement becomes even more critical.

Figure 5. Optimal control results of Scipy's SLQSP in terms of speed and heat control. Heat control starts at a high value and then decreases as the stress increases. The stress increases and reaches a maximum, and the control tends to bring the stress as close to the maximum as possible since the heat in the part can then increase as fast as possible. In the penultimate time step, both the heat control and the stress decrease to meet the boundary conditions in the final time steps while not exceeding a maximum stress threshold.

5.2. Physical interpretation

In the present study, we observe a tendency in the optimal control strategy wherein the heating process is prioritized initially, followed by a delayed initiation of rotation. This sequential approach appears to be reasonable due to the nature of the influence exerted by rotational speed on von Mises stress, which is perceived as an offset. Conversely, the influence of temperature on stress diminishes once the component is fully heated. The optimal control tries for control to get a constant von Mises stress as the objective function is only influenced by the maximum over time. This results in an optimal use of time if the stress is nearly constant.

The position of the maximum stress is nearly at the same place. This is also one major advantage for gradient optimization when minimizing a maximum because the the objective function will not take big jumps as response and remain continuous, which is not directly given by this objective function. The stress positions are illustrated in Figure 6. Moreover, our proposed solution leads to a reduction in maximum von Mises stress and an increase in internal heat. Specifically, the von Mises stress is lowered from $830\frac{\mathrm{N}}{{\mathrm{mm}}^2}$ (see Figure 3, last time step) to $750\frac{\mathrm{N}}{{\mathrm{mm}}^2}$ (see Figure 5), resulting in a decrease of approximately 10%. Additionally, we establish that minimizing rotation is desirable to mitigate the occurrence of unavoidable von Mises stress. The break in von Mises stress and heat control in the penultimate time step is either explainable due to the evolution of the optimization (see Figure 4) where a ‘overshoot’ can happen and never be reversed. This can occur due to the objective function which only states the maximum von Mises stress and therefore never change non-maxima. Another explanation is, that the constraints of this optimization problem at the last time step need to be satisfied. The algorithm needs to reduce heat control in penultimate time step in order to satisfy constraints and not increase the maximum von Mises stress in the last time step.

Figure 6. Position of the maximum von Mises stress in timestep 20. Only a small partition of the blade is under high mechanical stress. The highest stress is marked by a pink dot and the node number in yellow.

At first we did not bound the gradient of the rotational speed which resulted in a jump in the speed at the last timestep, which is technically not possible.

6. Summary and outlook

The paper addresses the startup of a turbine through a transient finite element model incorporating thermoelasticity, with the heat equation, and rotational stress over 20 timesteps with two controls that can be varied over time: the rotational speed and the external heat. The aim is to minimize peaks in the von Mises stress through a control strategy of first heating and then rotating it as late as possible. This results in a convenient use of material and time in the sense that as much heat as possible is saved in the material before the rotational stress is applied. The slope of the rotational speed was limited by a constraint and a physically reasonable control was found. The von Mises stress peaks have been decreased by 10 % which can lead to longer lifetime or safer components overall. Up to the knowledge of the authors, for the first time, optimal control of an activation process with two control variables was successfully realized in a three-dimensional domain, leading to a 10% reduction in peak von Mises stress.

There are two types of limitations to this analysis. The first one concerns the abstracted physical behaviour: we assumed the convective heat coefficient h to be fixed and the inertia force is also neglected in this study which is not true in reality. This can be implemented in future work. The second limitation is the need for faster simulation models when extending the optimization from one blade to a full turbine model. Creating a surrogate model can drastically decrease the computational time for one simulation. The full FEM model, or single FEM outputs, or the objective can be trained in such an approach. The surrogate model must be trained globally as well as in the optimal region. One can use latin hypercube sampling for global input parameters and an iterative approach with Bayesian optimization [25] to get new samples in the optimal region. Finally, also new methods such as physics-informed neural networks [26] could be used.

Disclosure statement

The authors declare that no conflict of interest regarding financial or personal relationships exists that could have influenced this work.

Data availability statement

A repository containing the optimal control script including the FEM solver and a mesh preprocessing can be found at: https://github.com/nifri004/OC_Thermoelasticity

Additional information

Funding

This paper was produced as part of the joint research project ML-Real-Time, which is funded by the German Federal Ministry of Education and Research (BMBF) through the VDI Technologiezentrum as project management agency of the BMBF for the funding program FHprofUnt 2018 (funding code 13FH174PX8) based on a resolution of the German Bundestag and by Siemens Gas and Power GmbH & Co. KG financially supported H.G. also thanks the ‘Forschungsvereinigung der Arbeitsgemeinschaft der Eisen und Metall verarbeitenden Industrie e.V.’ (AVIF No. A316) for their financial support.

Notes

1 A refined analysis should model h as a function of the fluid boundary layer, which depends on temperature and mass flow of the working fluid.

Appendix

A.1. Workflow and used coefficients

Table A1. Workflow for optimizing a control a 3d model in FEM. Listed are the required work steps, the resulting file formats and the used programmes and tools.

Workstep	Used Program / Tool	File Format
3D-CAD Data Generation	Autodesk Inventor	.stl
Selecting Boundaries And calculation Area	gmsh	.geo
General Meshing	gmsh	.mesh
Convert Mesh For FEniCS	meshio (python)	.xmdf, .h5
Run FEM Simulation	FEniCS (python)
Run Optimization	Scipys SQP (python)
Show 3D Results	Paraview	.xmdf, .h5

Table A2. Used constant and coefficients.

Name	Abbreviation	Value and Unit
E-Modulus	E	210e9 $\frac{\mathrm{N}}{{\mathrm{m}}^2}$
Lame coefficient 1	λ	$\frac{\nu }{1-2\nu }\cdot \frac{E}{1+\nu }$
Lame coefficient 2	µ	$\frac{1}{2}\cdot \frac{E}{1+\nu }$
Poisson number	ν	0.3
Density	ρ	$8050\frac{\mathrm{kg}}{{\mathrm{m}}^3}$
Thermal expansion coefficient	α	13.5e-6 $\frac{1}{K}$
Specific heat	$c_p$	$420\frac{\mathrm{J}}{\mathrm{kGK}}$
Thermal conductivity	k	$36\frac{\mathrm{W}}{\mathrm{mK}}$
Convective heat coefficient	h	$20\frac{\mathrm{W}}{{\mathrm{m}}^2\mathrm{K}}$

A.2. Additional information about stress location

Table A3. Location of maximal van Mises stress.

Timestep	1	2	3	4	5	6	7	8	9	10
Node number	53	53	53	53	4	4	4	4	4	4
Timestep	11	12	13	14	15	16	17	18	19	20
Node number	4	4	4	4	4	4	4	53	53	53

Figure A1. Von Mises stress results for timestep 4 and 5. Location of maximal stress is moving between this points. This is important for our gradient based optimizer, otherwise it could result in jumps in stress. (a) Timestep 4 with maximal von Mises stress marked with dot and node number. (b) Timestep 5 with maximal von Mises stress marked with dot and node number.