diff --git a/README.md b/README.md index 30683b3e..f6bb3996 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,6 @@ # su2code.github.io [Link to Website](https://su2code.github.io/) + +## For Developers + +In order to make any changes in the documentation of SU2, the files in [_docs_v7/](_docs_v7) should be touched. This is the documentation which is currently online. diff --git a/_data/tutorials.yml b/_data/tutorials.yml index 9c604bd8..76c6f2a0 100644 --- a/_data/tutorials.yml +++ b/_data/tutorials.yml @@ -17,6 +17,7 @@ - Unsteady_NACA0012 - UQ_NACA0012 - NICFD_nozzle + - Aachen_Turbine - title: Incompressible Flow tutorials: @@ -29,6 +30,7 @@ - Inc_Streamwise_Periodic - Inc_Species_Transport - Inc_Species_Transport_Composition_Dependent_Model + - Inc_Von_Karman - Inc_Turbulent_Bend - title: Structural Mechanics @@ -55,6 +57,7 @@ - Multi_Objective_Shape_Design - Unsteady_Shape_Opt_NACA0012 - Species_Transport + - Inc_Turbulent_Bend_Opt - title: Workflow Setup tutorials: diff --git a/_docs_v7/Developing-SU2-on-GitHub-(Internal-Developers).md b/_docs_v7/Developing-SU2-on-GitHub-(Internal-Developers).md index 97b90519..f14cfdd6 100644 --- a/_docs_v7/Developing-SU2-on-GitHub-(Internal-Developers).md +++ b/_docs_v7/Developing-SU2-on-GitHub-(Internal-Developers).md @@ -5,10 +5,10 @@ permalink: /docs_v7/Developing-SU2-on-GitHub-(Internal-Developers)/ The repository for SU2 is being hosted here on GitHub. As you are likely aware, GitHub is simply an online project hosting service with a very useful web interface and additional tools to aid code development with Git as its backbone. Git is a version control system (VCS) that is similar to SVN, Mercurial, etc., and it helps organize the development of code over time by tracking changes. -To get started, you need to create a personal user account on GitHub (free) and follow the [basic setup instructions](https://help.github.com/articles/set-up-git). These instructions include how to get Git installed on your local machine. To sync up your local settings with GitHub, change the user.email and user.name variables for your local git configuration with +To get started, you need to create a personal user account on GitHub (free) and follow the [basic setup instructions](https://help.github.com/articles/set-up-git). These instructions include how to get Git installed on your local machine. To sync up your local settings with GitHub, change the `user.email` and `user.name` variables for your local git configuration with ``` -git config --global user.email "your_email@domain.com" -git config --global user.name "Your Name" +git config --local user.email "your_email@domain.com" +git config --local user.name "Your Name" ``` Note that the email address should be the one associated with your GitHub account. @@ -22,6 +22,7 @@ After cloning, you should have a new SU2/ folder in your current working directo ``` git log ``` +To setup the local copy of SU2 for development purposes, one must follow the steps mentioned in [Build-SU2-Windows](_docs_v7/Build-SU2-Windows.md) and [Build-SU2-Linux-MacOS](_docs_v7/Build-SU2-Linux-MacOS.md). ## Typical Workflow with Git diff --git a/_tutorials/compressible_flow/Aachen_Turbine/Aachen_Turbine.md b/_tutorials/compressible_flow/Aachen_Turbine/Aachen_Turbine.md new file mode 100644 index 00000000..c59324f2 --- /dev/null +++ b/_tutorials/compressible_flow/Aachen_Turbine/Aachen_Turbine.md @@ -0,0 +1,392 @@ +--- +title: Aachen turbine stage with Mixing-plane +permalink: /tutorials/Aachen_Turbine/ +written_by: alecappiello +for_version: 8.1.0 +revised_by: joshkellyjak +revision_date: +revised_version: 8.1.0 +solver: RANS +requires: SU2_CFD +complexity: advanced +follows: compressible flow tutorials +--- + +![Mid-span Rel. Mach](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/FullAnnulus_MachContour.png) + +## Goals +Upon completing this tutorial, the user will become familiar with 3D steady-state RANS calculations of multi-row axial turbines, elaborating viscous, compressible flows of air, modelled as an ideal gas. + + +The geometry chosen for the tutorial is a 1-1/2 axial turbine stage (stator 1 - rotor - stator 2), from Gallus, H. E., et al. "*Endwall and Unsteady Flow Phenomena in an Axial Turbine Stage.*" ASME. J. Turbomach. October 1995; 117(4): 562–570. https://doi.org/10.1115/1.2836568. + + +The solution will provide a flow field that can be compared against experimental data. However, to achieve a meaningful result, much more refined grids should be used. + +*Note that the original stator geometries have been modified, increasing the blade count from 36 to 41.* + + +The intent of this tutorial is to introduce a compressible flow problem involving an axial turbine to explain how turbo-specific features are used within SU2. +The following capabilities of SU2 will be showcased in this tutorial: +- Steady, 3D RANS equations with the Spalart–Allmaras (SA) turbulence model +- Giles non-reflective boundary conditions +- Non-reflective mixing-plane +- Hub and Shroud surface treatment +- Turbo Performance BCs +- Wall functions + + +## Resources + +You can find the resources for this tutorial in the folder [compressible_flow/Aachen_Turbine](https://github.com/su2code/Tutorials/tree/master/compressible_flow/Aachen_Turbine) in the [tutorial repository](https://github.com/su2code/Tutorials). +You will need the mesh file [Aachen_Turbine.su2](https://github.com/su2code/Tutorials/tree/master/compressible_flow/Aachen_Turbine/Aachen_Turbine.su2) and the main config file [Aachen_Turbine.cfg](https://github.com/su2code/Tutorials/tree/master/compressible_flow/Aachen_Turbine/Aachen_Turbine.cfg), together with the sub-config files: [stator1.cfg](https://github.com/su2code/Tutorials/tree/master/compressible_flow/Aachen_Turbine/stator1.cfg), [rotor.cfg](https://github.com/su2code/Tutorials/tree/master/compressible_flow/Aachen_Turbine/rotor.cfg), [stator2.cfg](https://github.com/su2code/Tutorials/tree/master/compressible_flow/Aachen_Turbine/stator2.cfg). + +*Note that, while the mesh used for this tutorial feature boundary layer refinements at the endwalls and on blade surfaces, +it does not include the rotor tip gap and it is rather coarse. For comparison of the results against literature, a finer meshes should be used.* + + + +## Tutorial + +The following tutorial will walk you through the steps required when solving state-state flows through multi-row turbines using SU2. It is assumed you have already obtained and compiled SU2_CFD. If you have yet to complete these requirements, please see the [Download](/docs_v7/Download/) and [Installation](/docs_v7/Installation/) pages. + +### Background + +This example uses a 3D one and a half turbine stage encompassing one stator, one rotor and a downstream stator, equal to the +first one. Consequently, the case requires multiple frame of reference to account for the rotation of the rotor. + + +### Problem Setup + +This tutorial will solve the flow through the Aachen turbine for the following boundary conditions: +- Working fluid is air +- Inlet Stagnation Temperature = 308.26 K +- Inlet Stagnation Pressure = 158245.38 Pa +- Inlet Flow Direction, unit vector (Flow dir-norm, Flow dir-tang, Flow dir-span) = (1.0, 0.0, 0.0) +- Outlet Static Pressure = 110050.96 Pa +- Rotor Rotational Speed, vector component (x y z) = 0.0 0.0 -366.52 rad/s + + +### Geometry Description + +To reduce the computational cost, the geometrical periodicity of the system is exploited reducing the computational domain to one blade passage per row +(with the blade at the centre of the mesh). As it will be explained later, this requires the use of periodic boundary conditions. + +![Full Annulus Domain](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/FullAnnulus.png) +![Single Passage Domain](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/SinglePassage.png) +Figure (1): Computational domain reduction. + +### Mesh Description +The mesh is composed of hexahedral elements, with 53634 elements for stator 1 passage, 67935 elements for the rotor passage and 52374 for stator 2 passage. + +*Upon preparing the mesh, it is important to make sure that the machine axis is the Z-axis and that the flow proceeds in the positive z-direction.* + + +### Multizone Config File +The current case features a multizone domain. As a consequence, the `MULTIZONE` option must be used, and the list of sub-config files (one per zone) must be provided via `CONFIG_LIST`. An example of the necessary syntax is provided in the following box. + +``` +MULTIZONE= YES + +% List of config files for zone-specific options +CONFIG_LIST=(stator1.cfg, rotor.cfg, stator2.cfg) +``` +#### Sub-config Files +The most important information that must be specified in the sub-config files concerns the grid movement. To activate it, the `GRID_MOVEMENT` option must be specified, together with the `MOTION_ORIGIN` and angular velocity (`ROTATION_RATE`) components about the rotation axis.\ +For the solver initialization the motion Mach number is specified via the `MACH_MOTION`. \ +Finally, a smoother convergence can be achieved via the `RAMP_ROTATING_FRAME`. Although not used in the present case, it can be activated similarly to the outlet pressure ramp, describe later in this tutorial. + +These options are shown for the rotor sub-config file in the following box. + +``` +% Type of dynamic mesh (NONE, ROTATING_FRAME) +GRID_MOVEMENT = ROTATING_FRAME + +MACH_MOTION= 0.35 + +MOTION_ORIGIN= 0.0 0.0 0.0 + +% Angular velocity vector in rad/s +ROTATION_RATE= 0.0 0.0 -366.52 + +RAMP_ROTATING_FRAME= NO + +% Parameters of the rotating frame ramp: + starting rotational speed, + updating-iteration-frequency, + total number of iterations for the ramp +RAMP_ROTATING_FRAME_COEFF= (0.0, 150.0, 2000) +``` + +### Boundary Conditions +To setup a turbomachinery CFD calculation, several kinds of boundary conditions are needed. The boundary conditions relevant to the Aachen turbine test case are discussed in the following, together with the necessary syntax to activate them. + +#### Adiabatic, no-slip boundary conditions: + +This kind of boundary condition is applied to all solid walls, highlighted in grey in Figure 2, such as HUB, SHROUD +AND BLADE belonging to each row. To make a wall adiabatic the ```MARKER_HEATFLUX``` can be used specifying a zero heatflux. This will also ensure that a no-slip condition is enforced on that wall. +The following box shows the syntax to be used to enforce this boundary condition at all solid walls present in the meshes. + +``` +% Format: (marker, 0.0) + +MARKER_HEATFLUX = (BLADE1, 0.0, BLADE2, 0.0, BLADE3, 0.0, HUB1, 0.0, SHROUD1, 0.0, HUB2, 0.0, SHROUD2, 0.0, HUB3, 0.0, SHROUD3, 0.0) +``` +![Adiabatic, no-slip bc](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_NoSlip.png) +Figure (2): Hub, shroud and blade surfaces. + + +#### (Rotational) Periodic boundary conditions: +This kind of boundary condition is applied to pairs of surfaces, such as the sides of the domain +(highlighted in Figure 3 in brown and purple). To enforce it, ```MARKER_PERIODIC``` should be used, providing the name of a "periodic marker" (purple surface), a "donor marker" (brown surface). +In cases involving rotational periodicity (in contrast to translational periodicity) the coordinates of the rotation centre (x,y,z) and the rotation angle components (x-axis, y-axis, z-axis) should be provided. +The following box shows the syntax to be used to enforce this boundary condition at all sold walls present in the meshes. +``` + +% Format: ( periodic marker, donor marker, rot_cen_x, rot_cen_y, rot_cen_z, rot_angle_x-axis, rot_angle_y-axis, rot_angle_z-axis, translation_x, translation_y, translation_z) + +MARKER_PERIODIC = (PER1_STATOR1, PER2_STATOR1, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7804878, 0.0, 0.0, 0.0, PER1_ROTOR, PER2_ROTOR, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7804878, 0.0, 0.0, 0.0, PER1_STATOR2, PER2_STATOR2, 0.0, 0.0, 0.0, 0.0, 0.0, 8.7804878, 0.0, 0.0, 0.0) + +``` +![Periodic boundary conditions](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Periodics.png) +Figure (3): Periodic and donor surfaces. + +#### Non-reflective boundary conditions +Non-reflective boundary conditions can be enforced by means of the +```MARKER_GILES``` boundary. They are used for both inlet-outlet, as well as mixing-plane boundaries. Furthermore, under-relaxation factors can be provided both for the average and Fourier components at any boundary where the Giles boundary condition is used. For additional details about their implementation and usage, the reader is referred to the paper by [Vitale et al.](https://doi.org/10.2514/1.B37685) + +##### Inlet and outlet boundaries +The ```TOTAL_CONDITIONS_PT``` option allows to enforce the total inlet condition on a chosen boundary (identified by its marker). These include pressure and temperature, together with the incoming flow direction, expressed via its local coordinates (normal, tangential and radial to the chosen boundary). +The ```STATIC_PRESSURE_1D``` can be used to enforce the outlet static pressure value on a boundary where the flow is in radial equilibrium. +These two options are applied to the inlet and outlet surfaces of the domain, respectively, highlighted in light blue in Figure 4. +The following box presents the corresponding syntax. + +``` +% Non reflecting boundary condition for inflow and outflow: + +% Format inlet: +( marker, TOTAL_CONDITIONS_PT, Total Pressure , Total Temperature, Flow dir-norm, Flow dir-tang, Flow dir-span, under-relax-avg, under-relax-fourier) + +(INFLOW_STATOR1, TOTAL_CONDITIONS_PT, 158245.38, 308.26, 1.0, 0.0, 0.0, 0.3, 0.0) + + +% Format outlet: +(marker, STATIC_PRESSURE_1D, Static Pressure value, -, -, -, -, under-relax-avg, under-relax-fourier) + +(OUTFLOW_STATOR2, STATIC_PRESSURE_1D, 110050.96, +0.0, 0.0, 0.0, 0.0 , 1.0, 0.0) +``` + +![Non-reflective inlet and outlet boundary conditions](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Giles.png) +Figure (4): Non-reflective inlet and outlet boundary conditions + +###### Mixing-plane surfaces +For mixing plane surfaces, as those highlight in light red in Figure 5, ```MIXING_IN``` and ```MIXING_OUT``` tags are used, depending on whether the flow enters or exits through that boundary. +The following box shows the corresponding syntax. + +``` +% Non reflecting boundary condition for mixing-plane: +% Format mixing-plane in and out: +( marker, MIXING_IN or MIXING_OUT, -, -, -, -, -, -, under-relax-avg, under-relax-fourier) + +(OUTFLOW_STATOR1, MIXING_OUT, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, +INFLOW_ROTOR, MIXING_IN, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0) +``` + +Additionally, to activate the transfer of information between zones at the mixing planes, couples of boundaries forming a mixing-plane must be specified via the ```MARKER_MIXINGPLANE_INTERFACE``` and ```MARKER_ZONE_INTERFACE``` in stream-wise order, as shown in the following box. + +``` +MARKER_MIXINGPLANE_INTERFACE= (OUTFLOW_STATOR1, INFLOW_ROTOR, OUTFLOW_ROTOR, INFLOW_STATOR2) + +MARKER_ZONE_INTERFACE= (OUTFLOW_STATOR1, INFLOW_ROTOR, OUTFLOW_ROTOR, INFLOW_STATOR2) +``` + + + +The kind of interpolation method to be used at the mixing-planes can be specified via the ```MIXINGPLANE_INTERFACE_KIND```. Possible options are: ```LINEAR_INTERPOLATION```, ```NEAREST_SPAN```, ```MATCHING```. + +``` +MIXINGPLANE_INTERFACE_KIND= LINEAR_INTERPOLATION +``` +To make the mixing-plane turbulent, the following option must be used: +``` +TURBULENT_MIXINGPLANE= YES +``` + +![Non-reflective mixing-plane](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Giles_MixingPlane.png) +Figure (5): Non-reflective mixing-plane + +To turn on the non-reflective behaviour, the `SPATIAL_FOURIER` should be set to `YES`. However, for the current tutorial it is not activated, as in the following box. +``` +SPATIAL_FOURIER = NO +``` +If needed, extra under relaxation factor for the Giles BC at the hub and shroud can be provided as follows: +``` +GILES_EXTRA_RELAXFACTOR = (0.05, 0.05) +``` + +Finally, the kind of average process for linearizing the Navier-Stokes equation at inflow and outflow boundaries, including mixing-plane(s), can be controlled by ```AVERAGE_PROCESS_KIND```. Possible options are: ```ALGEBRAIC```, ```AREA```, ```MASSSFLUX```, ```MIXEDOUT```. The default one is ```AREA```, however, only the ```MIXEDOUT``` option guarantees that the mixing plane treatment is conservative.\ +The following box shows the setting to be used for the current case. + +``` +AVERAGE_PROCESS_KIND= MIXEDOUT +``` + + +#### Shroud boundary condition +The ```MARKER_SHROUD``` allows to specify the mesh boundaries belonging to the shroud, as those highlighted in dark red in Figure 6. +When a boundary is tagged as part of the shroud, its grid velocity components are set to 0, mimicking the physical behaviour of a shroud wall. Only rotating zones, such as a rotor, are affected. +The following box shows the syntax to enforce this functionality. +``` +MARKER_SHROUD = (SHROUD1, SHROUD2, SHROUD3) +``` + +![Shroud boundary condition](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_shroud.png) +Figure (6): Shroud boundary condition + +#### Wall functions +To guarantee a low computational cost, a particularly coarse mesh, featuring large $y^{+}$ values, was used. Therefore, to achieve convergence of the present calculation, wall functions are needed. These are activated by the ```MARKER_WALL_FUNCTIONS```, followed by the ```MARKER_TAG``` of the boundaries to which the wall functions will be applied. +The following box shows the basic syntax to active the wall functions for the current case. +``` +% format: +% MARKER_WALL_FUNCTIONS = (MARKER_TAG, STANDARD_WALL_FUNCTION) + +MARKER_WALL_FUNCTIONS= ( BLADE1, STANDARD_WALL_FUNCTION, + BLADE2, STANDARD_WALL_FUNCTION, + BLADE3, STANDARD_WALL_FUNCTION, + HUB1, STANDARD_WALL_FUNCTION, + SHROUD1, STANDARD_WALL_FUNCTION, + HUB2, STANDARD_WALL_FUNCTION, + SHROUD2, STANDARD_WALL_FUNCTION, + HUB3, STANDARD_WALL_FUNCTION, + SHROUD3, STANDARD_WALL_FUNCTION) + +WALLMODEL_KAPPA= 0.41 +WALLMODEL_B= 5.5 +WALLMODEL_MINYPLUS= 5.0 +WALLMODEL_MAXITER= 200 +WALLMODEL_RELFAC= 0.5 +``` + + +#### Ramps + +To ease the convergence of turbomachinery cases, particularly when the pressure ratio of the machine is large, the use of ramps can be beneficial. These can be applied to both rotational speed and outlet pressure via ```RAMP_OUTLET_PRESSURE``` and ```RAMP_ROTATING_FRAME```. \ +For the sake of completeness, the use of the outlet pressure ramp is demonstrated presenting in the following box the syntax to activate it. + +``` +% Specify ramp option for Outlet pressure (YES, NO) default NO +RAMP_OUTLET_PRESSURE= YES +% +% Parameters of the outlet pressure ramp (starting outlet pressure, updating-iteration-frequency, total number of iteration for the ramp) +RAMP_OUTLET_PRESSURE_COEFF= (140000.0, 150.0, 2000) +``` + + + + +#### Turbo Performance and Analysis +To have SU2_CFD computing the typical turbomachinery performance indexes, the turbomachinery layout and machine kind must be specified by means of ```TURBOMACHINERY_KIND``` and ```TURBO_PERF_KIND```. The former lets you specify whether the machine is axial, centripetal, centrifugal, etc., while the latter lets you specify whether the machine is a compressor or a turbine. These data must be specified for each zone of which the machine is made up of. + +Additionally, the inlet and outlet boundary to each zone must be specified in stream-wise order via the ```MARKER_TURBOMACHINERY```, as well as the inlet and the outlet boundaries of the full machine via the ```MARKER_ANALYZE```. + +Finally, the performance averaging method, i.e. AREA, MASSSFLUX, MIXEDOUT, etc., must be specified via the ```PERFORMANCE_AVERAGE_PROCESS_KIND```. In case the MIXEDOUT option is selected, the parameters for the Newton method used for the mixed out process can be specified via ```MIXEDOUT_COEFF```. + +The following box summarizes the setup for the Turbo Performance and Analysis of the current case. + +``` +% Specify kind of architecture (AXIAL, CENTRIPETAL, CENTRIFUGAL, CENTRIPETAL_AXIAL) +TURBOMACHINERY_KIND= AXIAL AXIAL AXIAL + +TURBO_PERF_KIND= (TURBINE, TURBINE, TURBINE) + +MARKER_TURBOMACHINERY= (INFLOW_STATOR1, OUTFLOW_STATOR1, + INFLOW_ROTOR, OUTFLOW_ROTOR, + INFLOW_STATOR2, OUTFLOW_STATOR2) + +MARKER_ANALYZE = (INFLOW_STATOR1, OUTFLOW_STATOR2) + +% (ALGEBRAIC, AREA, MASSSFLUX, MIXEDOUT) default AREA +PERFORMANCE_AVERAGE_PROCESS_KIND= MIXEDOUT + +%Parameters of the Newton method for the MIXEDOUT average algorithm (under relaxation factor, tolerance, max number of iterations) +MIXEDOUT_COEFF= (1.0, 1.0E-05, 15) +``` + + +#### Plotting + +You can choose which boundaries to include in the surface flow file written by SU2_CFD via ```MARKER_PLOTTING```. +The following box shows the setting to export the three blade surfaces to the surface flow solution file. + +``` +% Marker(s) of the surface in the surface flow solution file +MARKER_PLOTTING= (BLADE1, BLADE2, BLADE3) +``` + + +### Configuration File Options +The initialization of the flow solver can be performed from thermodynamic quantities, directly with `INIT_OPTION=TD_CONDITIONS`. It is recommended that you always confirm the resulting initialization state in the console output during runtime of SU2 that is reported just before the solver begins iterating. +The initialization options are specified in the following box. + +``` +% Free-stream pressure (101325.0 N/m^2 by default, only Euler flows) +FREESTREAM_PRESSURE= 140000.0 +% +% Free-stream temperature (273.15 K by default) +FREESTREAM_TEMPERATURE= 300.0 +% +% Free-stream temperature (1.2886 Kg/m3 by default) +FREESTREAM_DENSITY= 1.7418 +% +% Free-stream option to choose if you want to use Density (DENSITY_FS) or Temperature TEMPERATURE_FS) to initialize the solution +FREESTREAM_OPTION= TEMPERATURE_FS +% +% Free-stream Turbulence Intensity +FREESTREAM_TURBULENCEINTENSITY = 0.025 +% +% Free-stream Turbulent to Laminar viscosity ratio +FREESTREAM_TURB2LAMVISCRATIO = 100.0 +% +% +%Init option to choose between Reynolds (default) or thermodynamics quantities for initializing the solution (REYNOLDS, TD_CONDITIONS) +INIT_OPTION= TD_CONDITIONS +``` + +#### Final Remarks + +To have the residuals exported from each volume zone `WRT_ZONE_HIST` should be set to `YES`. + +### Running SU2 + +To run this test case, the standard procedure can be followed. At a terminal command line: + +1. Move to the directory containing the config files (Aachen_Turbine.cfg, stator1.cfg, rotor.cfg, stator2.cfg) and the mesh (Aachen_Turbine.su2) file. Make sure that the SU2 tools were compiled, installed, and that their install location was added to your path. +2. Run the executable by entering in the command line: + + ``` + $ SU2_CFD Aachen_Turbine.cfg + ``` + +3. SU2 will print residual updates with each iteration of the flow solver, and the simulation will terminate after reaching the specified convergence criteria. +4. Files containing the results will be written upon exiting SU2. The flow solution can be visualized in ParaView (.vtk) or Tecplot (.dat for ASCII), depending on the user requirements in the config file. + +### Results + +Figure 7a and b show the relative Mach number distributions on a cylindrical cut slicing the flow domain at approximately mid-span. The resulting cylindrical surface is presented in Fig. 7a, while it 2D planar counter part is shown in Fig.7b for an easier visualization. + +To compute it the following steps can be take in a post-processing tool, such as Tecplot or ParaView: + +1. Compute Radius $R$, angualr coordinate $\theta$, and their product $R \cdot \theta$ + +2. Compute the relative velocity components and the relative Mach number + +3. Slice the flow volume at $R = const$ and plot the relative Mach number on the corresponding iso-surface at $R = const$ (Fig. 7a) + +3. Extract the slice + +4. Plot the relative Mach number contour of the extracted slice in a 2D space $(R \cdot \theta$ vs $z)$ and reverse the Z-axis, Fig. 7b. + +![Midspan Relative Mach number](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/Full_View_MachContour.png) +![2D Relative Mach number blade-to-blade view](../../../tutorials_files/compressible_flow/Aachen_Turbine/images/2D_MidSpan_RelMach.png) +Figure (7): (a) Relative Mach number contour on a cylindrical cut (blade-to-blade view) at midspan of the turbine flow domain; (b) planar (blade-to-blade) view of the cylindrical cut. diff --git a/_tutorials/design_features/Inc_Turbulent_Bend_Opt/Inc_Turbulent_Bend_Opt.md b/_tutorials/design_features/Inc_Turbulent_Bend_Opt/Inc_Turbulent_Bend_Opt.md new file mode 100644 index 00000000..bf43128f --- /dev/null +++ b/_tutorials/design_features/Inc_Turbulent_Bend_Opt/Inc_Turbulent_Bend_Opt.md @@ -0,0 +1,218 @@ +--- +title: Adjoint design optimization of a turbulent 3D pipe bend +permalink: /tutorials/Inc_Turbulent_Bend_Opt/ +written_by: Nijso Beishuizen +for_version: 8.0.1 +revised_by: +revision_date: +revised_version: +solver: INC_RANS +requires: SU2_CFD, FADO +complexity: advanced +follows: Inc_Turbulent_Bend +--- + + + +![bend](../../../tutorials_files/design_features/Inc_Turbulent_Bend/bend_with_FFD.png "Figure (1): The 90 degree bend with the FFD box.") +#### Figure (1): The 90 degree bend with the FFD box. + + +## Goals + +This tutorial closely follows the design optimization setup of the 2D mixing channel, see the tutorial on design optimization for [Species Transport](/tutorials/Species_Transport/). We will use the python framework [FADO](https://github.com/su2code/FADO) to set up the optimization problem. In this tutorial however we will optimize the pressure drop of the 90 degree pipe bend. The CFD results were already discussed previously in [this tutorial](/tutorials/Inc_Turbulent_Bend/). Here we focus on the following aspects: +* Setup of the FFD box for a 3D problem using SU2_DEF. +* Setup of the FADO problem: + - using previous solutions of the CFD and adjoint solver to reduce computing time. + - automatically add the correct keywords for the FFD box (type and degrees of freedom). + - discussion of mesh quality issues and improvement using FFD constraints. +* Presentation of the final setup and results. + +If you have not done so already, please first have a look at the prerequisite tutorials for a better understanding. + +## Prerequisites + +Besides the python library FADO, you will need to compile SU2 with automatic differentiation support. Note that the script provided uses 8 cores, so you need to compile with mpi support as well as enable autodiff: +``` +./meson.py build --optimization=2 -Ddebug=false -Denable-autodiff=true -Dwith-mpi=enabled --prefix=/home/user/Codes/su2_github_develop/su2/ +``` + +## Resources + +You can find the resources for this tutorial in the folder [design/Inc_Turbulent_Bend_Wallfunctions](https://github.com/su2code/Tutorials/tree/master/design/Inc_Turbulent_Bend_Wallfunctions) and the respective subfolders. Note that the setup used in this tutorial uses a pipe bend with a shorter pipe length before and after the bend to reduce the computing time. We have also merged all wall boundaries and all symmetry planes into one wall boundary and one symmetry plane. The gmsh file is provided in the repository so you can create your own pipe bend(s) with it. + +## 1. Basic setup of FFD box and FADO +Usually, designs are created with a CAD tool. These designs are then discretized into a computational mesh for CFD analysis. When optimizing a design, we then only have the discretized mesh available. We could manipulate the mesh nodes directly but this does not lead to very smooth deformations. Instead we modify our mesh using FFD boxes. The nodes of the FFD box are moved according to the design sensitivities, and the mesh nodes inside the FFD box are then smoothly deformed using Bezier curves (default) or B-splines. +Figure (1) above shows the setup that we will be using. + +### Creation of the FFD box +The FFD box can be created and added to the .su2 mesh file using SU2_DEF. The important parameters to add to the configuration file are: + +``` +FFD_DEFINITION= (BOX, \ +-0.06, -0.009, -0.001, \ + 0.209, -0.009, -0.001, \ + 0.209, 0.060, -0.001, \ +-0.06, 0.060, -0.001, \ +-0.06, -0.009 ,0.27, \ + 0.209 ,-0.009, 0.27, \ + 0.209, 0.060, 0.27, \ +-0.06, 0.060, 0.27 ) + +FFD_DEGREE= (5,5,5) +DV_KIND= FFD_SETTING +DV_MARKER= ( wall, symmetry) +DV_PARAM= (1.0) +``` +Our FFD box is 5x5x5 cells, or 6x6x6=216 nodes. With each 3 degrees of freedom for the x-,y- and z-direction, we get a total of 648 d.o.f. The box is slightly larger than our original bend, but most importantly the symmetry plane is completely inside the FFD box. +We run the command + +``` +$ SU2_DEF sudo_0_add_FFD_box.cfg +``` + +We will only use the file *sudo_0_add_FFD_box.cfg* to create the FFD box. This command will create a new .su2 mesh called *mesh_out.su2* that has the definition of the FFD box added to the file. Note that at this stage we need to provide the boundaries inside the FFD box that are allowed to deform using the keyword **DV_MARKER**. The nodes on the boundaries inside the FFD box are part of the information that is added with the FFD box to the mesh_out.su2 mesh file. + + +### Setup the FADO script optimization.py + +Previously the python script *set_ffd_design_var.py* was used in the tutorial on the optimization of the mixing channel. We can however create the correct entries with a couple of simple python commands and add them directly to the main python script that FADO will use, *optimization.py*. We already used *optimization.py* before in the tutorial [Species Transport](/tutorials/Species_Transport/) and we will highlight some changes here. + +For the optimization, we need to modify 2 things in our config file: **DV_KIND** and **DV_PARAM**. The keyword **DV_PARAM** contains N entries of the form *( BOX, 0, 0, 0, 1.0, 0.0, 0.0 );* + Note that the meaning of the entries are *( FFD_BoxTag, i_Ind, j_Ind, k_Ind, x_Disp, y_Disp, z_Disp )* , meaning that after the keyword **BOX**, we get the 3 indices i,j,k of the FFD box, followed by the allowed displacement of that index in the x-,y- and z-direction. Mesh nodes in the symmetry plane only move in the symmetry plane, so symmetry is preserved. + +The list of FFD control points can be created using: +```python +s = "FFD_CONTROL_POINT" +ffd_string = s +for i in range((DX**NDIM)*NDIM - 1): + ffd_string = ffd_string + ", " + s +``` +In the sudo.cfg file we have a placeholder which has the form: +``` +DV_KIND= __FFD_CTRL_PTS__ +``` +And we automatically replace the placeholder **\_\_FFD_CTRL_PTS\_\_** in the fado script during run-time using the command +``` +replace_dv_kind = Parameter([ffd_string], LabelReplacer("__FFD_CTRL_PTS__")) +``` + +For the **DV_PARAM** string, we can use a similar piece of python code: + +```python +dv_param_string="" +for idim in range(NDIM): + xdim = ydim = zdim = "0.0" + if (idim==0): xdim="1.0" + elif (idim==1): ydim="1.0" + elif (idim==2): zdim="1.0" + for k in range(DX): + for j in range(DX): + for i in range(DX): + s = "( BOX, " + str(i) + ", " + str(j) + ", " + str(k) + ", " + xdim + ", " + ydim + ", " + zdim + " );" + dv_param_string += s +``` +which is replaced using: +```python +replace_dv_param =Parameter([dv_param_string], LabelReplacer("__FFD_PARAM__")) +``` + +The only thing we need to take care of now is to define the correct number of FFD nodes DX=6 and the correct dimension of the problem, NDIM=3. You could even get these values from the .su2 mesh file if you want! + +In our config file sudo.cfg, we also use the following setting: +``` +FFD_CONTINUITY= 1ST_DERIVATIVE +``` +Which means that when the sides of the FFD box cuts through the mesh, then at the cut, we want the mesh deformation to stay first order continuous. If you do not do this, you might get sharp jumps in your mesh at the location of the intersection. + +Unfortunately, if we use this unconstrained setup, the mesh deformation will be as shown in Figure (2) below: + +![bend](../../../tutorials_files/design_features/Inc_Turbulent_Bend/mesh_deformation.png "Figure(2): bad mesh for unconstrained FFD box.") +#### Figure (2): The unconstrained deformation leads to mesh cells that have collapsed on the symmetry plane. This leads to convergence issues at the next design iteration. + +As you can see, the row of cells just above the symmetry plane has collapsed onto the symmetry plane, leading to a very bad mesh. In the next design iteration, the simulations do not converge anymore and the optimization procedure stops. So we need some additional constraints on the movement of the nodes of the FFD box. + +### constrained FFD deformation + +The problem here is that the mesh nodes on the symmetry plane are not allowed to move vertically, they only move horizontally outward inside the symmetry plane. Unfortunately, the FFD deformation is such that mesh nodes just above the symmetry plane are moved down, almost on the symmetry plane. To improve the quality of the mesh deformation, we disallow the vertical movement of the FFD box nodes on the nodes in the bottom plane of the FFD box, with j-index 0 and vertical displacement (0,1,0). In Figure (3), the plane with index j=0 is the bottom plane, indicated in yellow. So we remove entries in the **DV_PARAM** list of the form *(BOX, i_Ind, 0, k_Ind, 0,1,0)*. Additionally, we also disallow the vertical movement of the nodes in the plane j=1. +Since we disallow vertical movement in 2x(6x6)=72 nodes in the planes with $$j=0$ and $j=1$$, The total degrees of freedom is then $$648 - 72 = 576$$ d.o.f. + + +![Pressure bend](../../../tutorials_files/design_features/Inc_Turbulent_Bend/opt_iter0_pressure_ffdbox.png "Pressure distribution") +#### Figure (3): The 90 degree bend with the FFD box. + + +The number of design variables needs to be reduced by 2x6x6 and we need to remove 2x6x6 strings from the ffd_string: +```python +nDV = nDV - 2*DX*DX +ffd_string.replace(s+", ","",2*DX*DX) +``` + +Additionally, we need to remove the vertical movement of the nodes in the plane with j=[0,1]: + +```python +jlist = [0,1] +dof = "0.0, 1.0, 0.0" + +for j in jlist: + for k in range(DX): + for i in range(DX): + remove_dof = "( BOX, " + str(i) + ", " + str(j) + ", " + str(k) + ", " + dof + " )" + print("removing ", remove_dof) + dv_param_string = dv_param_string.replace(remove_dof+";", "", 1) + # in case the plane is at the end, the string does not have a semicolon + dv_param_string = dv_param_string.replace(remove_dof, "", 1) +``` + +Another modification to the configuration file that we would like to add is to restart from the solution of the previous design. SU2 will automatically interpolate the solution to the nearest neighbor if the mesh coordinates do not match. Since we might not have an initial solution to start with, we would like to switch from **RESTART= NO** to **RESTART= YES** after the first iteration. + +We can do this with a simple *sed* command that searches and replaces the string in config.cfg in the working directory and then copies the config file back to the base dir: + +```python +restart_yes="sed -i 's/RESTART_SOL= NO/RESTART_SOL= YES/' config.cfg && cp " + configCopy + " ../../" +``` + +The easiest way to perform the update is to simply add it as another ExternalRun process: +```python +update_restart = ExternalRun("UPDATE_RESTART",restart_yes,False) # True means sym links are used for addData +update_restart.addData(configCopy) +``` +This means that we will call this function every design iteration, but only in the first design iteration will it have an effect. + +FADO copies the restart file and all other necessary files from the base directory (where you start your script) into the working directory. Here, the working directory is called *OPTIM*. Inside *OPTIM* we have subdirectories for the different runs, i.e. *OPTIM/DEFORM*, OPTIM/DIRECT*, *OPTIM/ADJOINT* and *OPTIM/DOT*. When all runs in this directory are done, they are archived in *DSN_001*, *DSN_002*, etc. + +So what we will do now is every time after we run the primal solver, we copy the restart file from the working directory *OPTIM/DIRECT* back to the base dir. + +```python +cfd_command = "mpirun -n " + ncores + " " + su2_run + "SU2_CFD " + configMaster + "&& cp restart.csv ../../solution.csv" +``` + +Note that the primal solver by default saves the file restart.csv (keyword: **RESTART_FILENAME= restart**), but reads the file solution.csv (keyword: **SOLUTION_FILENAME= solution**). Note that we save ASCII files here, which is determined by the option **OUTPUT_FILES= RESTART_ASCII**. When reading ASCII restart files, we also need **READ_BINARY_RESTART= NO**. +Please go through the *optimization.py* script now and check the settings. By default, the above procedure to restart from a previous solution is not active. You can activate it now if you wish. + +## 2. FADO optimized result + +We run the optimization script: + +``` +$ python optimization.py +``` + +We have set it to run 5 design iterations. Please note that the optimization was set up for a parallel run. + +The objective values for the 5 design iterations are given in the table below. + + +| iteration | $$\Delta P$$, total| $$\Delta P$$, bend| gain, total| gain, bend | +| --------|--------|--------|--------|--------| +|0|89.0 [Pa]|20.7 [Pa]| -| - | +|1|84.9 [Pa]|16.6 [Pa]|4.6 \%|19.8 \%| +|2|82.1 [Pa]|13.8 [Pa]|7.8 \%|33.3 \%| +|3|81.6 [Pa]|13.3 [Pa]|8.3 \%|35.7 \%| +|4|79.9 [Pa]|11.6 [Pa]|10 \%|44 \%| + + +We see that the global pressure drop between the inlet and the outlet reduces from 89 Pa to 80 Pa, a reduction of more than 10 \%. However, since we are optimizing only the bend, we should subtract the pressure drop of the straight parts and only consider the pressure drop of the bend for a more fair comparison. In paraview, we can integrate the pressure in a 2D slice at the start and end of the bend. The pressure drop of the bend comes to $$\Delta P = 20.7 Pa$$. That means that the reduction of pressure drop in the bend is actually 44 \% ! + +![Optimized bend](../../../tutorials_files/design_features/Inc_Turbulent_Bend/opt_iter5_pressure_ffdbox.png "Optimized bend") +#### Figure (4): Optimized bend after 5 design iterations, with the deformed FFD box. diff --git a/_tutorials/incompressible_flow/Inc_Von_Karman/Inc_Von_Karman.md b/_tutorials/incompressible_flow/Inc_Von_Karman/Inc_Von_Karman.md new file mode 100644 index 00000000..7cf450b8 --- /dev/null +++ b/_tutorials/incompressible_flow/Inc_Von_Karman/Inc_Von_Karman.md @@ -0,0 +1,207 @@ +--- +title: Unsteady von Karman vortex shedding +permalink: /tutorials/Inc_Von_Karman/ +written_by: Nijso Beishuizen +for_version: 8.0.1 +revised_by: +revision_date: +revised_version: +solver: INC_NAVIER_STOKES +requires: SU2_CFD +complexity: intermediate +--- + +![von_karman_street](../../../tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder.gif) + +Figure (1): impression of the von Karman vortex shedding. A high quality mp4 video can be found [here](../../../tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder.mp4) + +## Goals + +In this tutorial we will simulate the laminar but unsteady vortex shedding behind a 2D circular cylinder, also known as von Karman vortex shedding. The unsteady pattern behind the cylinder is known as a von Karman sheet. In this tutorial we will touch upon the following aspects: +- Setting up an unsteady simulation in SU2 +- Choosing appropriate settings for the unsteady simulation +- Inspecting the effects of our choices for the numerical setup +- Creating movies with paraview +- Extracting the shedding frequency with paraview +- visualize the forces on the cylinder with paraview + + +## Resources + +The resources for this tutorial can be found in the [incompressible_flow/Inc_Von_Karman](https://github.com/su2code/Tutorials/tree/master/incompressible_flow/Inc_Von_Karman_Cylinder) directory in the [tutorial repository](https://github.com/su2code/Tutorials). You will need the configuration file ([unsteady_incomp_cylinder.cfg](https://github.com/su2code/Tutorials/tree/master/incompressible_flow/Inc_Von_Karman_Cylinder/unsteady_incomp_cylinder.cfg)) and the mesh file ([cylinder_wake.su2](https://github.com/su2code/Tutorials/tree/master/incompressible_flow/Inc_Von_Karman_cylinder/cylinder_wake.su2)). Additionally, the Gmsh geometry is also provided so you can recreate the mesh yourself: [cylinder_wake.geo](https://github.com/su2code/Tutorials/tree/master/incompressible_flow/Inc_Von_Karman_Cylinder/cylinder_wake.geo). + + +### Background + +When the Reynolds number $$Re=\rho \cdot V \cdot D / \mu$$ is low (Re < 40), the flow around a circular cylinder is laminar and steady. At around Re=49, the flow becomes unsteady and a periodic shedding of vortices forms in the wake of the cylinder, known as vortex shedding. The frequency of this vortex shedding is usually expressed in terms of the Strouhal number $$St= f \ cdot D / U_{\infty}$$, with f the shedding frequency, D the diameter of the cylinder and U the far-field velocity. Important experimental work can be found in the paper of Williamson, *Vortex Dynamics in the Cylinder Wake*, Annual Review of Fluid Mechanics (1996) [doi](https://doi.org/10.1146/annurev.fl.28.010196.002401). After around Re=180, a second frequency is observed experimentally, in the longitudinal direction. This frequency can only be observed in a 3D simulation. The increase in number of frequencies continues until after around Re=1000 the flow is considered fully turbulent. + +### Problem Setup + +The configuration is a circular cylinder of 5 mm surrounded by a far field at $$L = 30 D$$ and a rectangular wake region of $$X = 150 D$$. The far-field velocity is $$U_{\infty} = 0.12 m/s$$. With a viscosity of $$\mu=1.0 \cdot 10^{-5}$$ and a density of $$\rho = 1 kg/m3$$, the Reynolds number is $$Re=120$$. + +![von_karman_mesh](../../../tutorials_files/incompressible_flow/Inc_Von_Karman/images/mesh.png) +Figure 2: Computational domain for the von Karman vortex shedding. + +The solution is initialized with uniform flow conditions. After an initial phase of around 4 seconds, the flow becomes periodic with a specific frequency depending on the Reynolds number. + + +### Mesh Description + +The mesh consists of a structured mesh with 70k cells and 75k points. The mesh was created using Gmsh and the configuration file to create the mesh can be found here: [cylinder_wake.geo](https://github.com/su2code/Tutorials/tree/master/incompressible_flow/Inc_Von_Karman/cylinder_wake.geo). The only thing you need to do to create a mesh from the geometry is start Gmsh, and then load the .geo file. You will then see the geometry in the Gmsh visualization window. If you click on *Mesh->2D* the 2D mesh will be generated. You can then export the mesh as a .su2 file by choosing *File->Export*. The mesh will automatically be saved in su2 format when the filename has the extension .su2. In general, you should not choose *save all elements* because this will also save additional points that were used to construct the geometry but are not part of the final mesh, like for example the center of a circle. + + +### Configuration File Options + +Several of the key configuration file options for this simulation are highlighted here. First, we activate the turbulence model: + +``` +SOLVER= INC_NAVIER_STOKES +INC_NONDIM= DIMENSIONAL +KIND_TURB_MODEL= NONE +``` + +We use the incompressible NAVIER-STOKES without turbulence model, since the flow is unsteady but laminar. + +``` +TIME_DOMAIN = YES +TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER +TIME_STEP= 0.005 +MAX_TIME=25.00 +TIME_ITER=2500 +INNER_ITER= 20 +``` + +The **TIME_DOMAIN** keyword activates transient simulations. We use a second order timestep, with a step size of $$\Delta t = 0.01 s$$. The maximum time is set to $$t_{max} = 25 s$$ or 2500 iterations, whichever is reached first. We use 20 inner iterations per timestep, which sufficiently converges the residuals. This is important to check, insufficient convergence per time step quickly leads to large errors. + +``` +INC_VELOCITY_INIT= ( 0.12, 0.0, 0.0 ) +MARKER_HEATFLUX= ( cylinder, 0.0 ) +MARKER_FAR= ( farfield_in,farfield_side,farfield_out ) +MARKER_MONITORING= ( cylinder ) +``` + +There is no heat transfer, so we set zero heatflux boundary conditions on the walls and we impose a far-field velocity of 0.12 m/s. The **MARKER_MONITORING** keyword is used to monitor forces at the cylinder surface. We are interested in the lift and drag coefficient. + +``` +% monitoring points in the cylinder wake +CUSTOM_OUTPUTS= 'velocity : Macro{sqrt(pow(VELOCITY_X, 2) + pow(VELOCITY_Y, 2) )};\ + probe1 : Probe{$velocity}[0.15, 0.0]' +HISTORY_OUTPUT= (TIME_ITER, CUR_TIME, RMS_RES, LIST, DRAG, LIFT, SURFACE_STATIC_PRESSURE, CUSTOM) +``` + +We define a monitoring point called *probe1* at location $$(x,y)=(0.15, 0.0)$$ downstream of the cylinder and at the centerline. In this probe, we monitor the velocity magnitude $velocity$, that we compute from the existing velocity components. We write the probe information to the history output using the keyword **CUSTOM**. + + +### Running SU2 + +If possible, always use a parallel setup to reduce computational time (wall clock time). Run the SU2_CFD executable in parallel using MPI and 4 nodes by entering: + + $ mpirun -n 4 SU2_CFD unsteady_incomp_cylinder.cfg + + + +### Results + +![particles_streamlines](../../../tutorials_files/incompressible_flow/Inc_Von_Karman/images/particles.png) +Figure (3): Developed von Karman street with injected particles. + +The figure above shows the von Karman vortices. A movie of the von Karman Vortex shedding as shown in this tutorial can be made using the paraview statefile here: [statefile_with_particles.pvsm](https://github.com/su2code/Tutorials/tree/master/incompressible_flow/Inc_Von_Karman_Cylinder/statefile_with_particles.pvsm). The colored vorticity gradients are created in the wake immediately behind the cylinder, and transported downstream. In paraview, we inject particles and let them be transported downstream over streamlines. We can clearly see that the particles are not mixed homogeneously but are being transported in batches of red and green particles moving with their corresponding vortices. + + +![von_karman_street_arrow](../../../tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder_arrow.gif) + +A link to a high quality video can be found [here](https://github.com/su2code/su2code.github.io/tree/master/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder_arrow.mp4). + +The forces on the cylinder can also be visualized, see the above movie. We see the unsteady velocity magnitude, and on the cylinder surface we visualize the local force as the local pressure normal to the surface. For this we need to extract the surface data only and then compute the surface normal using a programmable filter. We can then visualize the glyph distribution. The total lift force is the integrated vertical component, whose size can be visualize as a colored arrow of a size proportional to the lift force using another programmable filter. Below the image we visualize the lift force in time using the PlotSelectionOverTime feature. + +![validation](../../../tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman.png) + +Figure (4): Strouhal number, comparison with experimental data from Williamson (2006). + +The Strouhal number can be computed from the lift coefficient that was stored in the history.csv file. A simple python file gives us the frequency from the FFT: + +```python +#%matplotlib inline +#import matplotlib.pyplot as plt +import pandas as pd +import numpy as np +import scipy.fftpack +import csv + +# Return value belonging to key in config.cfg +# (splits key= value for you) +def find_config_key_value(filename,config_key): + with open(filename, "r") as file: + for line in file: + line = line.split('=') + if line[0] == config_key: + print(line[-1].strip() ) + return(line[-1].strip()) + raise ValueError('key not found:',config_key) + +# diameter of cylinder +D = 0.01 + +# read history, use comma as separator, with zero or more spaces before or after +df = pd.read_csv('history.csv', sep='\s*,\s*') + +T = float(find_config_key_value('unsteady_incomp_cylinder.cfg','TIME_STEP')) +U = find_config_key_value('unsteady_incomp_cylinder.cfg','INC_VELOCITY_INIT') +N = len(df.index) +print('timestep=',T) +print('samplepoints=',N) +print('velocity=',U) +U = float(U.replace('(','').replace(')','').split(',')[0].strip()) +print('velocity=',U) + +# assign data +x = df['"Cur_Time"'] +y = df['"CL"'] + +# compute DFT with optimized FFT +xf = np.linspace(0.0, 1.0/(2.0*T), N//2) +yf = np.fft.fft(y) +yf2 = 2.0/N * np.abs(yf[:N//2]) + +#fig, ax = plt.subplots() +#plt.xlim(0,5.0) +#ax.plot(xf, yf2 ) +#plt.show() + +print("index of max = ",np.argmax(yf2)) +freq = xf[np.argmax(yf2)] +print("frequency of max = ",freq) +St = freq * D / U +print("strouhal number = ",St) +``` + + +When we compare the Strouhal number with the experimental data from Williamson, we see in Figure 4 that the frequency is underpredicted. We will vary some numerical settings to investigate if we can improve the prediction of the Strouhal number. + + + +### Numerical variations + +![validation](../../../tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman_variation.png) + +Figure (5): Comparison of different numerical settings + +In Figure 5, we see the effect of different numerical settings on the prediction of the Strouhal number. The second order scheme predicts a Strouhal number of $$St = 0.1734$$, slightly over predicting the experimental value of $$St_{exp} = 0.170$$. Note that our predictions of the Strouhal frequency depends on the number of samples and sampling rate that we provide to the FFT. We took 2500 timesteps of 0.01 s which contains enough cycles for an accurate frequency prediction using an fft. +When switching from second order in time to first order, the Strouhal number is under predicted by 6 \% compared to the experimental value. Also, when increasing the time step from 0.01 s to 0.02 seconds, the St decreases by 2 \%. When increasing the time step even further to $$ \Delta t = 0.04 s%%, St is under predicted by 8 \%. The period of the dimensional frequency is $$f \approx 0.5 s$$, so with a timestep of 0.01 s we have 50 time steps per period, we have 25 time steps when $$\Delta t = 0.02 s$$, and only 12 time steps when $$\Delta t = 0.04 s$$. It is clear that 12 time steps per period is not sufficient. + +It is also known that the size of the computational domain influences the results, so we reduce the domain by half, $$L = 15 D$$ and $$X = 75 D$$. The Strouhal then increases to $$St = 0.1768$$, an increase of 2 \%. It seems that a far-field that is 15D away from the cylinder is sufficient. + +As a final test, the testcase can be executed for varying Reynolds numbers, ranging from Re=60 to Re=180, giving the result in Figure (6). + +![validation](../../../tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman_Re.png) + +Figure (6): Comparison of numerical Strouhal numbers with experiments for varying Reynolds numbers. + +We get a pretty good agreement compared to the experimentally measured values. + +### Final notes + +The paraview statefile to create the movie can be found here: [statefile_with_particles.pvsm](https://github.com/su2code/Tutorials/blob/master/incompressible_flow/Inc_Von_Karman/statefile_with_particles.pvsm) +and here: +[statefile_movablearrow_timeseries.pvsm](https://github.com/su2code/Tutorials/blob/master/incompressible_flow/Inc_Von_Karman/statefile_movablearrow_timeseries.pvsm) +Note that you have to select your own, local files when you load the statefile. \ No newline at end of file diff --git a/_tutorials/index.md b/_tutorials/index.md index 66bb006f..778e8d80 100644 --- a/_tutorials/index.md +++ b/_tutorials/index.md @@ -50,6 +50,8 @@ Simulation of unsteady, external, viscous flow around an airfoil. Perform uncertainty quantification of errors arising due to assumptions inherent in turbulence models. * [Non-ideal compressible flow in a supersonic nozzle](/tutorials/NICFD_nozzle/) Simulation of compressible flow in a nozzle using non-ideal thermodynamic models. +* [Turbomachinery: Aachen Turbine stage with mixing plane](/tutorials/Aachen_Turbine/) +Simulation of compressible flow of the Aachen turbine demonstrating turbomachinery application. #### Incompressible Flow @@ -71,6 +73,8 @@ Simulation of internal, turbulent, incompressible flow in a unit cell of a 2D pi Simulation of internal, turbulent, incompressible flow through a mixing channel. * [Species Transport Composition Dependent Model](/tutorials/Inc_Species_Transport_Composition_Dependent_Model/) Simulation of internal, turbulent, 3D incompressible flow through a Kenics static mixer. +* [Vortex Shedding](/tutorials/Inc_Von_Karman/) +Simulation of unsteady laminar vortex shedding behind a circular cylinder. * [Turbulent Bend](/tutorials/Inc_Turbulent_Bend/) Simulation of turbulent flow in a 90 degree pipe bend using wall functions. @@ -114,3 +118,5 @@ Perform an optimal shape design with multiple objectives and a penalty function Shape optimization of an 2D airfoil in unsteady flow conditions. * [Species Transport](/tutorials/Species_Transport/) Optimization of internal, turbulent, incompressible flow through a mixing channel. +* [Optimizing pressure drop of a pipe bend](/tutorials/Inc_Turbulent_Bend_Opt/) +Optimization of the pressure drop (also known as head loss), of a pipe bend. diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/2D_MidSpan_RelMach.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/2D_MidSpan_RelMach.png new file mode 100644 index 00000000..538c5e14 Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/2D_MidSpan_RelMach.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/FullAnnulus.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/FullAnnulus.png new file mode 100644 index 00000000..30f62553 Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/FullAnnulus.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/FullAnnulus_MachContour.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/FullAnnulus_MachContour.png new file mode 100644 index 00000000..f154cb92 Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/FullAnnulus_MachContour.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/Full_View_MachContour.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/Full_View_MachContour.png new file mode 100644 index 00000000..a7c2f262 Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/Full_View_MachContour.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Giles.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Giles.png new file mode 100644 index 00000000..0ac98b4e Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Giles.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Giles_MixingPlane.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Giles_MixingPlane.png new file mode 100644 index 00000000..1e5da5b6 Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Giles_MixingPlane.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_NoSlip.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_NoSlip.png new file mode 100644 index 00000000..42801fb8 Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_NoSlip.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Periodics.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Periodics.png new file mode 100644 index 00000000..2cbb1e4f Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_Adiabatic_Periodics.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_shroud.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_shroud.png new file mode 100644 index 00000000..8033b78b Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/SingleMesh_shroud.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/SinglePassage.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/SinglePassage.png new file mode 100644 index 00000000..1abf94cd Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/SinglePassage.png differ diff --git a/tutorials_files/compressible_flow/Aachen_Turbine/images/SinglePassageMesh.png b/tutorials_files/compressible_flow/Aachen_Turbine/images/SinglePassageMesh.png new file mode 100644 index 00000000..e4538d03 Binary files /dev/null and b/tutorials_files/compressible_flow/Aachen_Turbine/images/SinglePassageMesh.png differ diff --git a/tutorials_files/design_features/Inc_Turbulent_Bend/bend_with_FFD.png b/tutorials_files/design_features/Inc_Turbulent_Bend/bend_with_FFD.png new file mode 100644 index 00000000..93c5f856 Binary files /dev/null and b/tutorials_files/design_features/Inc_Turbulent_Bend/bend_with_FFD.png differ diff --git a/tutorials_files/design_features/Inc_Turbulent_Bend/mesh_deformation.png b/tutorials_files/design_features/Inc_Turbulent_Bend/mesh_deformation.png new file mode 100644 index 00000000..daee5aca Binary files /dev/null and b/tutorials_files/design_features/Inc_Turbulent_Bend/mesh_deformation.png differ diff --git a/tutorials_files/design_features/Inc_Turbulent_Bend/opt_iter0_pressure_ffdbox.png b/tutorials_files/design_features/Inc_Turbulent_Bend/opt_iter0_pressure_ffdbox.png new file mode 100644 index 00000000..498c766e Binary files /dev/null and b/tutorials_files/design_features/Inc_Turbulent_Bend/opt_iter0_pressure_ffdbox.png differ diff --git a/tutorials_files/design_features/Inc_Turbulent_Bend/opt_iter5_pressure_ffdbox.png b/tutorials_files/design_features/Inc_Turbulent_Bend/opt_iter5_pressure_ffdbox.png new file mode 100644 index 00000000..f0bffb98 Binary files /dev/null and b/tutorials_files/design_features/Inc_Turbulent_Bend/opt_iter5_pressure_ffdbox.png differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/mesh.png b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/mesh.png new file mode 100644 index 00000000..7f35d19c Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/mesh.png differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/particles.png b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/particles.png new file mode 100644 index 00000000..c56f57a0 Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/particles.png differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman.png b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman.png new file mode 100644 index 00000000..e4f861bf Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman.png differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman_Re.png b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman_Re.png new file mode 100644 index 00000000..631dded9 Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman_Re.png differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman_variation.png b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman_variation.png new file mode 100644 index 00000000..8c1cef0a Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/strouhal_cylinder_karman_variation.png differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder.gif b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder.gif new file mode 100644 index 00000000..0a90f4d0 Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder.gif differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder.mp4 b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder.mp4 new file mode 100644 index 00000000..0fcfc6bf Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder.mp4 differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder_arrow.gif b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder_arrow.gif new file mode 100644 index 00000000..aa292f80 Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder_arrow.gif differ diff --git a/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder_arrow.mp4 b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder_arrow.mp4 new file mode 100644 index 00000000..15f39cec Binary files /dev/null and b/tutorials_files/incompressible_flow/Inc_Von_Karman/images/unsteady_cylinder_arrow.mp4 differ