Tutorial¶
This tutorial will walk through the process of writing a new user-defined mean-line solver and using it to run CFD on the Cambridge Wilkes3 cluster. We first need to download the turbigen code:
$ git clone https://gitlab.developers.cam.ac.uk/jb753/turbigen.git
$ cd turbigen
$ source setup.sh
Problem statement¶
Suppose we wish to design a rotor-only axial fan. We shall assume a constant axial velocity. The inlet state is specified as fixed values of \(T_{01}\) and \(p_{01}\) with no inlet swirl, \(\alpha_1=0\). We can then parametrise the aerodynamics of the stage using the following design variables:
Total pressure rise, \(\Delta p_0\)
Mass flow rate, \(\dot{m}\)
Flow coefficient, \(\phi=V_x/U\)
Loading coefficient, \(\psi= \Delta h_0/U^2\)
Hub-to-tip ratio, \(\mathit{HTR}=r_\mathrm{hub}/r_\mathrm{cas}\)
Total-to-total isentropic efficiency guess, \(\eta_\mathrm{tt}\)
To proceed with the annulus and blade shape design, turbigen requires the following quantities to be calculated from the above information. At the inlet and exit of the blade row:
Mean radii, \(r_\mathrm{rms}\)
Annulus areas, \(A\)
Absolute frame velocity vectors, \((V_x, V_r, V_\theta)\)
Static thermodynamic states, \((P, T)\) for example
Rotor shaft speed, \(\Omega\)
\(r_\mathrm{rms}\), \(A\), \(\Omega\) and the thermodynamic states should all be arrays of length 2. The velocity vectors should be of shape (3,2).
Mean-line design equations¶
We need to combine conservation of mass, momentum, and energy with definitions of our design variables to solve for the flow in the turbomachine. This will yield a set of equations that we will later implement numerically in the forward function.
The specified total pressure rise and guess of total-to-total efficiency allow calculation of the compressor work \(\Delta h_0 = h_{02}-h_{01}\)
where \(h_{02s}=h(p_{01}+\Delta p_0, s_1)\) is the ideal exit stagnation enthalpy for isentropic compression, i.e. enthalpy evaluated at the real exit stagnation pressure and inlet entropy.
The blade speed \(U\) is then given from the definition of loading coefficient
The definition of flow coefficient yields the axial velocity \(V_x\)
Assuming no inlet swirl, \(V_{\theta 1}=0\) and the Euler work equation yields the rotor exit circumferential velocity
We now have stagnation thermodynamic states and velocity vectors at inlet and exit of the rotor. If we knew the equation of state of the working fluid, we could evaluate static thermodynamic states and quantities such as density.
Conservation of mass gives the annulus area
and further specifying a hub-to-tip ratio fixes the mean radius
Finally, the shaft angular velocity is simply
Setting up skeleton files¶
To integrate our new mean-line design into turbigen, we have two
functions to write: a forward function which takes our design variables as
inputs and returns a turbigen.meanline.MeanLine object; and an
inverse function that recalculates the design variables from an input
turbigen.meanline.MeanLine object. Now that we know what input and
output data are required, and have the equations for our algorithm, we can start writing the functions. In a new file
called fan.py, copy and paste these definitions:
import turbigen.flowfield
import numpy as np
def forward(So1, DPo, mdot, phi, psi, htr, etatt):
"""Caluclate mean-line from inlet and design variables."""
# Insert code to calculate rrms, A, Omega, Vxrt, states
# ...
raise NotImplementedError
# Return assembled mean-line object
return turbigen.flowfield.make_mean_line(
rrms, # Mean radii
A, # Annulus areas
Omega, # Shaft angular velocity
Vxrt, # Velocity vectors
S # Thermodynamic states
)
def inverse(ml):
"""Calculate design variables from a mean-line object."""
# The output should be a dictionary keyed by the args to forward
return {
'So1': ml.stagnation[0],
# 'DPo': ...,
# 'mdot': ...,
# 'phi': ...,
# 'psi': ...,
# 'htr': ...,
# 'etatt': ...,
}
So1 is a fluid object that encapsualtes the inlet stagnation thermodynamic state. All
thermodynamic properties can be accessed as attributes, and there are functions
to manipulate the state to new values, described fully in turbigen.fluid .
We also need a minimal configuration file to test our mean-line functions. Create a new config.yaml with the following content:
# All files relating to the case are held in a working directory
workdir: runs/fan
# Perfect gas inlet state
inlet:
Po: 1e5
To: 300.
cp: 1005.
mu: 1.8e-5
gamma: 1.4
# Mean-line design
mean_line:
type: fan.py # Path to the mean-line module we are writing
# Our chosen design variables (args to forward)
DPo: 2000.
mdot: 5.
phi: 0.5
psi: 0.4
htr: 0.8
etatt: 0.9
At this point, running the config.yaml file through turbigen generates a NotImplementedError because the body of the forward function is missing.
Implementing the algorithm¶
We can now start to add the Mean-line design equations to the forward function inside fan.py.
The first task is to calculate the ideal exit enthalpy \(h_{02s}=h(p_{01}+\Delta p_0, s_1)\)
in Eqn. (1). Mean-line design functions should be written to make
no assumptions about the working fluid equation of state — this is accomplished
using the fluid modelling abstractions in turbigen.fluid. We take a
copy of the inlet state, and set its pressure and entropy to the required
values.
def forward(So1, DPo, mdot, phi, psi, htr, etatt):
"""Caluclate mean-line from inlet and design variables."""
# Get the ideal exit state
So2s = So1.copy() # Duplicate the inlet state
So2s.set_P_s(So1.P + DPo, So1.s) # Set pressure and entropy
# ...
We can now calculate the compressor work by reading off enthalpy values from our two state objects So1 and So2s.
# ...
# Work from defn efficiency Eqn. (1)
Dho = (So2s.h-So1.h)/etatt
# ...
Proceeding straightforwardly to calculate blade speed and velocity vectors
# ...
# Blade speed from defn psi Eqn. (2)
U = np.sqrt(Dho/psi)
# Axial velocity from defn phi Eqn. (3)
Vx = phi*U
# Circumferential velocity from Euler Eqn. (4)
Vt2 = Dho/U
# Assemble velocity vectors
# shape (3 directions, 2 stations)
Vxrt = np.stack(
(
(Vx, Vx), # Constant axial velocity
(0., 0.), # No radial velocity
(0., Vt2), # Zero inlet swirl
)
)
# ...
Next, we need to calculate the static thermodynamic states. As we know stagnation states and velocity vectors everywhere, this is most straightforward to do by evaluating the static enthalpy \(h=h_0-\frac{1}{2}V^2\). The static and stagnation states have the same entropy. In code, this looks like:
# ...
# Outlet stagnation state from known total rises
So2 = So1.copy().set_P_h(So1.P + DPo, So1.h + Dho)
# Assemble both stagnation states into a vector state
So = So1.stack((So1,So2))
# Get static states using velocity magnitude and same entropy
Vmag = np.sqrt(np.sum(Vxrt**2,axis=0))
h = So.h - 0.5*Vmag**2 # Static enthalpy
S = So.copy().set_h_s(h , So.s)
# ...
Now that the static states are known, the density can be used in the conservation of mass equation to continue with evaluating areas, the RMS radius, and the shaft angular velocity. The completed function is:
def forward(So1, DPo, mdot, phi, psi, htr, etatt):
"""Caluclate mean-line from inlet and design variables."""
# Get the ideal exit state
So2s = So1.copy() # Duplicate the inlet state
So2s.set_P_s(So1.P + DPo, So1.s) # Set pressure and entropy
# Work from defn efficiency Eqn. (1)
Dho = (So2s.h-So1.h)/etatt
# Blade speed from defn psi Eqn. (2)
U = np.sqrt(Dho/psi)
# Axial velocity from defn phi Eqn. (3)
Vx = phi*U
# Circumferential velocity from Euler Eqn. (4)
Vt2 = Dho/U
# Assemble velocity vectors
# shape (3 directions, 2 stations)
Vxrt = np.stack(
(
(Vx, Vx), # Constant axial velocity
(0., 0.), # No radial velocity
(0., Vt2), # Zero inlet swirl
)
)
# Outlet stagnation state from known total rises
So2 = So1.copy().set_P_h(So1.P + DPo, So1.h + Dho)
# Assemble both stagnation states into a vector state
So = So1.stack((So1,So2))
# Get static states using velocity magnitude and same entropy
Vmag = np.sqrt(np.sum(Vxrt**2,axis=0))
h = So.h - 0.5*Vmag**2 # Static enthalpy
S = So.copy().set_h_s(h , So.s)
# Conservation of mass for annulus area, Eqn. (5)
A = mdot/S.rho/Vx
# Mean radius from HTR Eqn. (6)
rrms = np.sqrt(A/np.pi/2.*(1.+htr**2)/(1.-htr**2))
# Shaft angular velocity
Omega = U / rrms
# Return assembled mean-line object
return turbigen.flowfield.make_mean_line(
rrms, # Mean radii
A, # Annulus areas
Omega, # Shaft angular velocity
Vxrt, # Velocity vectors
S # Thermodynamic states
)
This concludes the forward function — all the required quantities have been evaluated and can be returned for further processing.
Inverse function¶
If we run turbigen on the config.yaml file now, it will complete mean-line design successfully using the forward function, but raise an Exception because the inverse function is incomplete.
The inverse function serves as a verification check that the mean-line
matches the design intent, and also to extract design variables from a
mixed-out CFD solution. We add the design variables as keys in the output
dictionary, using the attributes of the turbigen.meanline.MeanLine
class to calculate them:
def inverse(ml):
"""Calculate design variables from a mean-line object."""
So1 = ml.stagnation[0]
So2s = So1.copy().set_P_s(ml.Po[-1], ml.s[0])
ho2s = So2s.h
return {
"So1": So1,
'DPo': ml.Po[-1] - ml.Po[0],
'mdot': ml.mdot[0],
'phi': ml.Vx[0]/ml.U[0],
'psi': (ml.ho[-1]-ml.ho[0])/(ml.U[0])**2,
'etatt': (ho2s-So1.h)/(ml.ho[-1]-ml.ho[0]),
'htr': ml.rhub[0]/ml.rtip[0]
}
Running CFD¶
We have now finished the mean-line design. To create blade shapes and run a computational fluid dynamics simulation, we can add some extra code to the config.yaml. These options are described fully in Configuration file format.
# All files relating to the case are held in a working directory
workdir: runs/fan
# Perfect gas inlet state
inlet:
Po: 1e5
To: 300.
cp: 1005.
mu: 1.8e-5
gamma: 1.4
# Mean-line design
mean_line:
type: fan.py # Path to the mean-line module we are writing
# Our chosen design variables (args to forward)
DPo: 2000.
mdot: 5.
phi: 0.5
psi: 0.4
htr: 0.8
etatt: 0.9
# ADD annulus configuration
annulus:
AR_gap: [1.0, 1.0] # Span to inlet/exit boundary distance
AR_chord: 3. # Span to chord
nozzle_ratio: 0.9 # Exit nozzle contraction
# ADD blade shapes
blades:
- DFL: 0.45 # Set number of blades using Lieblein
sections: # One blade section at midheight
- spf: 0.5
q_thick: [0.05, 0.12, 0.3, 0.02, 0.02, 0.18]
qstar_camber: [0., 0., 1.0, 1.0, 0.0]
# ADD mesh generation
mesh:
type: h # Mesh topology
yplus: 30.0 # Non-dimensional wall distance
resolution_factor: 0.5 # Use a coarse mesh
# ADD CFD solver
solver:
type: ts3 # Use Turbostream 3
nstep: 20000
nstep_avg: 5000
ilos: 1 # Mixing-length turbulence model
fmgrid: 0.4 # Full multigrid
# ADD control mass flow using a PID on exit pressure
operating_point:
mdot_pid: [0.5, 0.1, 0.0]
If we now run turbigen on our config.yaml using the shell command, we can quickly obtain a CFD solution for our newly designed fan.
$ turbigen config.yaml
TURBIGEN v1.5.1
Starting at 2024-01-29T13:57:10
Working directory: ...
Inlet: PerfectState(P=1.000 bar, T=300.0 K)
Designing a fan.py...
MeanLine(
Po=[1. 1.02] bar,
To=[300. 301.8913] K,
Ma=[0.099 0.127],
Vx=[34.5 34.5],
Vr=[0. 0.],
Vt=[ 0. 27.6],
Vt_rel=[-68.9 -41.4],
Al=[ 0. 38.66],
Al_rel=[-63.43 -50.19],
rpm=[2182. 2193.],
mdot=[5. 5.] kg/s
)
Checking mean-line conservation...
Checking mean-line inversion...
Annulus(
xmid=[-7.5253e-07 2.2099e-02],
rmid=[0.2999 0.2983],
span=[0.0666 0.0663]
)
Re_surf/10^5=[2.]
Nblade=[54], s_cm=[1.57], tip=[0.]
Generating an H-mesh...
Mesh Npts/10^6=0.11
Applying boundary conditions...
Setting intial guess...
Calculating wall distance...
Running solver ts3...
Using 1 GPUs on 1 nodes, 1 per node.
Checking convergence over last 5000 steps...
mdot drift = 0.0%, mdot error = -0.3%, eta_drift = -0.0%
Post-processing...
Mixed-out CFD result:
MeanLine(
Po=[0.9997 1.0115] bar,
To=[300.0521 301.2397] K,
Ma=[0.099 0.113],
Vx=[34.5 34.6],
Vr=[ 0.4 -0.8],
Vt=[ 0.8 18.2],
Vt_rel=[-68.1 -50.7],
Al=[ 1.35 27.72],
Al_rel=[-63.17 -55.69],
rpm=[2182. 2193.],
mdot=[5. 5.] kg/s
)
Design variable Nom CFD
------------------------------
DPo 2000.0 1172.7
etatt 0.9000 0.8433
htr 0.8000 0.8000
mdot 5.0000 4.9961
phi 0.5000 0.4997
psi 0.4000 0.2511
eta_tt=0.843, eta_ts=0.203
Elapsed time 1.40 min.
Creating and running designs with different velocity triangles is as simple as changing a line or two in the mean-line section of config.yaml. This allows us to explore a new design space very quickly.
Iterating the design¶
The table at the end of the program output compares the nominal mean-line design variables to actual values calculated using cuts from the three-dimensional CFD solution (the cuts are mixed out at constant area). Inspecting the output for our new fan, we can identify several problems:
Design variable Nom CFD
------------------------------
DPo 2000.0 1172.7 # Pressure rise too low
etatt 0.9000 0.8433 # Guessed efficiency too high
htr 0.8000 0.8000
mdot 5.0000 4.9961
phi 0.5000 0.4997
psi 0.4000 0.2511 # Loading too low
The root cause of the lack of pressure rise is that we have not allowed for deviation in designing the blade shapes, hence the flow is underturned. Assuming a guess of efficiency was neccesary to complete mean-line design, but its value should be updated so that the annulus areas are compatible with the intended velocity triangles.
Although it is not evident from the table, the inlet flow is not precisely aligned with the inlet metal angle, leading to unwanted accelerations around the leading edge. We should locate the stagnation point on the nose of the aerofoil to yield the smoothest pressure distributions.
turbigen has the capability to correct for all these issues. Adding an iterate key to the config.yaml will cause the program to repeatedly run the CFD, updating the efficiency guess and recambering the leading and trailing edges as needed:
# ...
# ADD new section for iterative corrections
iterate:
mean_line: # Correct efficiency guess
match_tolerance:
etatt: 0.01 # Efficiency to within 1%
relaxation_factor: 0.5 # Change is half CFD minus nominal
deviation: # Correct exit flow angles by TE recamber
clip: 5.0 # Maximum recamber in one step
relaxation_factor: 0.8 # Multiplier on changes to metal angle
tolerance: 1.0 # Absolute tolerance for termination in degrees
incidence: # Correct incidence by LE recamber
clip: 5.0 # Maximum recamber in one step
relaxation_factor: 0.2 # Multiplier on changes to metal angle
tolerance: 2.0 # Absolute tolerance for termination in degrees
Running the extended input file gives:
$ turbigen config.yaml
TURBIGEN v1.5.1
Starting at 2024-01-29T09:58:58
Working directory: ...
Iterating for max_iter=20 iterations
Min Inc DInc Dev DDev etatt Detatt
-------------------------------------------
1.401 6.271 1.254 -5.49 4.398 0.844 -0.002
1.401 4.064 0.812 -2.56 2.051 0.893 0.0232
1.401 2.502 0.500 -1.24 0.996 0.907 0.0182
1.401 1.491 0.298 -0.60 0.482 0.911 0.0115
1.401 0.885 0.177 -0.32 0.258 0.913 0.0067
1.401 0.497 0.099 -0.11 0.093 0.915 0.0041
Design variable Nom CFD
------------------------------
DPo 2000.0 1932.2
etatt 0.9111 0.9153
htr 0.8000 0.8000
mdot 5.0000 4.9941
phi 0.5000 0.4995
psi 0.4000 0.3832
eta_tt=0.915, eta_ts=0.383
Iteration finished in 8.4 min.
The corrections applied, DInc, DDev, and Detatt, decrease with each iteration indicating stable convergence. When the iteration terminates, the mixed-out CFD solution corresponds closely to the design intent. A new configuration file has been written out in the working directory runs/fan/config_conv.yaml for the converged solution. Inspecting this file:
# ...
qstar_camber:
- 3.142689298065078
- 8.280434755893827
- 1.0
- 1.0
- 0.0
# ...
Under the qstar_camber key that defines the camber line, we see that 3.1 degrees of recamber was required to align the stagnation point, and the deviation was 8.3 degrees. The efficiency has also been updated to 91.5%.
Extensions¶
This tutorial has demonstrated some of the functionality of turbigen. With the current choice of parameterisation, any change to the design is just an edit to the config.yaml, as described in Configuration file format.
Increase the number of blades by changing DFL
Increase the grid density under mesh
Control camber and thickness distributions by changing qthick and qstar_camber
Specify blade sections at multiple spanwise locations
Change the aspect ratio AR_chord
With a compatible CFD solver, change the working fluid to a real gas under inlet
To change the mean-line design, edit the forward and inverse functions in fan.py. For example: relax the assumption of constant axial velocity by adding a velocity ratio as one of the arguments to forward, replace specification of loading coefficient with a de Haller number, or specify an inlet Mach number instead of mass flow rate.
To add a stator, extend forward to take additional design variables and perform the necessary calculations. The output data should be at the inlet and exit of both blade rows, e.g. A an array of length 4, the velocity vectors should be of shape (3,4).