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

Automation

Scripts, Journals, UDF, Field Functions

A collection of scripts, journals and macros in CFD simulations to automate some tasks as well as enhance simulation capabilities. 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 properties and interface/periodicity information can be written out in a log file for self and peer review.

Definitions and concepts
  • 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.

Programs and Scripting/Macro Language
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;
Note:
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 character. 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 access 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
In STAR-CCM, it is achieved as:
  • Magnitude: mag($$Velocity)
  • X-component: $$Velocity[1]
  • Y-Component: $$Velocity[1]...
  • The directions X,Y,Z are represented by [0],[1],[2] respectively.
An example of a journal script which can be used in FLUENT to set-up a solver is as follows. This journal just needs to be edited for user's need. Note that there should be no empty line till the end of file. An empty line is assumed to be end of input. Use the remark character ; instead. The journal can be further parametrized by defining parameter using syntax (define X 0.7). Similarly, mathematical operation can be performed a = a + 0.1 can be written as (define a (+ a 0.1)).
  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 orthotropic thermal conductivity which is otherwise not possible through GUI.
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

EXAMPLE OF UDF
Only velocities in Cartesian coordinates (Ux,Uy,Uz) are accessible through the UDF macros. Radial, tangential & axial velocities (Ur, Ut, Ua) within a fluid zone will need to be computed using the UDF itself. The UDF below can be used to define inlet velocity as a function of radius and time. It is assumed that the centre of inlet face is at origin.
  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 angular 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;
}

Summary of UDF commands
Exerpts from user manual: Think of a ‘zone’ as being the same as a ‘thread’ data structure when programming UDFs. Thus: node thread = grouping of nodes, face thread = grouping of faces, cell thread = grouping of cells. A domain is a data structure in ANSYS FLUENT that is used to store information about a collection of node, face threads, and cell threads in a mesh.
UDF thread domain hierarchy
  • 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
  • C_K_L(c, t) = Access cell material property "thermal conductivity". 'L', 'T' and 'EFF' suffix is used for 'Laminar', 'Turbulent' and 'Effective' properties. For example, C_MU_L, C_MU_T and C_MU_EFF refers to 'Laminar', 'Turbulent' and 'Effective' viscosities respectively.
  • 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 previous to 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 of respective mesh.
  • C_UDMI(c, t,0) = store result in user-defined memory, location index 0. The user-defined memory field gets saved to the data file and can be used to generate contours and any other post-processing activity.
  • boolean FLUID_THREAD_P(thread): Check if 'thread' is of type FLUID

More examples of FLUENT UDF
UDF for Temperature Dependent Viscosity
DEFINE_PROPERTY(visc_T, c, Tp)		{
  real mu; real a = C_T(c,t);  
  mu = 2.414e-05 * pow(10, 247.8/(Tp - 140));
  return mu;						
}
Compute area of a face zone:
 #include "udf.h"
 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);
Compute volume of a cell zone:
 #include "udf.h"
 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);
UDF: change time step value in Transient Simulation This example, UDF needs to operate on a particular thread in a domain (instead of looping over all threads) and the DEFINE macro DEFINE_DELTAT used in UDF does not have the thread pointer passed to it from the solver. Hence, Lookup_Thread is required in the UDF to get the desired thread pointer.
 #include "udf.h"
 DEFINE_DELTAT(timeStep, d) {
 real newTime = CURRENT_TIME; real oldT; real minT = 0.0; real maxT = 0.0;
 Thread *t; cell_t c; d = Get_Domain(1);

 int zID = 1; Thread *t = Lookup_Thread(d, zID);
 begin_f_loop(f, t) { /* Loop over all face elements*/
   oldT = F_T_M1(f, t); /* Get face temperature at previous time-step */
   if (oldT < minT || minT == 0.0) minT = oldT;
   if (oldT > maxT || maxT == 0.0) maxT = oldT;
 }
 end_f_loop(f, t)

 if(maxT < 100.0)
  timeStep = 0.5;
 else
  timeStep = 0.1;

 return timeStep;
 }

Density as function of temperature:

#include "udf.h"
DEFINE_PROPERTY(rho_T, c, Tp)		{
  real rho;
  /* Get temperature of the cell in K and convert into C */
  real Tp = C_T(c,t) - 273.11; 
  real a0 =  999.8396;
  real a1 =  18.22494;
  real a2 = -7.92221e-03;
  real a3 = -5.54485e-05;
  real a4 =  1.49756e-07;
  real a5 = -3.93295e-10;
  real b =   1.81597e-02;

  rho = a0 + a1*Tp + a2*Tp*Tp + a3*pow(Tp, 3) + a4*pow(Tp, 4) + a5*pow(Tp, 5);
  rho = rho / (1 + b*Tp);
  return rho;						
}

Tips and Tricks on UDF / UFF
  • 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(cell, thread). Similarly, C_DVDX(cell, thread) can be defined.
  • Similarly in FLUENT, C_T_G(cell, thread)[0] returns the x-component of the cell temperature gradient vector.

ICEM CFD: Tck/Tk Examples
For loop to copy a point entity:
  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}]
  }


Macro in STAR CCM+
The first step in writing a JAVA macro for STAR CCM+ is to import the relevant packages. For example:
package macro;  - similar to #include "udf.h"
import java.util.*;  - similar to header files in C
import star.turbulence.*; - import turbulence model data
import star.common.*;
import star.material.*;
import star.base.neo.*;
import star.vis.*;
import star.flow.*;
import star.energy.*;
import star.coupledflow.*;

// defineSim is the name of macro and the file name should be defineSim.java.
public class defineSim extends StarMacro {

  public void execute() {
    execute0();
  }
  
  //Get active simulation data
  Simulation getSim = getActiveSimulation();
  
  //Get domain named as FLUID and store it as 'cfd' - similar to Get_Domain in FLUENT UDF
  Region cfd = getSim.getRegionManager().getRegion("FLUID");
  
  //Set viscous model
  PhysicsContinuum pC0 = ((PhysicsContinuum) cfd.getContinuumManager().getContinuum("Physics 1"));
    pC0.enable(SteadyModel.class);
    pC0.enable(SingleComponentGasModel.class);
    pC0.enable(CoupledFlowModel.class);
    pC0.enable(IdealGasModel.class);
    pC0.enable(CoupledEnergyModel.class);
    pC0.enable(TurbulentModel.class);
    pC0.enable(RansTurbulenceModel.class);
    pC0.enable(KEpsilonTurbulence.class);
    pC0.enable(RkeTwoLayerTurbModel.class);
    pC0.enable(KeTwoLayerAllYplusWallTreatment.class);
	  
  //Get boundary named INLET and velocity specification CLASS
  Boundary b1 = cfd.getBoundaryManager().getBoundary("INLET");
  VelocityProfile vP1 = b1.getValues().get(VelocityProfile.class);
  
    //Specify velocity normal to boundary with specified "MAGNITUDE and DIRECTION"
    //Note the word scalar in ConstantScalarProfileMethod.class
    b1.getConditions().get(InletVelocityOption.class).setSelected(InletVelocityOption.MAGNITUDE_DIRECTION);
    vP1.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(5.0);
	
    //Inlet velocity by its COMPONENTS, note 'vector' in ConstantVectorProfileMethod.class
    //b1.getConditions().get(InletVelocityOption.class).setSelected(InletVelocityOption.COMPONENTS);
    //vP1.getMethod(ConstantVectorProfileMethod.class).getQuantity().setComponents(5.0, 0.0, 0.0);
  
    //Set turbulence parameters - TURBULENT INTENSITY and VISCOSITY RATIO at INLET boundary
    TurbulenceIntensityProfile TI = b1.getValues().get(TurbulenceIntensityProfile.class);
    TI.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(0.02);

    TurbulentViscosityRatioProfile TVR = b1.getValues().get(TurbulentViscosityRatioProfile.class);
    TVR.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(5.0);
	
    //Specify fluid temperature in [K] at INLET
    StaticTemperatureProfile Tin = b1.getValues().get(StaticTemperatureProfile.class);
    Tin.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(323.0);
	
  //Get boundary named OUTLET and pressure boundary CLASS
  Boundary b2 = cfd.getBoundaryManager().getBoundary("OUTLET");
  StaticPressureProfile sP0 = b2.getValues().get(StaticPressureProfile.class);
    
    //Specify static pressure at OUTLET boundary
    b2.setBoundaryType(PressureBoundary.class);
    sP0.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(0.0);
  
    //Specify back flow turbulence parameters at OUTLET boundary
    TurbulenceIntensityProfile TI2 = b2.getValues().get(TurbulenceIntensityProfile.class);
    TI2.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(0.01);

    TurbulentViscosityRatioProfile TVR2 = b2.getValues().get(TurbulentViscosityRatioProfile.class);
    TVR2.getMethod(ConstantScalarProfileMethod.class).getQuantity().setValue(2.0);
	
    //Other options for reverse flow specifications
    b2.getConditions().get(BackflowDirectionOption.class).setSelected(BackflowDirectionOption.EXTRAPOLATED);
    b2.getConditions().get(BackflowDirectionOption.class).setSelected(BackflowDirectionOption.BOUNDARY_NORMAL);
    b2.getConditions().get(ReversedFlowPressureOption.class).setSelected(ReversedFlowPressureOption.ENVIRONMENTAL);
    b2.getConditions().get(ReversedFlowPressureOption.class).setSelected(ReversedFlowPressureOption.STATIC);
    b2.getConditions().get(ReferenceFrameOption.class).setSelected(ReferenceFrameOption.LOCAL_FRAME);
    b2.getConditions().get(ReferenceFrameOption.class).setSelected(ReferenceFrameOption.LAB_FRAME);
    b2.getConditions().get(KeTurbSpecOption.class).setSelected(KeTurbSpecOption.INTENSITY_LENGTH_SCALE);
    b2.getConditions().get(KeTurbSpecOption.class).setSelected(KeTurbSpecOption.INTENSITY_VISCOSITY_RATIO);
	
  //Save SIM file by specifying full path - note double backslashes
  getSim.saveState(resolvePath("C:\\STAR_CCM\\PipeFlow.sim"));
  
  //Close the simulation scene
  getSim.close(true);
}

Contact us
Disclaimers and Policies

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.