post processing in CFD

 

post

post

Table of Contents

1 Post-processing

 

1.1 Reverse flow

It is common to encounter the regions of reversed flow at the initial stages of the simulation – This is normal.

However, if the “reversed flow” warnings do not disappear as the simulation progresses, then one needs to address the issue and move the outlet boundaries to a location where the inflow is no longer encountered.

Normally, it has to do with the outlet boundary condition.

if you were to use the "pressure outlet" boundary condition at the outlet, your outlet must be set far away from the object of interest.

In Fluent, I normally use the " outflow" boundary condition at the outlet and it does not give me reverse flow.

1.1.1 reasons

It is virtually impossible to prescribe correct values for varying turbulence characteristics, temperature and species concentrations in those cells where the backflow occurs.

  • The number of such “backflow” cells might vary from one iteration to another,

and any discontinuity in calculated turbulence, species and temperature values for the neighbouring “outflow” and “inflow” cells is very likely to have an adverse effect on convergence.

  • The accuracy of the solution can also be affected, especially if the area of interest is located close to the outlet with the “backflow”.

reversed_flow.png

1.1.2 solutions

  • if vortex formation/recirculation of flow near an outlet boundary, 1.1.2.

> check the mesh quality, improve it.

> use higher order scheme

> At last reduce the relaxation factor if necessary.

Calculation steps:

  • Set up the model and calculate a solution on the initial mesh
  • Please write the interpolation file first.
  • Close the existing Fluent session.
  • Open a new Fluent session and read the case file only.
  • Perform 1.1.2
  • Check the new zones and ensure that the prescribed outlet boundary conditions are set and named consistently.
  • Initialise the flow.
  • Read in the interpolation file.
  • Continue simulation.

  • The original data file cannot be read onto the new mesh.
  • The new case file need to be initialised.
  • Old data can still be recycled and used again by interpolation:
 file interpolate read-data original_ip_file_name.ip

Mesh/modify-zones/extrude-face-zone-delta 
#specify the face zone id/name i.e. inlet boundary name or zone ID 
distance delta 1 0.2[1]
#first extrusion layer edge length
distance delta 2 0.25 
distance delta 3 0.3 
distance delta 4 
#press “Enter” once a sufficient outlet section has been built[2]***.
  1. Reference

    • 5.9.8. Extruding Face Zones, user guide • File “prevent backflow at outlet.fss”

  2. Note:
    • This text command is not available in the parallel version of ANSYS Fluent.

    extrusion is not possible from boundary face zones that have hanging nodes.

    • Extruding face zones is not allowed on polygonal face zones.
    • Extruding face zones is only allowed in the 3D version of ANSYS Fluent.
  3. outlet Extrusion error

    Problem description:

    • mesh/modify-zones/extrude-face-zone-delta
    • Face zone id/name [] 27 #outlet face zone
    • Distance delta(1) [()] 0.1
    • Distance delta(2) [()]
    • Extrude face zone? [yes]
    • Warning: the base thread has multiple adjacent threads!
    • Please separate the extruded side threads manually, and set the boundary condition accordingly.
    • Error: received a fatal signal (Segmentation fault).

1.2 variable definitions– ch34, field function definitions, user guide, fluent

1.3 Residuals

residuals should below 1e-4

  • If you are having convergence difficulties, plot the residual value fields (for example, using contour plots) to determine where the high residual values are located

> 29.15.1.11. Postprocessing Residual Values

1.3.1 access to the residual values for each cell in my flow domain?

save the residuals for post processing (in TUI solver/set/expert ) and then to scrutinize the areas where the maximum values of the unscaled residuals are happening and to find out whether they are related to a local issue such as poorly meshed regions (poor metrics or low resolution), domain interface, local abrupt physic change (shock), etc. or to a global effect.

Answer:

These values can be activated by the following TUI command:

/solve/set/expert

note: The residual values that are made available for postprocessing are the unscaled values.

One option listed here is: "Save cell residuals for post-processing? [no]" This option has to be set to "yes" or "y".

example log: ** /solve/set/expert Save cell residuals for post-processing? [no] yes Keep temporary solver memory from being freed? [no] no Allow selection of all applicable discretization schemes? [no] no *

Afterwards you have to do one dummy iteration so that the specific data fields are created.

Now the residuals of all solved equations can be accessed and postprocessed within Fluent, see Screenshot1[2041045_Screenshot1.png]. > contours/ > contours of residuals

They can also be used for convergence monitors.

If you want to store the residuals in your data file, you have to adjust the data file quantities (File > Data File Quantities) see Screenshot2.[2041045_Screenshot2.png].

1.3.2 access the local residuals through UDF

#2052972 Product Family: Fluid Dynamics Last Updated: May 08 2018

Answer:

It might be useful to access the residuals of variables at the cell level.

In the TUI of Fluent, type: solve/set/expert and answer 'yes' to 'Save cell residuals for post-processing?'

You can now access the residuals in the post-processing.

It is also possible to access the residual values in a UDF. An example is provided in attachment.

In this example, a first DEFINEONDEMAND stores the residual values of Pressure (= Mass Imbalance) X,Y and Z velocity, and the temperature in 5 UDMs.

A second DEFINEONDEMAND use a general way to access any post processing variable with CPOSTVAR. In this second macro, "temperature-residual" can be replaced by the name you would enter in the TUI if you were to use TUI commands to produce e.g. a contour plot…

Disclaimer Note that UDF is not covered by the ANSYS support contract. The UDF in this solution is provided as it is without support. It was used by ANSYS internally with Fluent 190.. It might work with previous and future versions but there is no guarantee.

  1. localresidual.c

    /*Note that UDF is not covered by the ANSYS support contract. The UDF in this solution is provided as it is without support. It was used by ANSYS internally with Fluent 190.. It might work with previous and future versions but there is no guarantee.*/

    #include "udf.h" DEFINEONDEMAND(reportresiduals) { cellt c; Thread *ct; Domain *d=GetDomain(ROOTDOMAINID);

    if(RPGetBoolean("save-cell-residuals?")) { Message("\nSaving Residuals in UDMs…\n"); threadloopc(ct,d) { begincloop(c,ct) { CUDMI(c,ct,0) = CSTORAGER(c,ct, SVPR); CUDMI(c,ct,1) = CSTORAGER(c, ct, SVUR); CUDMI(c,ct,2) = CSTORAGER(c, ct, SVVR); CUDMI(c,ct,3) = CSTORAGER(c, ct, SVWR); CUDMI(c,ct,4) = CSTORAGER(c, ct, SVTR); } endcloop(c,ct) } } else { Message("\nResiduals are not set to be saved.\n"); } } DEFINEONDEMAND(alternative) { Domain *domain = GetDomain(1); #if !RPHOST Thread *t; cellt c; #endif CellFunctionValues(domain, "temperature-residual"); #if !RPHOST if (RPGetBoolean("save-cell-residuals?")) { Message0("\nSaving Residuals in UDMs…\n"); threadloopc(t, domain) { begincloop(c, t) { CUDMI(c, t, 4) = CPOSTVAR(c, t); } endcloop(c, t) } } else { Message0("\nResiduals are not set to be saved.\n"); } #endif }

1.3.3 How to interpret residuals in Explicit Density-Based turbulence analysis?

When using the Explicit solution method, the residual values may not indicate convergence in the usual way, because they are defined as time rate of change of the conserved variable . Therefore, in this case, as the flow field changes with time, the value of the residuals may potentially increase or stay relatively flat. However, this does not necessarily mean that the solution is divergent.

If you are having convergence difficulties, it is often useful to plot the residual value fields (for example, using contour plots) to determine where the high residual values are located

It should be noted that the residual values that are displayed for the Implicit method of Density-based analysis will indicate the convergence of the calculation at each time step

1.4 Checklists

  • residual <1e-4
    • if residual is high, plot the residual and find out the locations of high residuals
  • mass imbalance < 0.5%
  • y plus
    • y plus are different at different iteration stages
    • for wall function, (30-200)
    • enhanced wall, sst k-omega, y plus < 2
  • number of nodes inside BL
    • method 1: plot eddy viscosity (BL thickness) vs mesh at a plane normal to the wall boundary
    • method 2: plot velocity at a plane normal to the wall boundary
  • Cp
  • blade surface Pressure

1.5 Animation

  • Methods
    • Using Animations Sequences available under Calculation Activities > Solution Animations
    • Saving images periodically using TUI commands

(a movie can then be created using a third party software like Windows Movie Maker)

  • Solution Animation File

If you selected Metafile or PPM Image under Storage Type in the Animation Sequence panel, then FLUENT will save the solution animation file for you automatically. It will be saved in the specified Storage Directory, and its name will be the Name you specified for the sequence, with a .cxa extension (e.g., pressure-contour.cxa). In addition to the .cxa file, FLUENT will also save a metafile with a .hmf extension for each frame (e.g., pressure-contour0002.hmf). The .cxa file contains a list of the associated .hmf files, and tells FLUENT the order in which to display them. From http://jullio.pe.kr/fluent6.1/help/html/ug/node875.htm

1.5.1 mesh deformation

  • use journal file to save mesh motion for sliding mesh model

k

1.5.2 tutorials

streamlines animation using keyframe animation

1000220streamlineanimation.pdf 2039574vofanimation.pdf

  1. Floating Objects

    FluentFloatingObjects17.0 >> file://c:/akmkemin/Backup/ANSYS/Tutorial/Free_Surface/Floating_Objects_17.0/Fluent_Floating_Objects_17.0

    1. imagessetup.jou

      display set contours clip-to-range? no display set contours filled-contours? yes display set contours global-range? no display set contours n-contour 100 display set contours node-values? yes display set contours render-mesh? yes display set filled-mesh? yes display set render-mesh? yes display set lights lights-on? yes display set lights headlight-on? yes display set lights lighting-interpolation gouraud display set picture color-mode color display set picture driver png display set picture invert-background? yes display set picture x-resolution 1920 display set picture y-resolution 1080 display set windows scale format "%0.1f" display set windows scale alignment left display set windows scale clear? yes display set windows scale border? yes display set windows text visible? yes display set windows text date? no display set windows logo? no display set windows axes visible? no display set colors color-by-type? no

    2. imagessave.jou

      display set-window 5 display set contours surfaces _vof () display set mesh-surfaces wallsdof () display views camera projection perspective display set titles left-top "Free surface level View 1, top front" display contour mixture y-coordinate -1 1 display views restore-view vof-0 display save-picture ./zvof-0-%7t.png y display set titles left-top "Free surface level View 2, top front close" display contour mixture y-coordinate -1 1 display views restore-view vof-1 display save-picture ./zvof-1-%7t.png y display set titles left-top "Free surface level View 3, top back" display contour mixture y-coordinate -1 1 display views restore-view vof-2 display save-picture ./zvof-2-%7t.png y display set titles left-top "Free surface level View 4, top back" display contour mixture y-coordinate -1 1 display views restore-view vof-3 display save-picture ./zvof-3-%7t.png y display set titles left-top "Free surface level View 5, bottom back" display contour mixture y-coordinate -1 1 display views restore-view vof-4 display save-picture ./zvof-4-%7t.png y display set titles left-top "Free surface level View 6, bottom front close" display contour mixture y-coordinate -1 1 display views restore-view vof-5 display save-picture ./zvof-5-%7t.png y display set titles left-top "Free surface level View 7, top front/side far" display contour mixture y-coordinate -1 1 display views restore-view vof-6 display save-picture ./zvof-6-%7t.png y display views restore-view vof-7 display set titles left-top "Free surface level View 8, top front/side reversed" display contour mixture y-coordinate -1 1 display save-picture ./zvof-7-%7t.png y display views restore-view vof-8 display save-picture ./zvof-8-%7t.png y display set titles left-top "Free surface level View 9, top frot/side reversed close" display contour mixture y-coordinate -1 1 display set titles left-top "Mesh from side" display surface-mesh wallsdof _x=0 () display views restore-view mesh-x0 display views camera projection orthographic display update-scene select-geometry yes _x=0 no no display update-scene display yes no no yes yes yes yes yes no 83. 83. 83. 0. display save-picture ./zmesh-x0-%7t.png y display views restore-view mesh-x1 display views camera projection orthographic display save-picture ./zmesh-x1-%7t.png y display set titles left-top "Mesh from top view" display surface-mesh wallsdof _y=0 () display views restore-view mesh-y0 display views camera projection orthographic display update-scene select-geometry yes _y=0 no no display update-scene display yes no no yes yes yes yes yes no 83. 83. 83. 0. display save-picture ./zmesh-y0-%7t.png y display views restore-view mesh-y1 display views camera projection orthographic display save-picture ./zmesh-y1-%7t.png y display close-window 5

1.5.3 Animation/residual/monitor on Linux Cluster (TUI)

The cluster cannot plot any figures. If you do, errors transpire and you’ll have to delete the job (see Section 3.3). This includes animations, residuals and monitors.

If you want to make animations, write the data files and use TecPlot on your PC.

Your output file is where you will monitor the convergence of your simulation. From: http://web2.clarkson.edu/projects/cluster/FLUENT_cluster.pdf From http://www.cfd-online.com/Forums/fluent/90579-creating-animation-through-tui.html Here is a sample journal file for a 2D simulation in batch mode.

  1. journal.jou

    ;read case and data file/read-case single-48cm-2D-batchtest.cas ; ;initialize domain with air /solve/initialize/initialize-flow ; ;patch water to raceway depth of 3.5 cm /adapt/mark-inout-rectangle yes no 0 0.62875 -0.48 0.035 /solve/patch water () (0) mp 1 ; ;set up image output /display/set/contours/filled-contours yes /display/set/picture/driver png /display/set/picture/landscape yes /display/set/picture/x-resolution 960 /display/set/picture/y-resolution 720 /display/set/picture/color-mode color /views/restore-view front ; ;print front view of phases and velocity magnitude at t=0 /display/contour air vof 0 1 /display/save-picture airvof%t.png /display/contour mixture velocity-magnitude 0 0.5 /display/save-picture velmag%t.png ; ;set up display commands to print front view of phases and velocity magnitude every 20 time steps /solve/execute-commands/add-edit command-2 20 "time-step" "/display/contour air vof 0 1" /solve/execute-commands/add-edit command-3 20 "time-step" "/display/save-picture airvof%t.png" /solve/execute-commands/add-edit command-4 20 "time-step" "/display/contour mixture velocity-magnitude 0 0.5" /solve/execute-commands/add-edit command-5 20 "time-step" "/display/save-picture velmag%t.png" ; ;set up auto-save /file/auto-save data-frequency 500 /file/auto-save append-file-name-with time-step 6 ; ;iterate over 5000 time steps solve/set/time-step 0.001 solve/dual-time-iterate 5000 50 file/write-data single-48cm-2D-batchtest.dat

1.5.4 Fluent Animation video examples

file://c:/akmkemin/Backup/ANSYS/Post/Animation

author: cynthia ecklor

animations for A fan-stage >> C:\akmkemin\Backup\ANSYS\Post\Animation\animationFluent Fly in Flight, Vortex Shedding around a circular cylinder with LES Floating Box 3D. Wright Flyer 1 (Created by Advantage CFD using Ensight from CEI) and Wright Flyer 2. Airplane Wing 2D. Airplane Wing 3D. Spinning Rugby Ball, Train moving through a tunnel, and Passing cars.

Chapter 7 (3454.0K) In this zip file you will find animations for Rocket Launch.

Chapter16 (1778.0K) In this zip file you will find animations for Pipe Flow.

Chapter 20 (10882.0K) In this zip file you will find animations for Toy Airplane, Paper Airplanes 20 (Courtesy of Dartmouth College), and Oscillating NACA Wing.

All Others (51841.0K) In this zip file you will find animations for Rocket Separation, IC Engine, Rotary Engine 2D, and Rotary Compressor.

1.6 Questions

  • Interest of Quantites?
  • Am I using the correct turbulence model for the types of results I am looking for?
  • Do I have appropriate Y + value and a sufficient number of inflation layers?
  • Am I using the right wall function for my problem?

1.7 y plus value

  • set up reference values ( ρ, μ, velocity ) before show y plus contour

Reynolds number uses the reference length, density, and viscosity.

1.7.1 How is y plus calucated in Fluent?

Laminar Flow :

  • velocity gradient >> wall shear stress >> frictional velocity, uτ >> y plus
  • set reference values: density, viscosity

1.7.2 How to?

For a 3D Case

  1. Click Display-> Contour-> Turbulance-> Y Plus of the surface of interest.
  2. Turn off Node Value display.
  3. You can now check if the minimum and maximum Y+ forms an appropriate range

For a 2D Case

  1. Click Display-> Contour-> Turbulence-> Y Plus. You do not need to select any Surfaces
  2. Turn off Node Value display
  3. You can now check the value of Y+ (note to ignore the value of Y+ except for cells adjacent to the wall).

In a fine mesh, a recommended value of 1 for Y+ in a fine mesh, and 30 to 500 for a wall function.

1.7.3 Y plusfor different turbulence models

  • k -e model, (12, 300)
  • SST k -omega model, y+ < 2

For those who are interested, 11.06 is the y+ value where the linear velocity profile in the sublayer intersects with the logarithmic velocity profile in the log layer.

For the k-e model you should use y+ < 300. If y+ is below about 11 for the k-e model then it still works fine, but it doesn't make use of the fine near wall mesh, so it's a waste of mesh. Essentially the k-e model always uses the wall function approach. The wall function approach is not valid below y+ ~= 11/15, so we just ignore the mesh below y+ ~= 11.

If you want accurate boundary layer predictions, such as separation prediction, then you should use the SST model with y+ < 2.

In this case the SST model will switch to a low-Re formulation near the wall rather than the wall function formulation. SST will still work fine with 11 < y+ < 300, but the results will be fairly similar to the k-e model since it will be using wall functions. For 2 < y+ < 11 it will be a blend between wall functions and low-Re formulations.

In some situations, such as accurate boundary layer heat transfer predictions or when using the transition model, an even lower y+ of about 1 is recommended with the SST model.

1.8 Velocity profile along normal direction of airfoil surface

  • 10-15 nodes in the BL
  • Estimate BL thickness via plotting velocity normal to airfoil surface

CFD-Post does not know a priori the curvature of your airfoil so it cannot decide what is the normal direction at every point.

If you know the normal direction (because you have the equation of the airfoil or you have calculated the curvature using a custom function) then you can create the line in that direction and plot the velocity over that line. If you don't know the direction, you will need to use Pearl programming to find out neighboring faces to build up the curvature and normal (as you may imagine this requires programming skills).

Velocity Profile of a flat Plate From https://confluence.cornell.edu/display/SIMULATION/Flat+Plate+Boundary+Layer+-+Numerical+Results

1.9 Boundary Layer Characteristics

Flat plate Boundary Layer >> simcafe tutorial

1.9.1 the boundary layer thickness on a 3D external aerodynamics in post-processing

less accurate Method 1: Total Pressure loss (BL "Health" indicator) >> 2037906 boundarylayerthicknessSolution2037906 +software: EnSight can extract BL thickness accurate Method 2: local velocity profile plot velocity profile normal to the wall

1.9.2 References

2043492 BL profiles, thickness, separation, etc. https://www.ceisoftware.com/wp-content/uploads/2012/05/boundary_layer.pdf extraction of boundarylayer characteristics from Fluent results1.pdf

1.9.3 BL thickness

  • double size of max eddy visocosity ratio +max turbulent viscosity in the middle of the boundary layer
  • based on BL thickness of flat plate

Method 1: If you have the velocity field, just calculate the vorticity, omega, it decays exponentially near the boundary-layer edge because outside you have irrotational potential flow.

Then you may state that delta is where the vorticity has decayed to a small fraction of the maximum value that is closer to the wall.

Method 2 du/dy, and you integrate a new u*=Integral(omega times dy)

normal to the local wall, and that u* reaches a constant value outside the boundary layer because omega is zero there. Then use just the .99 criterion.

1.10 Streamline

steamline
the path that a particle of zero mass would take through the fluid domain
(no term)
steady-state flow is assumed, even with a transient simulation
(no term)
The path is calculated using a Runge-Kutta method of vector variable integration with variable timestep control

Therefore, when calculating streamlines for a solution for the first time, start by plotting a small number of streamlines and then increase the number of streamlines until the best generation time vs. detail ratio is found.

1000220streamlineanimation.pdf

1.11 Mark cells with high mass imbalance residuals?

go to Mark/Adapt –> Iso-value… -> Residuals -> Mass imbalance and click compute set min and max value >> mark manage display

Refer to >> How to display cells with high turbulent viscosity ratio.pdf

1.12 How to display Pressure and Velocity gradients?

Static Pressure and the X, Y, Z Velocity gradients can be displayed through Contour Plots (Postprocessing/Graphics/Contours). The gradient values are released from memory after being used in the plot calculation and thus cannot be accessed; however, some variables, such as Static Temperature will still be able to be accessed afterwards by executing the following text command: "solve/set/expert" and answering ‘Yes’ to the prompt “Keep temporary solver memory from being freed?”. In addition, immediately after displaying the contour plot of your desired variable, in the Contours dialog box, under the ‘Contours of’ dropdown menu, select Adaption…, then Adaption Space Gradient, then click Display. This will display the spatial gradient of the variable plotted immediately before.

1.13 Write pressure gradient history in a file

1119.pdf

1.14 Spanwise forces on an aircraft wing using UDF on Windows OS

#1366 FLUENT can report the total force that is acting on the wing. But sometimes, it is required to know the spanwise report of forces acting on that wing. This can be achieved by using the UDF included below.

Details on usage of the UDF are given in the UDF manual guide.

Steps:

  1. open case file, check the ID of solid sourface ( cylinder/blade surface)
  2. check the spanwise direction(such as z direction)
  3. check the force direction( such lift, y axis, drag, x axis)
  4. Make the appropriate changes in the UDF +(Axis direction, Force direction, no. of strips and ID of the wall zone)
  5. lauch Fluent from command line
    • Open the Start screen (press Windows button on your keyboard, Type 'VS2013 x86 x64 Cross Tools Command Prompt'
    • Navigate to your working folder, i.e. the folder where your case and data files are (.cas & .dat).
    • Type “fluent” in the command window to launch Fluent or “C:\Program Files\ANSYS Inc\v161\fluent\ntbin\win64\fluent.exe” (quotation signs included, in case of standard installation; fluent version: 16.1).
    • Make sure that on the Environment tab 'Setup Compilation Environment for UDF' is ticked. The default address is fine.
  6. Compile the UDF using Define->User-Defined->Functions->Compiled…
  7. Read in the case file
  8. Increase no. of User-Defined Memory(UDM) to 1
    • Define->User-Defined->Memory…
  9. Initialize the solution so that UDM is initialized
  10. Read in the data file
  11. Execute the demand function "force" +Define->User-Define->Execute On Demand…
  12. Plot contours of UDM-0 and check if the strips are marked correctly
  13. Fluent will report the spanwise distribution on the Fluent terminal
  14. A file called spanwise-force-report.txt will also be created

1.14.1 The UDF code

main macro: DEFINEONDEMAND 

/*******************************************************/
/* Used to obtain a spanwise force report on wings */
/*******************************************************/

#include "udf.h"

/* Spanwise direction */
#define AXIS_X 0
#define AXIS_Y 1
#define AXIS_Z 0

/* The direction along which the force is to be reported */
#define F_X 1
#define F_Y 0
#define F_Z 0

/* Number of strips into which the wing is to be divided */
/* Increasing this number increases the resolution of reporting */
/* It is however limited by the number of cells present along the spanwise direction */
/* Recommended value: */
/* For unstructured cells: no of cells along spanwise dir. divided by 3 */
/* For structured cells: no. of cells along spanwise dir. divided by 3 */
#define N_STRIPS 20

/* ID of the wing boundary zone */
/* You can get this from the Boundary Conditions Panel */
#define TID 15


/* Make sure you increase no. of User Defined Memory to 1 */
/* The UDM stores the strip number for each cell */
/* You can see contour plots of cell values of UDM-1 to check  if the strips are marked correctly */

DEFINE_ON_DEMAND(force)
{
Domain *d=Get_Domain(1);
int i;
real A[ND_ND];
real axis[3], force_dir[3];
face_t f;
Thread *t=Lookup_Thread(d,TID);
cell_t c;
Thread *tc;
real x[ND_ND];
real NV_VEC(presforce),NV_VEC(viscforce);
real NV_VEC(dforce_pres),NV_VEC(dforce_visc);
real pforce[N_STRIPS],vforce[N_STRIPS];
real press=0.0;
real pressure=0.0;
real visc=0.0;
real min_dist=1000000000;
real max_dist=-1000000000;
real dist=0;
real strip_length;
FILE *fp;

NV_S(A,=,0.0);
NV_S(dforce_pres,=,0.0);
NV_S(dforce_visc,=,0.0);
NV_S(presforce,=,0.0);
NV_S(viscforce,=,0.0);

axis[0]=AXIS_X;
axis[1]=AXIS_Y;
axis[2]=AXIS_Z;

force_dir[0]=F_X;
force_dir[1]=F_Y;
force_dir[2]=F_Z;

fp=fopen("spanwise-force-report.txt", "w+");

begin_f_loop(f, t)
{
F_CENTROID(x, f, t);
dist=NV_DOT(x, axis);
if(dist<min_dist) min_dist=dist;
if(dist>max_dist) max_dist=dist;
}
end_f_loop(f, t)

strip_length=(max_dist-min_dist)/N_STRIPS;

Message("Marking strips...");
begin_f_loop(f,t)
{
F_CENTROID(x,f,t);
c=F_C0(f, t);
tc=THREAD_T0(t);
dist=NV_DOT(x, axis);
C_UDMI(c,tc,0)=(int)((dist-min_dist)/strip_length);
}
end_f_loop(f,t)

Message("Donen");

for(i=0;i<N_STRIPS;i++) {
pforce[i]=0;
vforce[i]=0;
}

begin_f_loop(f,t)
{
F_AREA(A,f,t);
c=F_C0(f, t);
tc=THREAD_T0(t);
pforce[(int)C_UDMI(c, tc, 0)]+=F_P(f, t) * NV_DOT(A, force_dir);
vforce[(int)C_UDMI(c, tc, 0)]+=(-1*(NV_DOT(F_STORAGE_R_N3V(f,t,SV_WALL_SHEAR), force_dir)));
}
end_f_loop(f,t)

Message("nnPressure forcet Viscous force");
for(i=0;i<N_STRIPS;i++)
{
press+=pforce[i];
visc+=vforce[i];
Message("n %g t %g",pforce[i],vforce[i]);
fprintf(fp, "%d t %f t %f n", i, pforce[i], vforce[i]);
}
Message("n---------- -----------");
Message("n %gt%g",press,visc);
fclose(fp);
}

1.14.2 new code

I rewrote the UDF in order to be more robust and to be compatible with a parallel execution.

The UDF mark faces with strip number and store that number in user defined memory UDM-0.

There may be some problems if strip number exceeded the size of the array that was declared.

So I increased the size of the array pforce and vforce.

Basically, it transfer all the information from each partition to node 0 and write to the file or to standard output from node 0.

I try to add some comments but some knowledge of UDF is required to understand each line. Maybe you'll find his new UDF more difficult to read because I added parallel compatibility. It also work in serial.

  • FP(f,t) : The pressure force of an individual facet
  • FSTORAGERN3V(f,t,SVWALLSHEAR) : The wall shear stress vector

and direction is opposite to the direction in the solver. That"s the reason one multiply by -1.0.

We need to update the solution.

I tested my UDF on a case in serial and in parallel. And I checked that resulting text file was the same.

/*******************************************************/
/* Used to obtain a spanwise force report on wings */
/*******************************************************/
#include "udf.h"
#include "para.h"
/* Spanwise direction */
#define AXIS_X 0
#define AXIS_Y 1
#define AXIS_Z 0
/* The direction along which the force is to be reported */
#define F_X 0
#define F_Y 0
#define F_Z 1
/* Number of strips into which the wing is to be divided */
/* Increasing this number increases the resolution of reporting */
/* It is however limited by the number of cells present along the spanwise direction
*/
/* Recommended value: */
/* For unstructured cells: no of cells along spanwise dir. divided by 3 */
/* For structured cells: no. of cells along spanwise dir. divided by 3 */
#define N_STRIPS 20
/* ID of the wing boundary zone */
/* You can get this from the Boundary Conditions Panel */
#define TID 7

/* Make sure you increase no. of User Defined Memory to 1 */
/* The UDM stores the strip number for each cell */
/* You can see contour plots of cell values of UDM-1 to check */
/* if the strips are marked correctly */
DEFINE_ON_DEMAND(force)
{
#if !RP_HOST
	Domain *d=Get_Domain(1);
	int i;
	real A[ND_ND];
	real axis[3], force_dir[3];
	face_t f;
	Thread *t=Lookup_Thread(d,TID);
	cell_t c;
	Thread *tc;
	real x[ND_ND];
	real NV_VEC(presforce),NV_VEC(viscforce);
	real NV_VEC(dforce_pres),NV_VEC(dforce_visc);
	real pforce[N_STRIPS+1],vforce[N_STRIPS+1];
	real press=0.0;
	real pressure=0.0;
	real visc=0.0;
	real min_dist=1000000000;
	real max_dist=-1000000000;
	real dist=0;
	real strip_length;
	FILE *fp;
	NV_S(A,=,0.0);
	NV_S(dforce_pres,=,0.0);
	NV_S(dforce_visc,=,0.0);
	NV_S(presforce,=,0.0);
	NV_S(viscforce,=,0.0);
	axis[0]=AXIS_X;
	axis[1]=AXIS_Y;
	axis[2]=AXIS_Z;
	force_dir[0]=F_X;
	force_dir[1]=F_Y;
	force_dir[2]=F_Z;
	
#if PARALLEL
if (I_AM_NODE_ZERO_P)
	fp=fopen("spanwise-force-report.txt", "w+");
#else
    fp=fopen("spanwise-force-report.txt", "w+");
#endif
	begin_f_loop(f, t)
	if (PRINCIPAL_FACE_P(f,t))
	{
		F_CENTROID(x, f, t);
		dist=NV_DOT(x, axis);
		if(dist<min_dist) min_dist=dist;
		if(dist>max_dist) max_dist=dist;
	}
	end_f_loop(f, t)
	
#if PARALLEL
/* get global min and max on node 0*/ 
	min_dist=PRF_GRLOW1(min_dist);
	max_dist=PRF_GRHIGH1(max_dist);
#endif	
	strip_length=(max_dist-min_dist)/N_STRIPS;
	Message0("Marking strips...");
	
	begin_f_loop(f,t)
	{
		F_CENTROID(x,f,t);
		c=F_C0(f, t);
		tc=THREAD_T0(t);
		dist=NV_DOT(x, axis);
		/* assign strip number to UDM-0 */
		C_UDMI(c,tc,0)=(int)((dist-min_dist)/strip_length);
	}
	end_f_loop(f,t)
	Message0("Done\n");
	
	for(i=0;i<=N_STRIPS;i++)
	{
		pforce[i]=0;
		vforce[i]=0;
	}
 
	begin_f_loop(f,t)
	if (PRINCIPAL_FACE_P(f,t))
	{
		F_AREA(A,f,t);
		c=F_C0(f, t);
		tc=THREAD_T0(t);
		i=(int)(MIN(N_STRIPS,C_UDMI(c, tc, 0)));
		/* report force if thread id TID corresponds to a wall */
     if (THREAD_TYPE(t)==THREAD_F_WALL)
	 {
/* pressure force of the strip i along vector (F_X,F_Y,F_Z) */		
		pforce[i]+=F_P(f, t) * NV_DOT(A, force_dir);
/* viscous force of the strip i along vector (F_X,F_Y,F_Z) */				
		vforce[i]+=(-1*(NV_DOT(F_STORAGE_R_N3V(f,t,SV_WALL_SHEAR),force_dir)));	
		/*Message("CHECK: %d\n",sizeof(vforce) / sizeof(vforce[0]));*/	
	 }
	 else
	 {
	 if (f==0) Message0("No forces will be reported since, zone of ID %d is not a wall",TID);
	 }
	}
	end_f_loop(f,t)
	Message0("\n\nPressure force\t Viscous force");
	
#if PARALLEL
/* Get forces of each partition on node 0 */
   	for(i=0;i<=N_STRIPS;i++)
	{
		pforce[i]=PRF_GRSUM1(pforce[i]);
		vforce[i]=PRF_GRSUM1(vforce[i]);
	}
/* report force on node 0 in parallel- Node0 should have a disk system otherwise UDF has to be rewritten to exchange data through the host */	
if (I_AM_NODE_ZERO_P)
{
	for(i=0;i<=N_STRIPS;i++)
	{
		press+=pforce[i];
		visc+=vforce[i];
		Message0("\n %g \t %g",pforce[i],vforce[i]);
		fprintf(fp, "%d \t %f \t %f \n", i, pforce[i], vforce[i]);
	}
	Message0("\n---------- -----------");
	Message0("\n %g\t%g",press,visc);

	fclose(fp);
}
#else
    for(i=0;i<=N_STRIPS;i++)
	{
		press+=pforce[i];
		visc+=vforce[i];
		Message("\n %g \t %g",pforce[i],vforce[i]);
		fprintf(fp, "%d \t %f \t %f \n", i, pforce[i], vforce[i]);
	}
	Message("\n---------- -----------");
	Message("\n %g\t%g",press,visc);

	fclose(fp);
#endif
	Message0("\nUDF DONE\n");
#endif	/* !RP_HOST */
}

1.15 Spanwise forces on an aircraft wing using UDF on Linux OS

refer to running parallel UDF on linux os http://ansysofkemin.blogspot.co.uk/2017/10/running-parallel-udf-on-linux-os-hpc.html

  1. setup the directory structure
  2. build the UDF library
  3. load and Hook the UDF in the Journal file
  4. run the jurnal file

1.16 How to create pictures and videos on-the-fly in Fluent?

Answer: With a set of execute commands it is possible to write out pictures every x iteration/time step and create a video after simulation. This will save a lot of disc space because dat- and cdat-files are no longer needed for the purpose to create videos after simulation. Alternatively, one can create a "solution animation" on-the-fly. This allows the user to change the view/camera position for the final video - w/o saving a lot of dat- or cdat-files.

1.17 Visualize bubble separation

Ensight 3d wing https://www.youtube.com/watch?v=L_mOZz4uMvU

1.18 Why are my results inaccurate?

Firstly, is the simulation correctly set up? Does your model include all relevant physics? Has your solution converged to a reasonable value (see above)? For a simple analysis, an RMS residual of 1E-5 should be sufficient. Keep lowering your residual value until the solution no longer converges monotonically. Perform a sensitivity analysis on the relevant features of your simulation. In its simplest form, a sensitivity analysis involves altering one parameter (e.g., grid refinement) while keeping the other parameters constant and observing the effect this parameter has upon the simulation accuracy (e.g., by monitoring skin friction). The purpose of a sensitivity analysis is to obtain a simulation that is both accurate and cost effective. For just about every simulation a sensitivity analysis on mesh size and convergence tolerance is required. Depending on what you are modelling other sensitivity analyses may be required, covering things such as: Domain size (how close can the boundaries be to the region of interest?) Grid density Grid quality Grid type (e.g., structured versus unstructured), Spatial discretization (or advection) scheme Temporal discretisation scheme Turbulence model, Turbulence numerics, Turbulence intensity, Boundary condition type, Any empirical models used, such as heat/mass transfer coefficients, drag laws, free surface model, Two- versus three-dimensional simulations, and timestep size. Perform an error analysis: Calculate your discretization errors. A calculation should be performed on each quantitative parameter of interest (e.g., lift). The Journal of Fluids Engineering provides an explanation of predicting discretization errors here: Calculating Discretization Errors. In terms of errors, discretization errors are the greatest concern to a CFD modeler. Discretization errors are those that occur as a result of modeling the governing flow equations as algebraic expressions in a discrete space-time domain. Discretization errors should approach zero as the grid is resolved (i.e., as more grid points are generated and the solution becomes grid-independent). However, discretization errors also depend on grid quality (e.g., aspect ratio, orthogonality, skew). A well-refined grid is not necessarily a high quality grid and, thus, grid resolution and quality sensitivity analyses are recommended. Physical approximation errors can cause significant accuracy issues. This type of error is caused by physically approximating something to simplify the simulation (e.g., modeling a real fluid as an ideal gas). Computer round-off errors are caused when a computer stores floating point data. Depending on the accuracy your computer stores the data at, some round-off may occur. This type of error is usually negligible compared to the others and can be easily decreased by running the simulation in Double Precision mode. Iterative convergence errors are due to the resolution of the RANS equations. Since the Navier-Stokes equations are not solved directly, the solution iterates until a specified convergence criterion is reached (e.g., a residual value of a certain size is reached). This type of error can be decreased by allowing convergence to an acceptable level (see above). Usage errors are caused by user specifications. For example, the user may run a turbulent simulation as laminar or an unsteady simulation as steady. These errors can be monitored via a convergence study. And finally - consider the results are you comparing against. Are they accurate? How do you know? Are you modelling EXACTLY the same conditions as the results you are comparing against?

1.19 References

1.20 Are my results publishable?

The Journal of Fluids Engineering provides a good guide of recommendations regarding publishing your data. Journal of Fluids Engineering Numerical Accuracy

2 k in Park Model using Gnuplot

input

# One Rotor, front, eldad blade
# TSR 5, TI =15%
#X/D	X   half width,	Ux	U	Ux/U
0.5	0.23	0.275	0.363	0.6	0.605
1	0.46	0.93	0.353	0.6	0.588333333
2	0.92	0.316	0.383	0.6	0.638333333
3	1.38	0.334	0.41	0.6	0.683333333
4	1.84	0.349	0.441	0.6	0.735
5	2.3	0.365	0.462	0.6	0.77
6	2.76	0.379	0.48	0.6	0.8
7	3.22	0.39	0.497	0.6	0.828333333
8	3.68	0.405	0.509	0.6	0.848333333
9	4.14	0.42	0.52	0.6	0.866666667
10	4.6	0.431	0.529	0.6	0.881666667

code

# fit
% normalized u, f(x) = (U_inf - U_delta)/U_inf
% x, normalized distance, x=X/D
% d, diameter of rotor
C_T=0.7133
U_inf= 0.6
a = U_inf*(1-sqrt(1-C_T))
d=0.46
f(x) = (U_inf - a*d*d/((d+2*k*x)*(d+2*k*x)))/U_inf
fit f(x) 'area_averaged_axial_mean_velocity_TI_1.txt' using 1:6 via k
# plotting
set terminal postscript eps font 24
set out 'k_fit_ti_1_tsr5.eps'
set autoscale
unset log
unset label
unset pm3d
set xtic auto
set ytic auto
unset grid
# set title 'Normalized velocity recover in the wake'
set xlabel  "Normalized axial distance, X/D"
set xrange [*:12]
# r0 initial pulse
set yrange [0.5:1]
set ylabel "Normalized mean axial velocity, ~U{0.8-} / U{/Symbol \245}"
set style line 1 lt 1 lc rgb "black" lw 4 pt 1 ps 2
set style line 2 lt 2 lc rgb "black" lw 4  pt 3 ps 2
set style line 3 lt 3 lc rgb "black" lw 4  pt 5 ps 2
set style line 4 lt 4 lc rgb "black" lw 4  pt 7 ps 2
set style line 5 lt 5 lc rgb "black" lw 4
set style line 6 lt 6 lc rgb "brown" lw 4
k_value = sprintf("k = %.3f", k)
set label 1 at 1, 0.85  k_value
set key at graph 0.9, 0.3
set key spacing 1
plot 'area_averaged_axial_mean_velocity_TI_1.txt' using 1:6 ls 1 with points title 'RANS', f(x) lw 3 title "best fitted"

output

3 Pressure coefficient of a wing/blade

for airfoil Cp = (p-p∞)/q∞

  • q∞ = 0.5 ρ v∞2

for imcompressible flow,inviscid, based on Bernoull's eqs

Cp = 1-(v/v∞)2 (eq. 3.38, anderson, 6th)

  • this is derived from Bernoulli's equation

tidal turbine

  1. σ - local pressure coefficient (eq 6, lust, 2017)

σ = 2(p∞ - pl)/(ρ Urel2)

  • Pl : local pressure at a given point on the surface of the turbine blade
  • Uref : relative seen by blade profile
  • Uref = U∞2 + (ω r)2

Ronsten, G., (1991). ¡®Static pressure measurements in a rotating and a non-rotating 2.35 m wind turbine blade¡¯. Proceedings of the EWEC Conference

Definitionn of pressure coefficient for a blade profile ( foil)

\begin{equation} C_{pr} = 2(p- p_\infty)/(\rho U_{rel}^2) \end{equation}
  • P∞ : free stream pressure
  • Uref : relative seen by blade profile, unknown, you can use the following approximate:
\begin{equation} U_ref = U_\infty^2 + (\omega r)^2 \end{equation}

How to

input: fluent case, data files

software: CFD-POST ansys, Matlab

  1. get pressure distrition on a blade profile ( such as 0.9R station) – CFD-post a) open workbench, add fluent, import case and data files b) add cfd-post
    • use Macro calculator

    Tools > Macro Calculator

cp_polar_plot_cfd_post.png change y coordinate variable to pressure (default is pressure coefficient, Cpr) the default Cpr is not right, as 'velocity' defined in CFD-POST is 'relative velocity' you need to process it in Matlab

output shuld look like:

#X [ m ] Pressure [ Pa ] 0.0138593884 -125.609001 0.0138394814 -126.603676 0.0138003938 -128.664795

  1. pressure coefficient in Matlab

delete header and comments of ascii files,

Before:

#X [ m ] Pressure [ Pa ] 0.0138593884 -125.609001 0.0138394814 -126.603676 0.0138003938 -128.664795

After:

#X [ m ] Pressure [ Pa ] 0.0138593884 -125.609001 0.0138394814 -126.603676 0.0138003938 -128.664795

a) pressurecoeff.m – code for get pressure coefficient

function cp= pressure_coeff(p)
% v : velocity
% p :  pressure
% p0 :   reference pressure
% rho : density
% omega: angular velocity, rad/s
% r : radius
  r=0.207;
v = 0.6;
p0=0;
omega = 13.043;
rho = 998.2;
U = v^2 + (omega*r)^2;
cp=2*(p-p0)/(rho*U)

b) To get a file with pressure coefficient vs. x/c, use Cpr.m Matlab script to save you time

% plot and save 'pressure coefficient vs. x/c'
load pressure_9r_7.1m.out;
pressure=pressure_9r_7_1m(:,2);
x=pressure_9r_7_1m(:,1);
x=x+0.0139;
x_nor=x/0.0276;
cp=pressure_coeff(pressure);  % user defined function
plot(x_nor,cp);
cp_x = [x_nor cp];
save -ascii presure_coeff_7.1m.txt cp_x

Plotting with Gnuplot

#set terminal jpeg
#set output 'coeff_epp_gci.jpg'
#set terminal png
#set output 'pressure_tsr5.png'
set terminal postscript eps font 24
set out 'pressure_coeff_tsr5.eps'
#set terminal x11
set autoscale
unset log
unset label
unset pm3d
set key at graph .7, .4
set key spacing 1
set xtic auto
set ytic auto
set xlabel "x/c - Normalized chord length"
set xrange [*:*]
# r0 initial pulse
set yrange [*:*]
 
set ylabel "C_{pr} - Pressure coefficient"
set style line 1 lt 1 lc rgb "black" lw 4 pointtype 6
set style line 2 lt 2 lc rgb "black" lw 4 pointtype 2
set style line 3 lt 3 lc rgb "black" lw 4
set style line 4 lt 4 lc rgb "red" lw 4
set style line 5 lt 5 lc rgb "black" lw 4
set style line 6 lt 6 lc rgb "brown" lw 4
set pointsize 2
set bars 3
plot "presure_coeff_3.3m.txt" using 1:2  t  "3.5M" ls 1 with lines,\
"presure_coeff_5.4m.txt" using 1:2  t  "5.4M" ls 2 with lines,\
"presure_coeff_7.1m.txt" using 1:2  t  "7.1M" ls 3 with lines

4 Grid/mesh independence   GCI

keywords: Richardson's extrapolation, Grid convergence index a summary of Richardson's extrapolation is here

requirement: GCI < 5%

a summary of GCI from nasa web , local downloaded file is here ( print version is in BEM file folder)

4.1 Richardson extrapolation

4.2 grid refinement ratio

  • Hexa mesh –>> grid refinement ratio
    • double nodes along each coordinates (x, y,z)
  • Tetra mesh –>> effective grid refinement ratio

Definitions:

\[ r_{ij} = h_i/h_j \]

  • r: grid refinement ratio,
  • h: grid spacing
  1. effective grid refinement ratio

    For tetra mesh, the effective grid refinement ratio is defined as:

    \begin{equation} r_e =( \frac{N_1}{N_2}) ^{1/D} \end{equation}

    Where \( N \) is the total number of grid point and \( D \) is the dimension of the flow domain.

4.3 Example of Grid convergence study

The example is from here, The Fortran 90 program verify.f90 was written to carry out the calculations associated with a grid convergence study involving 3 or more grids The program is compiled on a unix system through the commands:

f90 verify.f90 -o verify

It reads in an ASCII file (prD.do) through the standard input unit (5) that contains a list of pairs of grid size and value of the observed quantity f.

input data format: 1.0 0.97050 2.0 0.96854 4.0 0.96178 verify < prD.do > prD.out

It assumes the values from the finest grid are listed first. The output is then written to the standard output unit (6) prD.out. The output from the of {\tt verify} for the results of Appendix A are:

#+BEGINEXAMPLE

— VERIFY: Performs verification calculations —

Number of data sets read = 3

Grid Size Quantity

1.000000 0.970500 2.000000 0.968540 4.000000 0.961780

Order of convergence using first three finest grid and assuming constant grid refinement (Eqn. 5.10.6.1) Order of Convergence, p = 1.78618479

Richardson Extrapolation: Use above order of convergence and first and second finest grids (Eqn. 5.4.1) Estimate to zero grid value, fexact = 0.971300304

Grid Convergence Index on fine grids. Uses p from above. Factor of Safety = 1.25

Grid Refinement Step Ratio, r GCI(%) 1 2 2.000000 0.103080 2 3 2.000000 0.356244

Checking for asymptotic range using Eqn. 5.10.5.2. A ratio of 1.0 indicates asymptotic range.

Grid Range Ratio 12 23 0.997980

— End of VERIFY —

#+END _EXAMPLE

4.4 calculation steps

  1. Complete at least 3 simulations (Coarse, medium, fine) with a constant refinement ratio, r, between them (in our example we use r=2)
  2. Choose a parameter indicative of grid convergence. In most cases, this should be the parameter you are studying. ie if you are studying drag, you would use drag.
  3. Calculate the order of convergence, p, using:
\begin{equation} p=ln(\frac{f_3 - f_2}{f_2- f_1}) / \ln (r) \end{equation}

where \( f_i \) is the solution at different meshes, f1 is fine grid, \( r \) is grid refinement ratio.

  1. Perform a Richardson extrapolation to predict the value at h=0
\begin{equation} f_e = f_1 + \frac{f_1 -f_2 }{r^p - 1} \end{equation}

fe, exact numerical value ( continuum value at zero grid spacing)

  1. Calculate grid convergence index (GCI) for the medium and fine refinement levels
\begin{equation} GCI_{fine} = \frac{F_s \vert \epsilon \vert }{r^p - 1} \end{equation}

where \( F_s \) is a safety factor. the recommended value is 3 for two grids comparisons and 1.25 for three or more grids comparisons.

  1. Ensure that grids are in the asymptotic range of convergence by checking: \frac{GCI2,3}{rp × GCI1,2} \approxeq 1

4.6 References

Roache, P. J. Perspective: A Method for Uniform Reporting of Grid Refinement Studies, Journal of Fluids Engineering, Vol. 116, 1994; 405-413.

Roache, P. J. Quantification of Uncertainty in Computational Fluid Dynamics, in Annual Review of Fluid Mechanics

Roache, Patrick J. Verification and validation in computational science and engineering. Vol. 895. Albuquerque, NM: Hermosa, 1998.

5 Cp calculation - post processing

input: unscaled moment of one blade output: power coefficient of a 3-blades wind/tidal turbine

function cp=coeff(m) % v : velocity % m : unscaled moment of one blade % rho : density % cp : power coefficient of 3 blades % omega: angular velocity, rad/s v = 0.6; A = 0.166; omega = 11.087; rho = 998.2; cp=(3*m*omega)/(0.5*rho*v*v*v*A)

6 Energy Spectrum for LES simulations in CFX:

Step 1: Achieve a statistically stable or averageable state, where turbulent instabilities due to the initial condition effects are convected out of the system. This might be a minimum of 1 through flow cycle.

Step 2: Create the following expressions for TKE.

C = 0.18 #(for Smagorinski, C = 0.18 / WALE, C = 0.5)

delta = (Volume of Finite Volumes)(1/3)

resTKE = 0.5*((u-u.Trnavg)2 + (v-v.Trnavg)2 + (w-w.Trnavg)2)

unresTKE = ((C*delta*Shear Strain Rate)2)/0.3

totalTKE = Total TKE = resTKE + unresTKE

Create an Additional Variable > TKE > Variable Type: Unspecified; Units: [m2 s-2]; Tensor Type: Scalar.

Inside the domain panel: Domain > Fluid Models > Additional Variable > TKE > Option: Algebraic Equation; Add. Var. Value > totalTKE

Step 3: (In CFX-Pre) Create monitor points at a few locations in the flow field (at specific regions of interest). Monitor TKE at those points.

Step 4: (In CFX-Solver) Monitor TKE at the all points (separately). (i) Workspace > New Monitor > OK; (ii) Plot Lines > USER POINT > TKE > Monitor Point Location > OK; (iii) Range Settings > Plot Data By > Simulation Time.

Step 5: (In CFX-Solver) Run for sufficiently long time so that you get a sharper energy spectrum plot. Once the simulation is complete, right click on the screen to save the monitor plot as a csv file.

Step 6: (In CFD-Post) (i) Insert > Chart > OK; (ii) Data Series > File (Select the saved csv file); (iii) General > Fast Fourier Transform; (iv) X Axis > Frequency > Log Scale; (v) Y Axis > Magnitude > Log Scale; (vi) Line Display > FFT Line Type > Lines; (vii) Apply. This will give you the energy spectrum.

Hope this helps. Let us know if you have any issues.

Rahul https://studentcommunity.ansys.com/thread/how-to-make-the-spectrum-for-turbulent-kinetic-energy/

7 coefficient of pressure–Tecplot

This video compares a CFD solution with experimental data. We will be creating a classic Coefficient of Pressure (Cp) plot. This process will include slice extraction, new frame creation, and XY line manipulation. Experimental data was gathered using pressure taps along the wing in a wind tunnel.

This third video in our External Flow Getting Started series, Comparing a CFD Simulation with Experimental Data, uses the ONERA M6 wing data set. You can find the link to the data in the description below. The ONERA M6 wing was tested by NASA in a wind tunnel at four different Mach numbers and various angles of attack. This is now a classic CFD validation case for external flows because it has a simple geometry, complex flows, and experimental data to validate against.

Thanks for watching! https://www.tecplot.com/2015/10/30/external-flow-comparing-cfd-simulation-with-experimental-data/

8 2D Streamline

 

8.1 side view streamline

streamline_sideview_tsr5_4D_eldad_mrf.png a recirculation bubble downstream of the rotor Effects of Circulation Control on Power Production for Large Scale Wind Turbines

  • qantities to describe vortex?

9 vortex vs cavitation

10 Velocity Vector at a blade station

Problem 0.5R_velocity_vector.png the vector near blade is overlapping with others answer:

Creating a plane with a hole around the blade is not really common. Why don't you simply reduce the symbol size of the vector? In case you really want to cut out a hole, you could use the iso-clip option using the wall distance as visibility parameter.

11 Seperation Point

  • friction coefficient is zero
  • slope of pressure coefficient is zero

12 Extract local angle of attack on wind turbine blades

keywords: local AOA, wind turbine

  1. The inverse BEM method [13, 14], which uses the pre-determined local forces to calculate

the local induction.

  1. Using CFD data to obtain the annular average of the axial velocity (and thereby the

induction a) at a given radial position in the rotor plane [15, 16, 17].

  1. Using CFD data to obtain the the axial velocity at a given radial position in the rotor plane

at the location of the blade (aB). This method is similar to method (2.).

  1. Determination of AOA by comparison of high-pressure-side CP distributions of a 3D case

with a 2D case with a known AOA 1. Methods (1.) and (4.) require sectional force coefficients and pressure distributions, respectively, making them suitable with experimental data where sectional force coefficients and pressure distributions are often the only data available. On the other hand, methods (2.) and (3.) require detailed flow field information upstream and downstream of the rotor, making them more suitable where full rotor CFD data is available. 2

12.1 References

2018 Eva Jost Angle of Attack (AoA), turbine, Extracting the angle of attack on rotor blades from CFD simulations 2017 Jost, Eva, et al. "Extracting the angle of attack on rotor blades from CFD simulations." Wind Energy (2017). 2018 Rahimi, Hamid, et al. "Evaluation of different methods for determining the angle of attack on wind turbine blades with CFD results under axial inflow conditions." Renewable Energy 125 (2018): 866-876.

Guntur2014An evaluation of several methods of determining the local angle of attack on wind turbine blades 2018 Eva Jost Angle of Attack (AoA), turbine, Extracting the angle of attack on rotor blades from CFD simulations 2017 Jost, Eva, et al. "Extracting the angle of attack on rotor blades from CFD simulations." Wind Energy (2017). 2018 Rahimi, Hamid, et al. "Evaluation of different methods for determining the angle of attack on wind turbine blades with CFD results under axial inflow conditions." Renewable Energy 125 (2018): 866-876.

13 Q-Criterion

example:

Qcriterion of a wind turbine QCriterion.png

isosurface
a 3D surface representation of points with equal values in a 3D data distribution. Is the 3D equivalent of a contour line.

13.1 Q-criterion– Hunt, Wray & Moin 1988

J. Fluid. Mech. 366, 87–108. Hunt, J. C. R., Wray, A. & Moin, P. 1988 Eddies, stream, and convergence zones in turbulent flows. Center for Turbulence Research Report CTR-S88.

velocity gradient, ∇ v ∇ v = S + Ω where, S is rate-of-strain tensor, and Ω is vorticity tensor

\begin{equation} S= 0.5[\nabla v + (\nabla v)^T] \Omega = 0.5[\nabla v - (\nabla v)^T] \end{equation}

Q-criterion: Q=0.5(||Ω||2 - ||S||2)

Q> 0

13.2 Q criterion-paraview

  1. import results
  2. gradient on an unstructured dataset

Q_criterion_setup_paraview.png

13.3 Q-Criterion with Tecplot

https://www.tecplot.com/blog/2017/05/15/calculate-q-criterion-tecplot-360/

>>> analysze/calculate variables/Q criterion http://download.tecplot.com/videos/3602017R2-QCriterion.mp4 Q-Criterion is an important calculation used to identify vortices. In this video we’ll show you how to calculate Q-Criterion, plot the results, and compare the performance of Plot3D and SZL file formats for this work.

In this example, we are using the final time-step of a transient simulation of a wind turbine. An overset mesh was used and is composed of 5863 zones with a total of 263 million elements. The output is in Plot3D file format and the grid and solution files combine for a total of 21 GB.

To calculate Q-Criterion:

>> Analyze > Calculate Variables and select Q Criterion from the list. The Calculate Variables dialog has a unique feature called Calculate-on-Demand – if you check this toggle, the formula will simply be registered and the calculation will not occur until it is required. This can save a significant amount of time particularly for a data set such as this which has many zones, because the calculation will only be performed on the zones required for the desired plot.

Use an iso-surface, at a positive value near zero, to view the resulting Q-Criterion calculation. Some adjustment of the iso-surface value will be required to see the vortices. Increasing Q reduces the complexity of the iso-surface, but too high of a value makes an iso-surface that is too sparse. It is important to find a value that results in an iso-surface that is neither too dense nor too sparse.

When run in batch mode, the total time to load the data, calculate Q and generate an image with an iso-surface at Q=0.01, without Calculate-on-Demand, was 517 seconds for the Plot3D file and 409 seconds for the SZL file. Using Calculate-on-Demand was 27% faster, taking 372 seconds for the Plot3D file and 299 seconds for the SZL file. Additionally the SZL file is 35% smaller at only 13GB https://www.youtube.com/watch?v=LmqQtTAsTiA

14 tangential wall shear in CFD-Post

Wall Shear is provided as a vector variable in CFD-Post. Taking the dot product with the wall normal gives the magnitude of the normal shear vector from which the vector components can be constructed. Subtracting the normal shear vector from the Wall Shear gives the tangential shear. Then a User Vector Variable can be defined using the X,Y,Z components of the tangential shear. An example state file that does this is attached.

14.1 tangentialshear.cst

LIBRARY: CEL: EXPRESSIONS: NormalShear = Wall Shear X * Normal X + Wall Shear Y * Normal Y + Wall \ Shear Z * Normal Z NormalShearX = NormalShear * Normal X NormalShearY = NormalShear * Normal Y NormalShearZ = NormalShear * Normal Z TangentialShear = (TangentialShearX2 + TangentialShearY2 + \ TangentialShearZ2)0.5 TangentialShearX = Wall Shear X - NormalShearX TangentialShearY = Wall Shear Y - NormalShearY TangentialShearZ = Wall Shear Z - NormalShearZ END END END

USER SCALAR VARIABLE:Normal Shear X Boundary Values = Conservative Calculate Global Range = Off Component Index = 1 Recipe = Vector Component Variable = Normal Shear END USER SCALAR VARIABLE:Normal Shear Y Boundary Values = Conservative Calculate Global Range = Off Component Index = 2 Recipe = Vector Component Variable = Normal Shear END USER SCALAR VARIABLE:Normal Shear Z Boundary Values = Conservative Calculate Global Range = Off Component Index = 3 Recipe = Vector Component Variable = Normal Shear END USER SCALAR VARIABLE:Tangential Shear X Boundary Values = Conservative Calculate Global Range = Off Component Index = 1 Recipe = Vector Component Variable = Tangential Shear END USER SCALAR VARIABLE:Tangential Shear Y Boundary Values = Conservative Calculate Global Range = Off Component Index = 2 Recipe = Vector Component Variable = Tangential Shear END USER SCALAR VARIABLE:Tangential Shear Z Boundary Values = Conservative Calculate Global Range = Off Component Index = 3 Recipe = Vector Component Variable = Tangential Shear END USER VECTOR VARIABLE:Normal Shear Boundary Values = Conservative Calculate Global Range = Off Recipe = Expression Variable to Copy = Pressure Variable to Gradient = Pressure X Expression = NormalShearX Y Expression = NormalShearX Z Expression = NormalShearZ END USER VECTOR VARIABLE:Tangential Shear Boundary Values = Conservative Calculate Global Range = Off Recipe = Expression Variable to Copy = Pressure Variable to Gradient = Pressure X Expression = TangentialShearX Y Expression = TangentialShearY Z Expression = TangentialShearZ END

15 wake

parameters for wake:

  • mean axial velocity vs axial distance (X)
  • velocity Contours
  • the shear region 
The wake was bordered by
a sharp transition zone distinguishable by a thin annular area with steep velocity gradient, the shear region.

(Heiner Schümann et al. / Energy Procedia 35 ( 2013 ) 285 – 296)

15.1 References

  • wake.doc
  • wake.org

16 boundary layer thickness on or around a 3D shape

Boundary layer thickness can be estimated for both CFX and FLUENT results using the Boundary layer "Health" indicator based on total pressure loss calculation.

This parameter is not directly available in CFD-Post but can be easily created via CEL expression as explained in the enclosed pdf.

17 Video/Animation

 

17.1 pressure contour for unsteady flow

keyword: post-processing, for loop, contour

input: N results file

package: Tecplot/Paraview

Method: read data 1 and slice, then save pressure, velocity contours, close data, then do loop for the aforementioned steps.

17.1.1 code

****************************************************************

#### import the simple module from the paraview
from paraview.simple import *
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()

for Num in range(36,47):
    # create a new 'EnSight Reader' and assign it to a variable, 'transientcase'
    transientcase = EnSightReader(CaseFileName='/home/kaiming/Documents/ZJU_Projects/Jet/data/transient_%i.case' %(Num+1))
    transientcase.PointArrays = [ 'density','v', 'pressure', 'temperature']

    # get active view
    renderView1 = GetActiveViewOrCreate('RenderView')
    # set a specific view size
    renderView1.ViewSize = [1022, 837]



    # get color transfer function/color map for 'density'
    densityLUT = GetColorTransferFunction('density')
    densityLUT.LockDataRange = 1


    # get opacity transfer function/opacity map for 'density'
    densityPWF = GetOpacityTransferFunction('density')


    # show data in view
    transientcaseDisplay = Show(transientcase, renderView1)
    # trace defaults for the display properties.
    transientcaseDisplay.ColorArrayName = ['POINTS', 'density']
    transientcaseDisplay.LookupTable = densityLUT
    transientcaseDisplay.GlyphType = 'Arrow'
    transientcaseDisplay.ScalarOpacityUnitDistance = 0.0016420380639339577

    # reset view to fit data
    renderView1.ResetCamera()

    # show color bar/color legend
    transientcaseDisplay.SetScalarBarVisibility(renderView1, False)

    # get opacity transfer function/opacity map for 'density'
    densityPWF = GetOpacityTransferFunction('density')

    # reset view to fit data
    renderView1.ResetCamera()



    #################
    ## slice
    ################
    # create a new 'Slice'
    slice1 = Slice(Input=transientcase)
    slice1.SliceType = 'Plane'
    slice1.SliceOffsetValues = [0.0]

    # init the 'Plane' selected for 'SliceType'
    slice1.SliceType.Origin = [-0.21849990739250558, 0.0, 0.0]

    # Properties modified on slice1.SliceType
    slice1.SliceType.Origin = [0.0, 0.0, 0.0]
    slice1.SliceType.Normal = [0.0, 0.0, 1.0]

    # show data in view
    slice1Display = Show(slice1, renderView1)
    # trace defaults for the display properties.
    slice1Display.ColorArrayName = ['POINTS', 'density']
    slice1Display.LookupTable = densityLUT
    slice1Display.GlyphType = 'Arrow'

    # hide data in view
    Hide(transientcase, renderView1)

    # show color bar/color legend
    slice1Display.SetScalarBarVisibility(renderView1, True)

    # set active source
    SetActiveSource(transientcase)

    # reset view to fit data
    renderView1.ResetCamera()


    # current camera placement for renderView1
    renderView1.CameraPosition = [-0.3656950276430585, -0.000908563692513454, 0.21027127790890924]
    renderView1.CameraFocalPoint = [-0.3656950276430585, -0.000908563692513454, 0.0]
    renderView1.CameraParallelScale = 0.25006857916835856


    # *****************
    # change legend layout, and its font color, position
    # *****************

    # get color legend for 'densityLUT' in view 'renderView1'
    densityLUTColorBar = GetScalarBar(densityLUT, renderView1)

    # Properties modified on densityLUTColorBar
    densityLUTColorBar.AutoOrient = 0
    densityLUTColorBar.RangeLabelFormat = '%.2f'
    ## legend orientation
    densityLUTColorBar.Orientation = 'Horizontal'

    ## legend normalized position
    densityLUTColorBar.Position = [0.3, 0.2]

    # change label color to 'black'
    densityLUTColorBar.LabelColor = [1.0, 1.0, 1.0]

    #  change titile color to 'black'
    densityLUTColorBar.TitleColor = [1.0, 1.0, 1.0]
    # ***************

    # set Background color as 'White'
    renderView1.Background =[1,1,1]


    # get layout
    viewLayout1 = GetLayout()
	# set a specific view size
    renderView1.ViewSize = [1022, 837]

    # current camera placement for renderView1
    renderView1.CameraPosition = [-0.3656950276430585, -0.000908563692513454, 0.21027127790890924]
    renderView1.CameraFocalPoint = [-0.3656950276430585, -0.000908563692513454, 0.0]
    renderView1.CameraParallelScale = 0.25006857916835856


    # save screenshot
    SaveScreenshot('/home/kaiming/Documents/ZJU_Projects/Jet/paraview/tem/d_%s.png' %(Num+1), layout=viewLayout1, magnification=1, quality=100)



    #################
    # pressure contour
    ##################

    # show color bar/color legend
    slice1Display.SetScalarBarVisibility(renderView1, False)

    # set active source
    SetActiveSource(slice1)

    # set scalar coloring
    ColorBy(slice1Display, ('POINTS', 'pressure'))

    # rescale color and/or opacity maps used to include current data range
    slice1Display.RescaleTransferFunctionToDataRange(True)

    # show color bar/color legend
    slice1Display.SetScalarBarVisibility(renderView1, True)




    # get color transfer function/color map for 'pressure'

    pressureLUT = GetColorTransferFunction('pressure')
    pressureLUT.LockDataRange = 1


    # get opacity transfer function/opacity map for 'pressure'
    pressurePWF = GetOpacityTransferFunction('pressure')




    # set active source
    SetActiveSource(transientcase)
    # *****************
    # change legend layout, and its font color, position
    # *****************

    # get color legend for 'pressureLUT' in view 'renderView1'
    pressureLUTColorBar = GetScalarBar(pressureLUT, renderView1)

    # Properties modified on vLUTColorBar
    pressureLUTColorBar.AutoOrient = 0
    pressureLUTColorBar.RangeLabelFormat = '%.2f'
    ## legend orientation
    pressureLUTColorBar.Orientation = 'Horizontal'

    ## legend normalized position
    pressureLUTColorBar.Position = [0.3, 0.2]

    # change label color to 'black'
    pressureLUTColorBar.LabelColor = [1.0, 1.0, 1.0]

    #  change titile color to 'black'
    pressureLUTColorBar.TitleColor = [1.0, 1.0, 1.0]
    # ***************


    # current camera placement for renderView1
    # current camera placement for renderView1
    renderView1.CameraPosition = [-0.3656950276430585, -0.000908563692513454, 0.21027125790890924]
    renderView1.CameraFocalPoint = [-0.3656950276430585, -0.000908563692513454, 0.0]
    renderView1.CameraParallelScale = 0.25006857916835856



    # save screenshot
    SaveScreenshot('/home/kaiming/Documents/ZJU_Projects/Jet/paraview/tem/p_%s.png' %(Num+1), layout=viewLayout1, magnification=1, quality=100)


    ######################
    # #temperature contour
    ######################

    # show color bar/color legend
    slice1Display.SetScalarBarVisibility(renderView1, False)

    # set active source
    SetActiveSource(slice1)

    # set scalar coloring
    ColorBy(slice1Display, ('POINTS', 'temperature'))

    # rescale color and/or opacity maps used to include current data range
    slice1Display.RescaleTransferFunctionToDataRange(True)

    # show color bar/color legend
    slice1Display.SetScalarBarVisibility(renderView1, True)




    # get color transfer function/color map for 'termperature'

    temperatureLUT = GetColorTransferFunction('temperature')
    temperatureLUT.LockDataRange = 1


    # get opacity transfer function/opacity map for 'temperature'
    temperaturePWF = GetOpacityTransferFunction('temperature')


    # ******
    # legend layout, and its font color, position
    # *****

    temperatureLUTColorBar = GetScalarBar(temperatureLUT, renderView1)

    # Properties modified on vLUTColorBar

    temperatureLUTColorBar.AutoOrient = 0
    temperatureLUTColorBar.RangeLabelFormat = '%.2f'
    ## legend orientation
    temperatureLUTColorBar.Orientation = 'Horizontal'

    ## legend normalized position
    temperatureLUTColorBar.Position = [0.3, 0.2]

    # label color to 'black'
    temperatureLUTColorBar.LabelColor = [1.0, 1.0, 1.0]

    #  change 'titile' color to 'black'
    temperatureLUTColorBar.TitleColor = [1.0, 1.0, 1.0]
    # *****************

    # set active source
    SetActiveSource(transientcase)

    # current camera placement for renderView1
    # current camera placement for renderView1
    renderView1.CameraPosition = [-0.3656950276430585, -0.000908563692513454, 0.21027127790890924]
    renderView1.CameraFocalPoint = [-0.3656950276430585, -0.000908563692513454, 0.0]
    renderView1.CameraParallelScale = 0.25006857917835856

    # save screenshot
    SaveScreenshot('/home/kaiming/Documents/ZJU_Projects/Jet/paraview/tem/t_%s.png' %(Num+1), layout=viewLayout1, magnification=1, quality=100)


    # *****************
    #  Velocity
    # ****************

    # show color legend
    slice1Display.SetScalarBarVisibility(renderView1, False)

    # set active source
    SetActiveSource(slice1)

    # set scalar coloring
    ColorBy(slice1Display, ('POINTS', 'v'))

    # rescale color and/or opacity maps used to include current data range
    slice1Display.RescaleTransferFunctionToDataRange(True)

    # show color bar/color legend
    slice1Display.SetScalarBarVisibility(renderView1, True)



    # get color transfer function/color map for 'v'

    vLUT = GetColorTransferFunction('v')
    vLUT.LockDataRange = 1



    # get opacity transfer function/opacity map for 'v'
    vPWF = GetOpacityTransferFunction('v')


    # *****************
    # change legend layout, and its font color, position
    # *****************

    # get color legend for 'densityLUT' in view 'renderView1'
    vLUTColorBar = GetScalarBar(vLUT, renderView1)

    # Properties modified on vLUTColorBar
    vLUTColorBar.AutoOrient = 0
    vLUTColorBar.RangeLabelFormat = '%.2f'
    ## legend orientation
    vLUTColorBar.Orientation = 'Horizontal'

    ## legend normalized position
    vLUTColorBar.Position = [0.3, 0.2]

    # change label color to 'black'
    vLUTColorBar.LabelColor = [1.0, 1.0, 1.0]

    #  change titile color to 'black'
    vLUTColorBar.TitleColor = [1.0, 1.0, 1.0]
    # ***************

    # set active source
    SetActiveSource(transientcase)

    # current camera placement for renderView1
    # current camera placement for renderView1
    renderView1.CameraPosition = [-0.3656950276430585, -0.000908563692513454, 0.21027125790890924]
    renderView1.CameraFocalPoint = [-0.3656950276430585, -0.000908563692513454, 0.0]
    renderView1.CameraParallelScale = 0.25006857916835856

    # save screenshot
    SaveScreenshot('/home/kaiming/Documents/ZJU_Projects/Jet/paraview/tem/v_%s.png' %(Num+1), layout=viewLayout1, magnification=1, quality=100)

    #### saving camera placements for all active views


	# show color bar/color legend
    slice1Display.SetScalarBarVisibility(renderView1, False)
    # current camera placement for renderView1
    renderView1.CameraPosition = [-0.3656950276430585, -0.000908563692513454, 0.21027127790890924]
    renderView1.CameraFocalPoint = [-0.3656950276430585, -0.000908563692513454, 0.0]
    renderView1.CameraParallelScale = 0.25006857916835856
	Disconnect()
	Connect()
    #### uncomment the following to render all views
    # RenderAllViews()
    # alternatively, if you want to write images, you can use SaveScreenshot(...).

18 Reference

> analysis.doc > Délery, Jean. Three-dimensional separated flows topology: singular points, beam splitters and vortex structures. John Wiley & Sons, 2013.

Footnotes:

1

DEFINITION NOT FOUND.

2

2014 S Guntur, N N Sorensen An evaluation of several methods of determining the local angle of attack on wind turbine blades

Created: 2019-06-19 周三 07:46

Emacs 25.3.1 (Org mode 8.2.10)

Validate

posted @ 2019-06-19 14:49  kaiming_ai  阅读(4062)  评论(0)  编辑  收藏  举报