- CFD, Fluid Flow, FEA, Heat/Mass Transfer
- +91 987-11-19-383
- amod.cfd@gmail.com

Automation

A coollection of scripts, journals and macros in CFD simulations to automate some tasks as well as enhance simulation capabilties. One may need such functions to apply special boundary conditions such as inlet velocity which is function of both time and space. Scripts or macros can be used to repeat the same simulation with changed boundary conditions or can be used to create a checking log where the summary of boundary conditions, solver setting , material propertis and interface/periodicity information can be written out in a log file for self and peer review.

- A STAR-CCM+ macro is a Java program that is compiled and executed within the STAR-CCM+ workspace. A FLUENT UDF is a C program that is compiled and executed within the FLUENT workspace.
- CFX uses a high level programming called CCL or CFX command language. Unlike UDF or JAVA macros, it does not need to be interpreted or compiled. However, for special post processing, commands in PERL and accessing solver program in FORTRAN is required.
- Note that Cp cannot be accessed or changed via UDF in FLUENT.

ICEM CFD | STAR-CCM+ | FLUENT | CFX | OpenFOAM |

Tck/Tk | Java | SCHEME, C | CCL, PERL, FORTRAN | C++ |

- A guide describing syntax used in SCHEME language can be found here (scheme programming).
- A summary of few basic yet key features of the programming languages mentioned above is tabulated below.

Feature | Tcl/Tk | FORTRAN | C | JAVA |

Case-sensitive | Y | N | Y | Y |

Comment | # | C, ! | /* ... */ | // |

End of Statement | ; | Newline character | ; | ; |

Line Continuation | \ | Any character | \ | Not required |

Variable Definition | set x = 1; | real x = 1 | real x = 1; | real x = 1; |

If Loop | if { x == y } { x = x + 1; } | if (x .EQ. y) then x = x + 1 endif | if (x = y) { x = x + 1; } | if (x = y) { x = x + 1; } |

For Loop | for {set j 1} {$j <= $n} {incr j} { ... } | - | for (i=0; i<= 10, i++) { ... } | for (i=0; i<= 10, i++) { ... } |

Arrays | $x(5); | x(5) | x[5]; | x[5]; |

File Embedding | source "common.dat"; | include "common.dat" | #include "common.h"; | import common.class; |

Java doesn't need a "line continuation character" since the end of the line has no significance so long a semi-colon is added at the end of a statement. It acts just as any other whitespace character. In other words, one can end the line immediately after an x = and continue on the assignment (say 10 in x = 10) on next line without any problems. For Tck/Tk used in ICEM CFD, the end of statement is also a newline charcter. However, if more than one statements are to be put on a line, they can be separated by a semi-colon.

Examples: In FLUENT, the in-built macros used to acccess components of the velocity are

- C_U(c,t): x-component of velocity
- C_V(c,t): y-component of velocity
- C_W(c,t): z-component of velocity

- Magnitude: mag($$Velocity)
- X-component: $$Velocity[1]
- Y-Component: $$Velocity[1]...
- The directions X,Y,Z are represented by [0],[1],[2] respectively.

define/models/solver/pressure-based yes ; define/models/viscous/ke-standard yes ke-realizable yes define/models/viscous/near-wall-treatment enhanced-wall-treatment yes ;-----------IDEAL-GAS ---------------- define/materials/change-create air air yes ideal-gas yes polynomial 3 1033.33 -0.196044 3.93365e-4 yes polynomial 3 6.24978e-6 9.73517e-5 -3.31177e-8 yes sutherland two-coeeficient-method two-coefficient-method 1.458e-6 110.4 yes 28.966 no no no no no ; define/boundary-conditions/zone-type inlet pressure-inlet define/boundary-conditions/pressure-inlet inlet no 0 no 0 no 300 no yes no yes 5 0.1 define/boundary-conditions/zone-type outlet pressure-outlet define/boundary-conditions/pressure-outlet outlet no 0 no 0 no 300 no yes no yes 5 0.1 ; define/operating-conditions operating-pressure 101325 define/operating-conditions reference-pressure-location 0 0 0 define/operating-conditions gravity no ; ;Define Discretization Scheme ;0-1st UDS, 1->2nd UDS, 2->Power Law, 4-> QUICK, 6->3rd Order MUSCL solve/set/discretization-scheme density 1 solve/set/discretization-scheme mom 1 solve/set/discretization-scheme k 1 solve/set/discretization-scheme epsilon 1 solve/set/discretization-scheme temperature 1 ; ;Press: 10->Std, 11->Linear, 12-> 2nd Order, 13->Body Force Weighted, 14-> PRESTO! solve/set/discretization-scheme pressure 12 ; ;Flow: 20->SIMPLE, 21->SIMPLEC, 22->PISO solve/set/p-v-coupling 21 ; ;Define Under-Relaxation Factors solve/set/under-relaxation body-force 0.8 solve/set/under-relaxation k 0.8 solve/set/under-relaxation epsilon 0.8 solve/set/under-relaxation density 0.8 solve/set/under-relaxation mom 0.4 (define rf 0.7) solve/set/under-relaxation pressure rf solve/set/under-relaxation turb-viscosity 0.8 solve/set/under-relaxation temperature 1.0 solve/set/limits 50000 150000 250 400 1e-14 1e-20 40000 ; solve/monitors/residual convergence-criteria 1.0e-6 1.0e-6 1.0e-6 1.0e-6 1.0e-6 1.0e-6 solve/monitors/residual plot y solve/monitors/residual print y

Similarly, the CCL file required to do similar operation in CFX is as follows. This file can also be used to set

COMMAND FILE: CFX Pre Version = 12.1 END # Define all the variables in a PERL file and include with 'require' statement. ! require "VarCFX_Pre.pl"; >load mode=new >update >gtmImport filename=$FCFX5, type= Generic, genOpt= -n, units=mm, \ nameStrategy=$Name > update >writeCaseFile filename=$CASE > update FLOW: Flow Analysis 1 &replace SOLUTION UNITS: Angle Units = [rad] Length Units = [m] #Options: [cm], [mm], [in], [ft], [yd], [mile] Mass Units = [kg] #Options: [g], [lb], [slug], [tonne] Solid Angle Units = [sr] Temperature Units = [K] #Options; [R] Time Units = [s] #Options: [min], [hr], [day], [year] END # SOLUTION UNITS: END # FLOW:Flow Analysis 1 > update # ---------------------------------------------------------------------+ LIBRARY: CEL: EXPRESSIONS: CpAir=1005.6 [J kg^-1 K^-1] + 5.6E-3 [J kg^-1 K^-2] * T END END END > update LIBRARY: CEL: EXPRESSIONS: kAir = -3.9333E-4[W m^-1 K^-1] + 1.0184E-4 [W m^-1 K^-2]*T \ -4.8574E-8 [W m^-1 K^-3]*T*T + 1.5207E-11 [W m^-1 K^-4]*T^3 END END END > update LIBRARY: CEL: EXPRESSIONS: MuAir=1.72E-5 [kg m^-1 s^-1] *(T / 293 [K])^0.742 END END END > update LIBRARY: CEL: EXPRESSIONS: AvHeatFlux=areaAve(Wall Heat Flux )@REGION:$WallName END END END > update # ---------------------------------------------------------------------+ LIBRARY: &replace MATERIAL: UserMat Material Description = User-defined Material Material Group = User Object Origin = User Option = Pure Substance Thermodynamic State = Gas PROPERTIES: Option = General Material DYNAMIC VISCOSITY: Dynamic Viscosity = MuAir Option = Value END # DYNAMIC VISCOSITY: EQUATION OF STATE: Molar Mass = $MolM [kg kmol^-1] Option = Ideal Gas END # EQUATION OF STATE: SPECIFIC HEAT CAPACITY: Option = Value Specific Heat Capacity = CpAir Specific Heat Type = Constant Pressure END # SPECIFIC HEAT CAPACITY: THERMAL CONDUCTIVITY: Option = Value Thermal Conductivity = kAir END # THERMAL CONDUCTIVITY: END # PROPERTIES: END # MATERIAL:UserMat END # LIBRARY: > update FLOW: Flow Analysis 1 &replace DOMAIN: Domain_Pipe Coord Frame = Coord 0 Domain Type = Fluid Location = $BlkName DOMAIN MODELS: BUOYANCY MODEL: Option = Non Buoyant END # BUOYANCY MODEL: DOMAIN MOTION: Option = Stationary END # DOMAIN MOTION: MESH DEFORMATION: Option = None END # MESH DEFORMATION: REFERENCE PRESSURE: Reference Pressure = $RefPress [atm] END # REFERENCE PRESSURE: END # DOMAIN MODELS: FLUID DEFINITION: Fluid 1 Material = Air Ideal Gas Option = Material Library MORPHOLOGY: Option = Continuous Fluid END # MORPHOLOGY: END # FLUID DEFINITION:Fluid 1 FLUID MODELS: COMBUSTION MODEL: Option = None END # COMBUSTION MODEL: HEAT TRANSFER MODEL: Option = Thermal Energy END # HEAT TRANSFER MODEL: THERMAL RADIATION MODEL: Option = None END # THERMAL RADIATION MODEL: TURBULENCE MODEL: Option = $TurbModel END # TURBULENCE MODEL: TURBULENT WALL FUNCTIONS: Option = Scalable END # TURBULENT WALL FUNCTIONS: END # FLUID MODELS: END # DOMAIN:Domain_Pipe END # FLOW:Flow Analysis 1 > update FLOW: Flow Analysis 1 DOMAIN: Domain_Pipe &replace BOUNDARY: $WallName Boundary Type = WALL Create Other Side = Off Interface Boundary = Off Location = $WallName BOUNDARY CONDITIONS: HEAT TRANSFER: Fixed Temperature = $WallT [K] Option = Fixed Temperature END # HEAT TRANSFER: MASS AND MOMENTUM: Option = No Slip Wall END # MASS AND MOMENTUM: WALL ROUGHNESS: Option = Smooth Wall END # WALL ROUGHNESS: END # BOUNDARY CONDITIONS: END # BOUNDARY: $WallName END # DOMAIN:Domain_Pipe END # FLOW:Flow Analysis 1 > update FLOW: Flow Analysis 1 DOMAIN: Domain_Pipe &replace BOUNDARY: Inlet Boundary Type = INLET Interface Boundary = Off Location = INLET BOUNDARY CONDITIONS: FLOW DIRECTION: Option = Normal to Boundary Condition END # FLOW DIRECTION: FLOW REGIME: Option = Subsonic END # FLOW REGIME: HEAT TRANSFER: Option = Static Temperature Static Temperature = $InletT [K] END # HEAT TRANSFER: MASS AND MOMENTUM: Mass Flow Rate = 0.1 [kg s^-1] Option = Mass Flow Rate END # MASS AND MOMENTUM: TURBULENCE: Option = Medium Intensity and Eddy Viscosity Ratio END # TURBULENCE: END # BOUNDARY CONDITIONS: END # BOUNDARY:Inlet END # DOMAIN:Domain_Pipe END # FLOW:Flow Analysis 1 > update FLOW: Flow Analysis 1 DOMAIN: Domain_Pipe &replace BOUNDARY: Outlet Boundary Type = OUTLET Interface Boundary = Off Location = OUTLET BOUNDARY CONDITIONS: FLOW REGIME: Option = Subsonic END # FLOW REGIME: MASS AND MOMENTUM: Option = $OutBC Relative Pressure = $OutBCValue [Pa] END # MASS AND MOMENTUM: END # BOUNDARY CONDITIONS: END # BOUNDARY:Outlet END # DOMAIN:Domain_Pipe END # FLOW:Flow Analysis 1 > update FLOW: Flow Analysis 1 &replace INITIALISATION: Option = Automatic INITIAL CONDITIONS: Velocity Type = Cartesian CARTESIAN VELOCITY COMPONENTS: Option = Automatic with Value U = $InitU [m s^-1] V = $InitV [m s^-1] W = $InitW [m s^-1] END # CARTESIAN VELOCITY COMPONENTS: STATIC PRESSURE: Option = Automatic with Value Relative Pressure = $InitP [Pa] END # STATIC PRESSURE: TEMPERATURE: Option = Automatic with Value Temperature = $InitT [K] END # TEMPERATURE: TURBULENCE INITIAL CONDITIONS: Option = Low Intensity and Eddy Viscosity Ratio END # TURBULENCE INITIAL CONDITIONS: END # INITIAL CONDITIONS: END # INITIALISATION: END # FLOW:Flow Analysis 1 > update FLOW: Flow Analysis 1 &replace SOLVER CONTROL: Turbulence Numerics = $TurbNum ADVECTION SCHEME: Option = $Advect END # ADVECTION SCHEME: CONVERGENCE CONTROL: Length Scale Option = Conservative Maximum Number of Iterations = $nMax Minimum Number of Iterations = $nMin Timescale Control = Auto Timescale Timescale Factor = 1.0 END # CONVERGENCE CONTROL: CONVERGENCE CRITERIA: Residual Target = $Norm Residual Type = $Type END # CONVERGENCE CRITERIA: DYNAMIC MODEL CONTROL: Global Dynamic Model Control = On END # DYNAMIC MODEL CONTROL: END # SOLVER CONTROL: END # FLOW:Flow Analysis 1 > update LIBRARY: CEL: EXPRESSIONS: VF=maxVal(Velocity u)@REGION:OUTLET / massFlowAve(Velocity \ u)@REGION:OUTLET END END END > update FLOW: Flow Analysis 1 &replace OUTPUT CONTROL: MONITOR OBJECTS: MONITOR BALANCES: Option = Full END # MONITOR BALANCES: MONITOR FORCES: Option = Full END # MONITOR FORCES: MONITOR PARTICLES: Option = Full END # MONITOR PARTICLES: MONITOR POINT: Mon1 Expression Value = VF Option = Expression END # MONITOR POINT:Mon1 MONITOR RESIDUALS: Option = Full END # MONITOR RESIDUALS: MONITOR TOTALS: Option = Full END # MONITOR TOTALS: END # MONITOR OBJECTS: RESULTS: File Compression Level = Default Option = Standard END # RESULTS: END # OUTPUT CONTROL: END # FLOW:Flow Analysis 1 > update # ---------------------------------------------------------------------+ > writeCaseFile > update >writeCaseFile filename=$DEF, operation=start solver interactive # operation="write def file" or "start solver batch" > update

Note that in order to access UDF parameters in post-processing, UDM (User Defined Memory) needs to be initialized using GUI path: Parameters & Customization → User Defined Memory → Number of User-Defined Memory Location where the UDM numbering starts from 0. In your UDF, assign your variable to a specific UDM. E.g. F_UDMI(face, thread, 0) = U_txyz; C_UDMI(cell, threat, 0) = vol_heat_source; '0' here refers to the UDM number.

#include "udf.h" DEFINE_PROFILE(U_txyz, thread, position) { /*position: Index of variable say 'U' to be set */ real x[ND_ND]; /* Position vector of nodes */ real xi, yi, zi, r; face_t f; real R = 0.050; /* Radius in [m] */ real t = CURRENT_TIME; /* Current time of simulation */ begin_f_loop(f, b) /* Loops over all faces in thread 'b'*/ { /* C_CENTROID(x,c,t): returns centroid in real x[ND_ND] */ F_CENTROID(x, f, b); /* Gets centroid of face f */ /* x[ND_ND]: Position vector of nodes*/ xi = x[0]; yi = x[1]; zi = x[2]; r = sqrt(xi*xi + yi*yi + zi*zi); F_PROFILE(f, b, position) = 2.0 + 0.5*sin(t/5)*sin(0.31416*r/R); } end_f_loop(f, thread) }

Another example of UDF which can be used to define a linear translation and angular displacement to a cylinder of plate with dynamic mesh motion setting is described below.

/* Rigid body motion of a cylinder: translation and rotations, can be used for a flapping plate if translation is set 0. */ #include "udf.h" /* ----- Define frequency of rotation / flapping in Hz. */ #define f 5.0 /* ----- Define angular velocity in [rad/s]. */ #define omega 2.0*M_PI*f /* ----- Define maximum angualar deflection in [rad] */ #define thetaMax M_PI/180 /* ----- Define linear translation in [m] */ #define xMax 0.01; DEFINE_CG_MOTION(TransRot, dt, cg_vel, cg_omega, time, dtime) { real angVel, linVel; linVel = xMax * omega * cos(omega*time); cg_vel[0] = linVel; cg_vel[1] = 0.0; cg_vel[2] = 0.0; /*cg_omega[0] -> wx, cg_omega[1] -> wy, cg_omega[2] - wz */ /*Axis of rotation is about origin and should be ensured. */ angVel = ThetaMax * omega * sin(omega*time); cg_omega[1] = angVel; cg_omega[2] = 0.0; cg_omega[3] = 0.0; }

- Data Types: Node is a structure data type that stores data associated with a mesh point, face_t is an integer data type that identifies a particular face within a face thread, cell_t is an integer data type that identifies a particular cell within a cell thread, Thread is a structure data type that stores data that is common to the group of cells or faces that it represents.
- NODE_X(node) = X-coordinate of node
- NODE_Y(node) = Y-coordinate of node
- NODE_Z(node) = Z-coordinate of node
- F_NNODES(f, t) = number of nodes in a face
- C_CENTROID(x, c, t) / F_CENTROID(x, f, t): = returns centroid in a 1x3 array real x[ND_ND]
- F_AREA(a, f, t): returns face area (normal) vector in real a[ND_ND]
- C_VOLUME(c, t) = volume of cell
- C_NNODES(c, t) = number of nodes
- C_NFACES(c, t) = number of faces
- C_FACE(c, t, i) = global face index, corresponding to local face index i
- C_R(c, t) / F_R(f, t) = Access cell / boundary face flow variable
**density** - C_P(c, t) / F_P(f, t) = Access cell / boundary and interior face flow variable
**pressure** - C_U(c, t) / F_U(f, t) = Access cell / boundary face flow variable
**x-velocity** - C_V(c, t) / F_V(f, t) = Access cell / boundary face flow variable
**y-velocity** - C_W(c, t) / F_W(f, t) = Access cell / boundary face flow variable
**z-velocity** - C_T(c, t) / F_T(f, t) = Access cell / boundary face flow variable
**temperature** - C_H(c, t) / F_H(f, t) = Access cell / boundary face flow variable
**enthalpy** - Similarly, K, NUT, D, O and YI can be added to C_ to access TKE, turbulent viscosity of Spalart-Allmaras model, TED (ε), specific dissipation rate (ω) and mass fraction respectively
- F_FLUX(f, t) = Access face flow variable
**mass flux** - d = Get_Domain(1) = Get domain using FLUENT utility. Here, domain_id = 1 is an integer for the mixture or fluid domain, but the values for the phase domains can be any integer greater than 1. E.g. porousDomain = Get_Domain(2);
- Derivative macro: C_UDSI(c, t, VORTICITY_Z) = C_DVDX(c,t) - C_DUDY(c,t) = define z-component of vorticity
- C_R_G(c, t), C_P_G(c, t), C_U_G(c, t) ... can be used to get gradient vector: {dp/dx dp/dy dp/dz)
- The M1 suffix can be applied to some of the cell variable macros to access value of the variable at the previous time step (t - Δt). These data may be useful in unsteady simulations. For example, C_T_M1(cell, thread); returns the value of the cell temperature at the previous time step
- Similart, M2 suffix can be applied to some of the cell variable macros to access value of the variable at second previous time step (t -2Δt). These data may be useful in unsteady simulations. For example, C_T_M2(cell, thread); returns the value of the cell temperature at the seond previous time step
- PRINCIPAL_FACE_P(f,t) = test if the face is the principle face - not required for serial operation, applicable to compiled UDFs only
- C_VOLUME = get the cell volume and accumulates it. The end result is total volume of each cell respective mesh.
- C_UDMI(c, t,0) = store result in user-defined memory, location index 0

real Ar1 = 0.0;

begin_f_loop(f, t)

if PRINCIPAL_FACE_P(f, t) {

/* compute area of each face */

F_AREA(area, f, t);

/*compute total face area by summing areas of each face*/

Ar1 = Ar1 + NV_MAG(area);

}

end_f_loop(f,t)

Ar1 = PRF_GRSUM1(Ar1);

Message("Area = %12.4e \n", Ar1);

real Vm1 = 0.0;

begin_C_loop(c, t)

/*compute total volume by summing volumes of each cell*/

Vm1 = Vm1 + C_VOLUME(c, t);

}

end_f_loop(c, t)

Vm1 = PRF_GRSUM1(Vm1);

Message("Volume = %12.4e \n", Vm1);

- To write output to a file from UDF operation / loop: open a file using statement: FILE *fileName; fileName = fopen("udfOutput.txt", "a"); fprintf(fileName,"%g %g\n", x[0], source); fclose(fileName); However, all the fopen / fprintf / fclose commands are incompatible with parallel operation.
- To send output to the console: Message("x = %g source = %g\n", x[0], source);)
- In STAR-CCM+ du/dy = grad($Velocity[0]/grad($Centroid[1]), in FLUENT UDF: C_DUDY(c, t). Similarly, C_DVDX(c, t) can be defined.

for {set i 1} {$i <= $NS} {incr i} {

ic_geo_duplicate_set_fam_and_data point ps$j ps[expr {$j+2}] {} _0

ic_move_geometry point names ps[expr {$j+2}] translate "0 $L 0"

ic_geo_duplicate_set_fam_and_data point ps[expr {$j+1}] ps[expr {$j+3}] {} _0

ic_move_geometry point names ps[expr {$j+3}] translate "0 $L 0"

set j [expr {$i*10+1}]

}

The content on CFDyna.com is being constantly refined and improvised with on-the-job experience, testing, and training. Examples might be simplified to improve insight into the physics and basic understanding. Linked pages, articles, references, and examples are constantly reviewed to reduce errors, but we cannot warrant full correctness of all content.

Copyright © 2017 - All Rights Reserved - CFDyna.com

Template by OS Templates