1. Introduction
To comply with the task of the reduction in carbon emissions, modern aircraft development currently follows (among others) two trends: on the one hand, the so-called more electric aircraft, which focuses on replacing hydraulic and/or pneumatic power distribution systems within the aircraft with an electric power distribution [
1], and on the other hand, high aspect ratio wings are used to reduce drag [
2]. For a fixed wingspan, a high aspect ratio reduces the wing’s thickness. Within the wing, complex systems to drive flaps, slats, and ailerons are positioned. The reduced space in the wings requires a highly compact design of these systems. Both trends result in the use of an electro-hydrostatic actuator (EHA) for primary flight control [
3].
Since an EHA is a compact system with high power density, thermal aspects are more relevant [
3,
4]. An EHA comprises a closed hydraulic loop. Therefore, compared to a conventional hydraulic actuator, heat cannot be dissipated from the system through the hydraulic lines. However, this closed hydraulic loop distributes the heat within the system [
5]. To develop an EHA for the application in primary flight control, it is therefore necessary to create a mechanical design that takes thermal aspects into account to ensure the full functionality of the EHA throughout all of the load cases [
6]. To leverage the full potential during the design phase, computer-aided optimization methods can support the developer to create an optimized design.
Computer-aided optimization methods are widely used in research and industry [
7]. There are different types of optimization methods, for example, parametric optimization [
8], genetic algorithms [
9], or shape and topology optimization [
10]. A general optimization problem can be formulated as follows:
where
f is the objective function,
g is an optional inequality constraint, and
h an optional equality constraint.
x denotes the design variables that are subject to variation and are in most cases bounded by set values. [
11] Many optimization algorithms for solving such problems work in iterations. This can lead to problems when each iteration requires a computationally expensive function evaluation [
12]. The flow of the hydraulic fluid within an EHA is inherently transient in most load cases. When applying iterative optimization algorithms to create designs for an EHA, the computational cost is extremely high as each function evaluation requires a time-consuming transient conjugate heat transfer (CHT) simulation.
CHT-simulation models are the state of the art when dealing with simulation based analysis of heat transfer in both fluid and solid regions [
13]. However, as already mentioned, they come at a high computational cost. Surrogate models can reduce the computational cost [
14]. They have been successfully applied to various engineering optimization problems, including but not limited to aerospace design [
15,
16] and heat exchanger optimization [
17,
18,
19], as well as coupled with neural networks for optimization problems [
20,
21].
The time scales of thermal conduction are much larger than those of convection [
22]. This results in a quasi-steady thermal state when considering periodic flows, such as those found in an EHA with oscillating movements. Many publications investigated the thermal behavior of oscillating flow. One of the fundamental works in the area of oscillating flow is a publication of Womersley, who investigated the oscillating flow in blood vessels [
23]. Named after him is the dimensionless Womersley number, defined as
where
is the oscillation’s angular frequency,
L is a characteristic length, and
is the fluid’s kinematic viscosity. The Womersley number gives insight into how the fluid behaves. For
, the oscillation frequency is relatively low. In this case, a flow has time to develop in both phases of the oscillation. For
, the frequency is relatively high, so that the fluid moves in a plug-like fashion within the pipe [
23].
Many researchers investigated the application of oscillating flow and heat transfer to stirling engines [
24,
25,
26]. There, the oscillation frequency was
and higher. For our application, frequencies around
and even lower are more common. Iwai et al. studied the heat transfer of a cylinder exposed to low-frequency oscillating flow [
27]. More recently, Bouvier et al. conducted an experimental study of a heated cylinder, in which air creates an oscillating flow [
28]. Choudari et al. conducted an extensive review of literature on oscillating flow heat transfer [
29].
In all these publications, the heat was entered either through the whole wall or through one end of the pipe in order to investigate and derive heat transfer coefficients along the pipe. In the present article, we investigate the heat distribution within zero-mean velocity oscillating flow coming from a local heat source.
The aim of this article is to derive a surrogate model to replicate the resulting quasi-steady thermal state of oscillating flow with a local heat source in a steady-state simulation model. This surrogate model is derived on a simplified geometry consisting of a straight pipe with oscillating flow and a local heat source within the fluid. In the context of an EHA, this resembles a periodic movement of the piston, while the hydraulic fluid oscillates through a heat source as, e.g., the pump. Applying this surrogate model to industrial thermal problems, where oscillating flow transports heat away from a heat source enables a time-efficient iterative optimization. This can support the product developer in generating suitable thermo-mechanical designs.
The article is structured into three sections.
Section 2 gives an overview of the proposed approach to derive a surrogate model for steady-state simulations.
Section 3 explains the details of the numerical simulation setup used within this article for both the transient as well as the steady-state simulations as well as the parametric optimization. In
Section 4, we present the results of the numerical simulations and evaluate the derived surrogate model.
2. Proposed Approach
The core idea of the proposed approach is to model the oscillating fluid’s effect on the temperature distribution throughout the pipe with a steady-state simulation model. The fluid moving through the volume where heat enters the system () transports this heat by means of advection as well as diffusion along the pipe. The periodic movement is inherently quasi-stationary, and after a certain time , the temperature distribution within the system is also quasi-stationary. The goal is to model the temperature distribution at the time averaged over one oscillation period.
To achieve this goal, we assume the fluid in the surrogate model to be stationary since the average fluid movement over one period equals zero. Afterwards, the fluid’s thermal conductivity and the length of the volume where heat enters the system are tuned through parametric optimization, so that the temperature distribution along the center line of the pipe equals the averaged values of the transient simulation. Results have shown that by using a single isotropic thermal conductivity in the fluid, the temperature distribution cannot be appropriately approximated. We therefore replaced the stationary fluid with a solid material with equal material properties (i.e., density, specific heat capacity, and thermal conductivity). This enables us to differentiate between a radial and an axial thermal conductivity and , respectively.
Figure 1 shows a graphical representation of the proposed approach, which will be explained in the following. A more detailed description of the individual steps taken for the practical application on the simplified pipe geometry is given in
Section 3.
At first, a case is defined for which a surrogate model shall be derived. This definition consists of the geometric dimensions of the pipe, specifically the total length , the inner (fluid) radius , the outer (solid) radius , and the width of the volume , where heat enters the system. In addition, the heat and the frequency f and amplitude A of the fluid movement need to be defined. This input data will further be referred to as input variables.
All input variables are used to set up a transient simulation of the oscillating fluid flow. After running the transient simulation, we extract the time series of the temperature along the pipe’s center line. This time series gets averaged over two oscillation periods to obtain a quasi steady-state temperature distribution from the transient simulation, which is written to a csv-file. The file is then imported into the basic simulation for the following parametric optimization as a reference for the optimization’s objective function. Since Simcenter STAR-CCM+ provides a Java interface, the extraction, averaging, and exporting to a file is done in a Java-macro that is executed at the end of the transient simulation. Importing the file into the basic simulation and adjusting the objective function is also done in a Java-macro.
The basic simulation, which contains the geometry data and heat source but not the frequency and amplitude since it is steady-state, is then imported into HEEDS to perform the parametric optimization. In addition to the basic simulation, bounds for the parameters to be optimized are needed. This is the last step where interaction from the user is necessary. The parametric optimization then runs until a stopping criterion is met. Within the parametric optimization, HEEDS runs a number of steady-state simulations with different values for
,
, and
. For each simulation, the resulting value of the objective function is compared to previous results. The parameters are then adjusted according to a proprietary algorithm (see [
31]) for a new number of steady-state simulations. The goal of the parametric optimization is to minimize the objective function given in Equation (
3) as presented in Equation (
1).
In this equation, the design variables consist of , , and . and describe the temperature at points along the pipe’s center line in the steady-state and transient simulation respectively. We investigated whether a normalization of the objective function is necessary. Results have shown that a normalization has no effect on the optimization result.
The result of using this approach is a set of parameters (, and , in the following referred to as output parameters). The resulting values for the orthotropic thermal conductivity can be interpreted as an effective thermal conductivity resulting from the combination of heat transfer through diffusion within the fluid and convective heat transfer from the fluid’s oscillating motion. To account for the limited extent of the volume where heat enters the system, the additional parameter needs to be adjusted. This set of parameters can then be applied by the product developer in a steady-state simulation to model the oscillating fluid’s influence on the temperature distribution in a quasi-stationary state. The resulting set of parameters is valid for one set of input variables. The method can, however, be applied to different sets of input variables to generate sets of output parameters for each individual case.
3. Simulation Setup
In this section, we present the setup of both transient and steady-state simulations. This comprises the geometric setup, meshing, boundary conditions, and involved physics. We also give a short overview of the parametric optimization and the used hard- and software.
All of the numerical simulations were carried out in Simcenter STAR-CCM+ version 2021.2.1. The optimization was carried out in HEEDS [
32]. Since our simulation models were relatively small, we were able to compute multiple simulations simultaneously on a 128-core 2 GHz server.
The points at which the temperature is evaluated for the objective function are equidistant points along the pipes center line with a distance of 1 cm.
Table 1 shows the input variables that were varied. Starting from a so called basic set of input variables marked in bold, each parameter was varied while keeping the other two at their base value, to investigate whether our proposed approach works for different geometries. This resulted in 13 separate sets of input variables that were investigated in this article. In addition to these varied variables, the remaining necessary input variables were kept constant throughout all cases; their values are as follows:
,
,
and
. The resulting Womersley number lies in the range of
. Note that
is part of the input variables as well as the output variables.
An overview of all combinations of input variables that are investigated is given in
Table 2. The resulting objective function values are discussed in
Section 4.
3.1. Geometry
Within this article, we derived surrogate models for oscillating flow within a straight pipe. We are evaluating our approach for different pipe geometries. The geometry, including the geometric parameters, is shown in
Figure 2. The geometry consists of the pipe itself and the contained fluid. The pipe has a total length
, an inner radius
, and an outer radius
. The pipe is filled with the fluid. In the center of the pipe, the region where heat enters the system is located. This region is a cylindrical shape with a length of
. The fluid radius
is part of the varied input variables. The corresponding values are listed in
Table 1.
3.2. Transient Simulations
The transient simulations have shown that a large physical time needs to be simulated in order to obtain a quasi-steady thermal state. Paired with the low enough time-step needed for a stable simulation, this results in high computational costs. In order to save computation time, the transient simulations were carried out on a 30 degree wedge of the fluid and the pipe. This assumption is valid since the case under consideration is of transverse symmetry. Periodic boundary conditions on the wedge’s sides ensure a proper representation of the pipe’s rotational symmetry. The boundary condition at the inlet has been adapted to replicate the full pipe’s mass flow.
The maximum Reynolds number was calculated for each set of input variables with the maximum velocity of the oscillation. With an overall maximum Reynolds number of , we are well below the laminar-turbulent transition region at around . Therefore, all simulations were set up with laminar flow.
The driving force of the oscillating flow is a mass flow boundary condition at one opening of the pipe. The second opening was set up as a pressure outlet, which enables a backflow into the pipe when a negative mass flow is prescribed at the first opening. The prescribed mass flow follows a sinusoidal function as defined in Equation (
4), where
is the density of the fluid and
t is the physical time.
This equation ensures that the oscillation of the fluid follows the given frequency and amplitude and oscillates around the pipes central region .
On the outer surface of the pipe, we prescribed an arbitrary fixed temperature of 20 °C. Due to constant material properties, this temperature serves as as a reference point only. The left and right surfaces of the pipe allow for no heat transfer out of the system. The inflowing fluid on both ends of the pipe also has a prescribed temperature of 20 °C.
The pipe and fluid both were meshed using an extruded mesh. Starting from one end of the pipe, a polyhedral mesh was extruded along the pipe and the fluid to obtain the full mesh. At the fluid’s boundary to the pipe, boundary layers were added. The cell size is approximately 1 mm perpendicular to the pipe’s direction and 2 mm along the pipe’s direction. The full mesh consists of around 38,500 cells for the transient basic simulation. For small values of , we reduced the cell size in order to properly discretize the region of the oscillating flow.
The time step was chosen to . This equals 100 time steps per oscillation period, which turned out to be sufficiently small for a stable simulation. Larger time steps resulted in non-physical simulation results. The stopping criterion was defined as the physical time of . This equaled 500 oscillation periods, which resulted in a quasi steady-state throughout all variations of the input variables.
3.3. Steady-State Simulations
The steady-state simulations were set up similarly to the transient simulations. If not mentioned otherwise in the following, the setup is the same as in the transient simulations. Even though the fluid part is modeled as a solid with the fluid’s altered material properties, we still refer to this region as the fluid region.
Compared to the transient simulations, the computational effort required to solve the steady-state simulations is relatively low (see
Figure 3). We therefore modeled the geometry of the full pipe. Since the fluid is modeled as a solid, its left and right sides were changed to walls. On these sides, a constant temperature of 20 °C is prescribed. The volume where heat enters the system is parameterized with the value
to allow for variation in the parametric optimization.
The pipe and the fluid were both meshed using a polyhedral mesh. The cell size was set to 2 mm, which results in an overall cell count of around 133,500 for the steady-state basic simulation.
The stopping criterion is a combination of the following three criteria. When all three criteria were met, the simulation was assumed to be converged:
The residual needs to be equal to or less than .
The standard deviation of the average temperature in the solid volume within the last 10 iterations needs to be equal to or less than .
The standard deviation of the average temperature in the fluid volume within the last 10 iterations needs to be equal to or less than .
Figure 3 shows the history of the residual and the volume-averaged temperatures of fluid and solid for the steady-state simulation of the basic simulation. The residual shows a logarithmic decrease, while the temperatures both converge. Due to the combination of a high thermal conductivity of the solid compared to the fluid and the fixed temperature boundary condition on the solid’s outer surface, the heat within the solid can not contribute to a large rise in temperature but is dissipated from the system through the boundary. Therefore, the solid shows almost no increase in average temperature.
As already mentioned in
Section 2, using an isotropic thermal conductivity for the surrogate models fluid did not yield adequate results. The oscillating fluid movement drastically increases the heat exchange throughout the fluid in the axial direction. Contrary to this, the radial heat exchange is barely affected by the fluid movement. Therefore, using an isotropic thermal conductivity results in either an overly large radial heat flux or an overly small axial heat flux with no value in between that satisfies our objective sufficiently. Since the fluid in the surrogate model is considered stationary, we can model this fluid as a solid part. This enables us to split the thermal conductivity into its axial and radial components. Both values can be tuned independently of one another to optimally parameterize the surrogate model.
3.4. Parametric Optimization
The optimization was carried out in HEEDS [
32]. HEEDS is a software for design space exploration and optimization. It includes the SHERPA optimization algorithm [
31]. The SHERPA algorithm is a proprietary algorithm that combines multiple optimization algorithms to find an optimal solution for a given optimization problem.
As mentioned in
Section 2, our goal is to find optimal parameters (
,
, and
) for a surrogate model to reproduce the same temperature distribution as in the underlying transient simulation model. The objective function was given in Equation (
3). Combining the objective function with the given parameter boundaries, the optimization formulation reads as follows:
For each set of input variables, the optimization was instructed to create 250 different sets of parameters to find an optimal solution. This amount resulted in a sufficient convergence of the optimization and a low value for the objective function. During the optimization, 15 sets of input parameters were run simultaneously before the next 15 sets of input parameters were defined. After the 250 simulation models are computed, the set of parameters with the lowest objective function value is chosen as the solution to the optimization problem.
4. Results and Discussion
In this section, we will show the results of applying our proposed approach to different sets of input variables. We will further discuss the results.
Figure 4 shows the temperature distribution from both the steady-state simulation as well as the transient simulation for the basic set of input variables. A symmetric temperature distribution can be seen. This is to be expected since the fluid is oscillating around the pipe’s center point at
. The curve is bell-shaped due to the fact that heat enters the system in the volume around the center point. The fluids’ oscillating movement as well as the heat conduction within the fluid then distributes this heat along the pipe. It is clearly visible that the temperature distribution resulting from the steady-state surrogate model almost perfectly replicates the time-averaged temperature distribution resulting from the transient simulation model. This shows that our proposed approach to construct a surrogate model for oscillating flow with a local heat source is feasible.
Table 2 shows the resulting values of the objective function for the considered cases. As can be seen, the final objective function value is generally low, except for cases 2 and 10. Case 2 has a higher objective function value, which is due to the small fluid radius of
. This small radius results in a Womersley number of
, meaning that the oscillating fluid movement exhibits characteristics of a fully developed flow during each half of the periodic cycle. Case 10 has a comparably small amplitude that results in a thermal hotspot. The current formulation of the surrogate model and its parameters is not capable of fully replicating this temperature profile.
Since the objective function is not normalized, the objective function value is generally a comparable indicator for the temperature deviations. The mostly low values of the objective function correspond to a good agreement between the averaged temperature distribution resulting from the transient simulation and the temperature distribution resulting from the surrogate steady-state simulation model.
Figure 5 and
Figure 6 show the resulting values for
over the frequency and the amplitude, respectively. Looking at the distribution of these discrete points, we suggest the assumption of a linear correlation of the radial thermal conductivity with the frequency and a quadratic correlation with the amplitude. These correlations are depicted by the dotted lines. Furthermore, it can be assumed that values of the radial thermal conductivity for frequencies and amplitudes that lie between the investigated points can be interpolated with the given linear and quadratic correlation. Care must be taken when interpolating near the extremes as, e.g.,
, which correlates to
, where the y-intercept of the trend line strongly affects the results of the surrogate model. Additionally, for
, the quadratic correlation for the thermal conductivity goes below zero, which is physically not possible. This must be taken into account when using these correlations for interpolations.
In all result plots, deviations between the discrete points and the trend line can be seen. This is due to the optimization process, which generally does not find the single optimal solution but rather a near-optimal solution.
Figure 7 and
Figure 8 show the results for
over the frequency and the amplitude, respectively. As in
Figure 5 and
Figure 6, the axial thermal conductivity shows a correlation with the amplitude. This time we suggest the assumption of an exponential correlation. When plotting
over the frequency, no clear correlation can be seen. The differences in axial thermal conductivity are quite large, with the maximum difference being
. A potential reason for this could be that the axial thermal conductivity is (nearly) independent of the frequency and has a small effect on the optimization’s objective function. This would theoretically result in a constant axial thermal conductivity over the frequency. Due to its partly evolutionary nature, the optimization algorithm, however, alters all parameters, including the axial thermal conductivity, and small numerical errors lead to certain values yielding slightly better optimization results (in our case lower objective function values). This is supported by the fact that the objective function value is sufficiently low for all cases, which indicates a successful optimization. We also observed a large change in the resulting axial thermal conductivity when slightly changing the models input parameters. However, the objective function value changes marginally, indicating a negligible influence of the axial thermal conductivity on the objective function and therefore on the resulting temperature distribution.
The same largely differing axial thermal conductivites can be observed when plotting the resulting parameter over the fluid radius
, as shown in
Figure 9.
Figure 10 seems to involve some form of correlation, although no low-order (first or second order) polynomial fits the points. Higher-order polynomials could create a fitting trend line, however overfitting might occur.
Two different regions can be identified: for , the radial thermal conductivity rises with an increasing fluid radius. For , the radial thermal conductivity decreases with an increasing fluid radius. This is due to a low Womersley number of around 1, which indicates that the oscillatory effects on the fluid flow are replaced by regular pipe flow in both phases of the oscillation.
Figure 11,
Figure 12 and
Figure 13 show the resulting
over the frequency, the amplitude, and the fluid radius, respectively. For
over
f and
, no clear correlation is recognizable. This is again due to the fact that the parameter
is theoretically independent of the frequency and the fluid radius. Numerical deviations in the simulation and optimization then prefer certain values over others, without a noticeable change of the surrogate model. Even though certain parameters show large variations following from small changes of the input variables, the surrogate model still replicates the temperature profile for the given set of input parameters.
For over A, however, a clear dependence can be observed. With an increase in amplitude, the value for increases as well. For , stagnates slightly below . This is due to our total pipe length set to . Therefore, cannot climb above that value. At a first glance, this might seem like an erroneous simulation and/or optimization setup as the given input parameters and boundaries restrict the optimization algorithm. However, the original goal of the surrogate model was to replicate the temperature distribution within the given geometry. When applying this method to create a surrogate model for industrial applications, one might encounter similar restrictions due to geometric or procedural parameters. The created surrogate model then replicates the resulting quasi-steady thermal state within the given restrictions of the oscillating transient flow.
5. Conclusions and Outlook
In this article, we presented an approach to derive a surrogate model for oscillating flow with a local heat source. We applied the approach to the generic case of a straight pipe with an internally oscillating fluid. This surrogate model can then be applied to steady-state simulations, where, due to the periodic oscillation, a quasi-steady thermal state develops. With the surrogate model, the temperature distribution resulting from the transient oscillating flow is reproduced. By using a steady-state simulation, the computational cost for thermal investigations comprising oscillating flow and a local heat source can drastically be reduced. By reducing the computational cost, these surrogate models can serve as the foundation for the iterative optimization of industrial applications.
Certain dependencies of the surrogate model’s parameters (axial and radial thermal conductivity and extent of the local heat source) on the input parameters (frequency, amplitude and fluid radius) could be observed:
The radial thermal conductivity depends on the fluid radius, the frequency, and the amplitude.
The axial thermal conductivity depends on the amplitude.
The extent of the local heat source depends on the amplitude.
We suggest the assumption that and do not depend on the fluid radius and the frequency and have a small effect on the objective function. This is due to the seemingly random variations over the fluid radius and frequency. These variations come from the (partly) evolutionary nature of the optimization algorithm. In the optimization process, all variables are varied simultaneously. Thereby, when an optimum is found, values that have a very small effect on the objective function can exhibit an arbitrary behavior.
Except for over the fluid radius, all of the dependencies can be approximated by a first-order polynomial or an exponential curve. By formulating these dependencies, combinations of input parameters that differ slightly from those that were investigated can be derived. However, for small values of the fluid radius and the frequency, the deviations between the low-order approximations and the calculated parameters increased. This is due to the fact that the Womersley number decreases to 1 and below, which indicates a fully developed flow in both phases of the oscillation. Small amplitudes of the oscillation also posed problems in finding optimal parameters for the surrogate model. For input parameters that differ largely from those investigated, it is advised to set up a new parameter fitting as presented in this article since geometrical aspects and effects based on too low or too high Womersley numbers affect the behavior of the parameters.
In future works, we will apply this method to the optimization of an EHA. Due to the EHA’s periodic movement, the internal hydraulic fluid creates an oscillating flow that distributes the heat from local heat sources, such as the motor or the pump, throughout the system. Difficulties might arise when applying the proposed approach to non-linear pipe geometries. Nevertheless, this should only increase the effort required to set up the simulation models. Based on the results from this article, it is reasonable to conclude that the proposed approach can successfully be applied to other geometries, such as the internal hydraulics of an EHA, as well. After deriving a surrogate model for the thermal behavior of the oscillating flow, optimization methods that require numerous function evaluations, such as the method of topology optimization or design space exploration, can then be applied with reasonable computational cost. This opens up further potential for the creation of mechanical designs that are thermally optimized.