Boundary Conditions: Last Updated on 24-Feb-2024
The boundary conditions of any problem is used to define the upper and lower limits of the field variables (albeit in absence of any source or sink). These are the operating conditions which govern both the micro- and macro behaviors of these variables. A suitable choice of boundary conditions is as good as a good test set-up! Intuitively, a boundary condition implies that "it is known what happens" on a particular boundary.
For flow and thermal simulation jobs, consultancy, training, automation and longer duration projects, reach me at fb@cfdyna.com. Our bill rate is INR 1,200/hour (INR 9,600/day) with minimum billing of 1 day. All simulations shall be carried out using Open-source Tools including documentations.
There are different (combination) of boundary conditions. For example, in a structural simulation, the number of boundary conditions can be varied to ensure the force- and moment balance of the entire system. This can be achieved by applying boundary condition at just one node or at 6 different nodes! Similarly, in any fluid problem, there must be an entry and an exit for the fluid (as an exception buoyancy-driven flow can be omitted for the time being). This most basic condition is termed as "Inlet" and "Outlet" boundary conditions in CFD parlance, though the choice of "field variables" such as velocity, pressure, temperature, mass flow rate, may vary as per problem set-up.
In any practical application of CFD simulations, the computational domain may comprise of many cells zones (fluid and solid zones) and boundary zones (walls, inlets and outlets). The engineer responsible for pre-processing may not be the one who creates solver file and post-processes the results. The reviewer(s) of the mesh and simulation set-up will certainly be not the engineer who created them. In order to convey the domain information seamlessly, a naming convention should be adopted, it can be a generic system applicable for large number of projects or a specific system for particular simulation set-up. An example is outlined below with following default setting: Newtonian, stationary, adiabatic, smooth boundaries or zones can be named arbitrarily though it is recommended to chose names and identifiers meaningfully.
References - [1]: "Fluid Mechanics, Fourth Edition, Frank M. White", [2]: "HEAT AND MASS TRANSFER by YUNUS A. ÇENGEL from University of Nevada and AFSHIN J. GHAJAR from Oklahoma State University", [3]: "Handbook of Hydraulic Resistance by I. E. Idelchik", [4]: "Chemical and Process Thermodynamics by B. G. Kyle", [5]A Heat Transfer Textbook, Fifth Edition by John H. Lienhard IV and John H. Lienhard V
Conversion Factors as per Reference [1] - Mass and Volume: 1 slug 32.174 = lbm 14.594 kg, 1 lbm = 0.4536 kg, 1 U.S. gal = 0.0037854 m^{3}, 1 U.S. fluid ounce = 2.9574E-5 m^{3}, 1 lbm/ft^{3} = 16.0185 kg/m^{3}, Pressure and Viscosity: 1 lbf/ft^{2} = 47.88 Pa, 1 poise = 1 g/cm.s = 0.1 Pa.s, 1 slug/ft-s = 47.88 Pa.s, Force and Energy: 1 BTU = 252 cal = 1055.056 J, 1 cal = 4.1868 J, 1 Btu/lb-°R = 4186.8 J/kg-K, 1 Btu/h-ft-°R = 1.7307 W/m-K Derived units: 1 BTUH = 0.29307 W, 1 gal/min or GPM = 0.06309 L/s, 1 CFM or ft^{3}/minute = 4.72E-4 m^{3}/s = 0.472 L/s
Validity of Continuum: It is assumed that the underlying assumption of continuum is valid for fluids. The assumptions of continuum breaks down at very low pressures and molecular dynamics takes over. Thus, application of CFD simulations to fluids under vacuum conditions should be taken after calculating Mean Free Path length which is calculate as λ = 1/[π√2] x 1/[d^{2}N_{A}] x R.T/p where d is the diameter [m] of the gas molecule, N_{A} is the Avogadro number, R is the specific gas constant [J/kg-K], T is temperature of gas in [K] and p is absolute pressure in [Pa]. A dimensionless parameter Knudsen number, Kn = λ / L, where λ / L, where λ is the mean free path and L is the characteristic length describes the degree of departure from continuum. Usually when Kn ≥ 0.01, the concept of continuum does not hold good. Beyond this critical range of Knudsen number, the flows are known as (1) slip flow if 0.01 < Kn ≤ 0.1, (2) transition flow if 0.1 ≤ Kn < 10 and (3) free-molecule flow when Kn ≥ 10.
The diameter of an oxygen molecule is 292 picometers or 2.92 angstroms or 2.92 x 10^{-10} [m] and that of nitrogen is 300 picometers = 3.00 x 10^{-10} [m]. The molecular diameter of CO_{2} is larger than that of O_{2}, with a value of 3.34 x 10^{-10} [m]. For nitrogen molecular at 10,000 [Pa] and 0 [°C], the mean free path is approx. 3.37E-05 [m], the value would be 0.0034 [m] at 100 [Pa]. Thus, for a device with characteristics length 0.10 [m] or 100 [mm], Kn would be 0.0034/0.10 = 0.034
The definition of material properties needs to be aligned with the expected variation in fluid properties. For example, if the variation in pressure and temperature is so high that is can cause variation in density of fluid ≥ 10% of free-stream conditions, the fluid should not be modelled as constant density. Same is the case for thermal conductivity and viscosity. The decision to use the constant value of thermal conductivity or temperature-dependent value should be based on expected temperature of the solid walls.
In some applications such as fans, there are significant reduction in pressure values (such as trailing edge and flow separation regions) as well as significant increase in pressure values (near leading edge). Hence, flow simulations dealing with gas and turbo-machine should use ideal gas law as final settings. Constant density can be used initially to get a convergence. The option of incompressible-ideal-gas available in ANSYS FLUENT makes the density a function of temperature only. This option also does not take into account change in density due to variation of pressure near the blade tip and separation regions.
To get a better accuracy, temperature dependent properties should be used. Some correlations available for gases are tabulated below. The Sutherland coefficients for some common gases are as follows. μ = μ_{0} * (273+C)/(T+C)*(T/273)^{1.5} where C is Sutherland coefficient in [K], T is fluid temperature in [K] and μ is viscosity in [Pa.s]. Reference [3].
Fluid | μ_{0} [Pa.s] | C [K] |
Air | 17.12 | 111 |
N2 | 16.60 | 104 |
O2 | 19.20 | 125 |
CO2 | 13.80 | 254 |
H2O-g | 8.93 | 961 |
H2 | 8.40 | 71 |
The variation in dynamic viscosity of air in the temperature range 0 to 100 [°C] is of the order of 25% which is a significant amount if the wall friction contributes more to the pressure drop in the system. Similarly, a 5% reduction in local pressure shall cause 5% reduction in density and hence the velocity or volume flow rate has to increase by 5% to ensure the mass balance.
Specific Heat Capacity - Reference [4]: Cp [kJ/kmol-K] = a + b.T + c.T^{2} + d.T^{3} where T is in [K]. For air, the coefficients are: a = 28.11 [kJ/kmol-K], b = 1.967E-03 [kJ/kmol-K^{2}], c = 4.802E-06 [kJ/kmol-K^{3}], d = -1.966E-09 [kJ/kmol-K^{4}]. Taking molecular weight of air to be 28.96 [kg/kmol], the coefficients to calculate Cp in [J/kg-K] are: a = 970.65 [J/kg-K], b = 0.06792 [J/kg-K^{2}], c = 1658E-04 [J/kg-K^{3}], d = -6.788E-08 [J/kg-K^{4}]. As you can notice, there is > 1% variation in specific heat capacity of air over the range of 100 [°C] and hence it can be assumed a constant value in calculations and simulations. Also note that for steady state simulations, the specific heat capacity and even density is not needed for solids. This can be inferred from the conduction equation for steady state conditions where only thermal conductivity appears as coefficient.
Thermal Conductivity - Reference [4] - for air: k [W/m-K] = a + b.T + c.T^{2} + d.T^{3} where T is in [K]. a = -3.933E-04 [W/m-K], b = 1.018E-04 [W/m-K^{2}], c = -4.857E-08 [W/m-K^{3}], d = 1.521E-11 [W/m-K^{4}]. There is 25% increase in thermal conductivity of air in the temperature range 15 ~ 100 [°C] and hence a constant thermal conductivity may yield lower heat transfer rate. For air at 50 [°C], k = 0.02735 [W/m-K] and at 100 [°C], k = 0.03095 [W/m-K]. In STAR: Field function for temperature dependent thermal conductivity: thCondAir = 1.52e-11 * pow($Temperature,3) - 4.86e-8 * pow($Temperature,2) + 1.02e-4 * $Temperature - 3.93e-4 where T is in [K]. As per the formulat, thCondAir = 0.0262 [W/m-K] at T = 300 [K]. Reference [2] - Soft Rubber: 0.13 [W/m-K], Glass Fiber: 0.043 [W/m-K], Urethane Rigid Foam: 0.026 [W/m-K]
Water: Reference [1] - Suggested curve fits for water in the range 0 ≤ T ≤ 100 [°C]: ρ (kg/m^{3}) ≈ 1000 - 0.0178 ( T[°C] - 4[°C] )^{1.7} ± 0.2%. Viscosity: ln(μ/μ_{0}) ≈ [-1.704 - 5.306.z + 7.003.z^{2}] where z = 273 [K]/T [K] and μ_{0} = 0 1.788 E-3 [kg/(m.s)] ≡ [Pa.s]. Reference [2] - Thermal conductivity of Olive and Peanut Oil: 0.168 [W/m-K], Rubber Natural: 0.28 [W/m-K], Rubber-vulcanized-soft: 0.13 [W/m-K], Rubber-vulcanized-hard: 0.16 [W/m-K]
Emissivities
Reference [2]: "HEAT AND MASS TRANSFER by YUNUS A. ÇENGEL from University of Nevada and AFSHIN J. GHAJAR from Oklahoma State University"
Material | Tempearature [K] | Emissivity |
Aluminium | ||
Polished | 300-900 | 0.04-0.06 |
Commercial sheet | 400 | 0.09 |
Heavily oxidized | 400-800 | 0.20-0.33 |
Anodized | 300 | 0.8 |
Steel | ||
Polished sheet | 300–500 | 0.08-0.14 |
Commercial sheet | 500-1200 | 0.20-0.32 |
Heavily oxidized | 300 | 0.80 |
Stainless Steel | ||
Polished sheet | 300–1000 | 0.17-0.30 |
Lightly oxidized | 600-1000 | 0.30-0.40 |
Heavily oxidized | 600-1000 | 0.70-0.80 |
Non-metals | ||
Teflon | 300–500 | 0.85-0.92 |
Rubber soft | 300–500 | 0.86 |
Paints black | 300 | 0.88 |
Wood | 300 | 0.90 |
Non-Newtonian Fluids: Most of the fluids behave such that shear stress is linear function of strain rate with zero initial threshold value of shear stress. However, there are some material where shear stress is non-linear function of strain rate with non-zero initial threshold value of shear stress (yield point) before which flow cannot start. Various type of such material also known as Bingham plastic and Herschel-Bulkley are described below.
Set Mass Flow Inlet boundary conditions using TUI: /define/b-c/mass-flow-inlet zNAME y y n 0.75 n 50.0 n 10000 n y n n y 2 5
zNAME | Name of the zone |
y | Reference Frame: Absolute? |
y | Mass Flow Specification Method: Mass Flow Rate? |
n | Use Profile for Mass Flow Rate? |
0.75 | Mass Flow Rate in [kg/s] |
n | Use Profile for Total Temperature? |
50 | Total Temperature [in unit selected i.e. K or °C] |
n | Use profile for Supersonic / Initial Gauge Pressure? |
10000 | Supersonic Initial Gauge Pressure (in [Pa]) |
n | Direction Specification Method: Direction Vector? |
y | Direction Specification Method: Normal to Boundary? |
n | Turbulence Specification Method: K and Epsilon? |
n | Turbulence Specification Method: Intensity and Length Scale? |
n | Turbulence Specification Method: Intensity and Viscosity Ratio? |
y | Use profile for reference frame Z-component of rotation axis? |
2 | Turbulent Intensity [%] |
5 | Turbulent Viscosity Ratio |
Some other considerations during application of Inlet B.C. is "Fully Developed Flow" Vs "Developing Flow". For example, if you are a beginner learning tips and trick of CFD by trying to simulate HTC and correlating it with Dittus-Boelter equation, make sure that the flow regime is fully developed. Sometimes, the inlet of the problem set-up is moved upstream the actual location to get the flow a bit developed. Specification of turbulence parameters (turbulent kinetic energy, TKE and turbulent eddy dissipation, TED) should be based on actual measurement of as far as possible. When there are any source of momentum such as centrifugal fan in the computation domain or sharp edges, the overall result gets affected by the turbulence set at the inlet. Followings are the methods to specify turbulence:
Boundary source: Inlet can also have a boundary source define to model heat sources such as solar radiation.
Recirculation Inlet: in many cases there is an internal source of momentum (such as fan) which circulates the air inside a closed domain. Such flow system does not have an inlet or outlet in the sense defined earlier. Such source of momentum can be modeled as fan boundary condition if its performance characteristics is known. However, in case a fixed flow boundary type has to be applied, ANSYS FLUENT has option to use Recirculation Inlet whereas Siemens STAR-CCM+ has option to use Interface of type "Fully Developed Flow".
Walls are required to store a liquid or contain the expansion (mixing) of gases. Since all the fluid flow has to be contained inside walls or at least in a channel, wall B.C. is natural extension into the numerical simulation process. Wall are not only the source of 'turbulence' that gets generated in the flow domain, its surface characteristics becomes important if certain assumptions gets violated.
In any CFD software it is not necessary to create 'named' 2D regions for the walls. This is because any faces of a 3D region which do not explicitly have a 2D region assigned to them, are automatically assigned to the default B.C. 'wall' having 'Adiabatic' condition. In case one wishes to create walls such as "Isothermal / Rotating / Heat Flux Wall", it must be created during the pre-processing.
Typically, there is no flow across the wall boundary conditions. However, in case of permeable or porous walls, flow does occur across the wall. Similarly, in case of suction or blowing (for example transpiration cooling in Gas Turbine Blades), the mass flow rate specifications are required on the wall boundary conditions.
In a stationary domain with rotating wall is valid only if wall motion is entirely tangential. If there is any normal component, rotating wall shall push the adjacent fluid element. Hence, it must be modelled as rotating fluid with stationary walls relative to rotating domain.
The setting for wall boundary conditions in ANSYS FLUENT is shown below: Typical classification of wall B.C. is:Surface Pair | Thermal Conductance [W/m^{2}K] |
Copper - Copper | 10,000 ~ 25,000 |
Aluminium - Aluminium | 2,200 ~ 12,000 |
Ceramic - Metals | 1,500 ~ 8,500 |
Stainless steel - Stainless steel | 2,000 ~ 3,7000 |
Ceramic - Ceramic | 500 ~ 3,000 |
Shell conduction can be used to account for thermal mass in transient thermal analysis problems such as thermal soaking (ramp-up) or thermal cool-down. It can also be used for multiple junctions and allows heat conduction through the junctions. Shell conduction can be applied on boundary walls as well as internal walls. Fluxes at the ends of a shell conducting wall are not included in heat balance reports. These fluxes are accounted for correctly in the ANSYS FLUENT solution, but are not listed in the flux report.
FLUENT creates a shadow wall by adding suffix ':external' to the walls on which shell conduction is switched ON.
Reference: "HEAT AND MASS TRANSFER by YUNUS A. ÇENGEL from University of Nevada and AFSHIN J. GHAJAR from Oklahoma State University"
Natural Convection (Buoyancy-driven Flows)
Boussinesq approximation: ρ_{EFF} = 1 - β(T_{WALL} - T_{REF}). Here: ρ_{EFF}= effective driving density β = coefficient of thermal expanison (equal to 1/T for ideal gas, T in [K]). Boussinesq approximation is only valid for β × (T_{WALL} - T_{REF}) << 1.0 say β × (T_{WALL} - T_{REF}) < 0.10. In other words, the temperature difference should be < 30 [K] for reference temperature of 300 [K].When the Boussinesq approximation is not used, in buoyancy-driven flows, the body force [Fz = (ρ - ρ_{OP})×g] acts in the momentum equation. Here, ρ_{OP} is operating density, default method to compute the operating density in ANSYS FLUENT is by averaging over all cells.
As per FLUENT User's Manual: "The boundary pressures that you input at pressure inlet and outlet boundaries are the redefined pressures as given by equation above. In general you should enter equal pressures, p'_{HS} at the inlet and exit boundaries of your ANSYS FLUENT model if there are no externally-imposed pressure gradients. For example, if you are solving a natural-convection problem with a pressure boundary, it is important to understand that the pressure you are specifying is p'_{HS}. Therefore, you should explicitly specify the operating density rather than use the computed average. The specified value should, however, be representative of the average value. In some cases the specification of an operating density will improve convergence behavior, rather than the actual results. For such cases use the approximate bulk density value as the operating density and be sure that the value you choose is appropriate for the characteristic temperature in the domain."
Strictly speaking, this is not a boundary condition. That is, any numerical simulation can proceed without it. However, this is a great tool to reduce the computational effort and resource if the flow can be envisaged to be symmetrical about a plane or pair of planes. It must be noted that there is a subtle difference between geometrical symmetry and periodicity. Periodic interfaces are treated as if one side of the interface has been translated or rotated to align with the second side of the interface. The periodic type determines the type of transformation (translational or rotational) used to map one side of the interface to the other.
Strictly speaking, this also is not a boundary condition. That is, any numerical simulation can proceed without it. However, this is a great tool to reduce the computational effort and resource if the flow can be envisaged to be symmetrical about a plane or pair of planes. It must be noted that the geometrical symmetry does not guarantee symmetry of the flow. Similarly, cases where micro-structure of flow eddies are being captured such "Large Eddy Simulation" or "DES – Detached Eddy Simulation", symmetry cannot be used owing to inherent 3D nature of the eddies.
The newer versions of program FLUENT and STAR-CCM+ have option to use named expression and field functions for customized boundary conditions. This eliminate need to write an UDF (User Defined Function). For example, following 5 expressions can be used to define inlet mass flow rate (kg/s) as function of static pressure (Pa) at outlet where mass flow rate is a parabolic function of pressure mdot [kg/s] = C_{0} [kg/s] + C_{1} [m.s] × pr + C_{2} [m^2 s^3 kg^-1] × pr^{2}. Here pr = MassFlowAve(StaticPressure, ['Outlet']).
/define/named-expressions add exprName definition "1.25 [m/s]" q: the last entry 'q' is needed to bring the TUI to root position.
In case of field values: /define/named-expressions add exprName definition "MassFlowAve(StaticPressure, ['outlet-1', 'outlet-1']" q - here 'outlet-1' and 'outlet-2' are names of the zone of type outlet.
The boundary conditions can be defined in terms of named expressions: define/b-c/set velocity-inlet inlet () vmag "namedExpr" q where "namedExpr" is the name of the Named Expression with double quotes. vmag is a built-in variable which tells that velocity magnitude is being set or updated. As summarized in the sample journal file above: a zone type can be changed by "define b-c zone-type zName mass-flow-inlet" where zName is the name of the zone. Only the new type of zone (Mass Flow Inlet here) is needed irrespective of the original type of that zone. Once a boundary is set as type m-f-i, use following lines to make additional changes: "define b-c set m-f-i zName () direction-spec n y q" where 'n' is response to "Direction Specification Method: Direction Vector?" and 'y' is response to "Direction Specification Method: Normal to Boundary?".
The heat transfer rate has unit of [W/m^{2}] or [cal/cm^{2}]. However, there is a similar unit [kWhr] - this is measure of energy and 1 [kWhr] = 1,000 [W] x 3,600 [s] = 3,600,000 [J] = 3.6 [MJ]. kWhr is usually used to denote electrical power consumption while dealing with large values. In many cases such as solar radiation the energy flux is represented as [kWhr/m^{2}/day] which will translate into average value of 1000 [W] × 3600 [s] / [24×3600] = 41.67 [W/m^{2}]. The use of kWhr is mainly because the users are charged for energy consumption and not the rate of energy [power] consumption. Thus, a device which consumes 2 [kW] and runs for 15 hours a day will consume 30 [kWhr] = 108 [MJ].
Except shell and tube heat exchangers, for most engineering problems it is impractical and sometimes impossible to correctly model individual fins and tubes of a heat exchanger core. In cross-flow heat exchangers such as automotive radiators, air is called primary fluid (the heat sink) and water or coolant as auxiliary (secondary) fluid (the heat source). There are multiple approaches to model such heat exchangers when combined with surrounding systems and sub-systems.
The heat exchanger models were designed for 'compact' heat exchangers, implying that the primary fluid side flow is unidirectional. The auxiliary fluid is assumed to flow through a large number of parallel tubes, which can optionally double back in a serpentine pattern to create a number of passes. One can independently choose the direction of principal auxiliary fluid flow, the pass-to-pass direction and the direction of external primary fluid flow. Pressure loss is modeled as a momentum sink in the momentum equation and heat transfer is modeled as a heat source in the energy equation.
Lumped Model | Macro Model | Dual Cell Model | |
Accounts for the pressure loss and auxiliary fluid heat rejection | Simple-effectiveness-model | Number-of-Transfer-Units (NTU) model | Number-of-Transfer-Units (NTU) model |
Core modeled as volumetric heat source | Compute auxiliary fluid inlet temperature for a fixed heat rejection or total heat rejection for a fixed auxiliary fluid inlet temperature | NTU method for heat transfer calculations | |
- | Auxiliary flow is modeled as 1D flow | Auxiliary flow modelled on a separate mesh (other than the primary fluid mesh) |
Heat Exchangers Models in STAR-CCM+
Heat Exchangers models are used to replicate the transfer of heat between two fluid streams - hot fluid and cold fluid. Applications include air conditioner evaporators, condensers, charge air coolers, radiators, electric heaters, and electronic devices. There are two options available:
Solver Behaviour | Inlet | Outlet |
Most Robust | Velocity or Mass Flow Rate | Static Pressure |
Somewhat Robust | Total Pressure | Velocity or Mass Flow Rate |
Sensitive of Guess (Initialization) | Total Pressure | Static Pressure |
Unreliable | Static Pressure | Static Pressure |
Not possible (divergence guaranteed) | Any | Total Pressure |
Solver Behaviour | Inlet | Outlet |
Most Robust | Velocity or Mass Flow Rate | Static Pressure |
Somewhat Robust | Velocity or Mass Flow Rate | Outflow or Outlet-vent |
Only for incompressible flows | Velocity Inlet | Outflow |
Not available | Any | Mass Flow Rate |
Not for compressible | Specified Velocity | Any |
A fan is considered to be infinitely thin, and the discontinuous pressure rise across it is specified as a function of the velocity through the fan. The relationship may be a constant, a polynomial - of the form a + b*x^{2} + ... , or piecewise-linear, or piecewise-polynomial function, or a user-defined function. In FLUENT, a zero thickness plane can be used to represent a fan. However, in CFX a volume is required to represent fan as momentum source coefficient. The general momentum source in CFX has unit of [kg/m^{3}-s^{2}]. The source can be linearlized for better convergence and stability which has unit of [kg/m^{3}-s].
Sometimes, the axis of the fan or MRF domain does not pass through origin. If the geometry or mesh was created in some other units (not in meters), scaling the domain in ANSYS FLUENT may change the location of axis. That is the coordinates of a point on axis of the fan may not be same when the mesh is scaled into meter in ANSYS FLUENT. This is bacause scaling of a 3D space needs some reference point: centroid of the geometry or orgin are two ideal options. The X, Y and Z coordinates on axis of a fan can be estimate by calculating Vertex Average of all the nodes on any axi-symmetric surface of the rotating domain. However, this is not a precise method unless the distribution of nodes can be assumed to be uniform on the surface. Another option is to find out the centre of a circle passing though 3 non-collinear points.
Let (x_{0}, y_{0}, z_{0}) is the centre of circle passing though 3 points (x_{1}, y_{1}, z_{1}), (x_{2}, y_{2}, z_{2}) and (x_{3}, y_{3}, z_{3}). From equality of distance (radius): (x_{0} - x_{1})^{2} + (y_{0} - y_{1})^{2} + (z_{0} - z_{1})^{2} = (x_{0} - x_{2})^{2} + (y_{0} - y_{2})^{2} + (z_{0} - z_{2})^{2}. This simplifies to following set of equations. Note that this is valid only when the points are truely 3D that is the circumcircle is not aligned with either of the axes.
2(x_{1}-x_{2}) x_{0} + 2(y_{1}-y_{2}) y_{0} + 2(z_{1}-z_{2}) z_{0} = (x_{1}^{2} - x_{2}^{2}) + (y_{1}^{2} - y_{2}^{2}) + (z_{1}^{2} - z_{2}^{2}) ... [1]2(x_{2}-x_{3}) x_{0} + 2(y_{2}-y_{3}) y_{0} + 2(z_{2}-z_{3}) z_{0} = (x_{2}^{2} - x_{3}^{2}) + (y_{2}^{2} - y_{3}^{2}) + (z_{2}^{2} - z_{3}^{2}) ... [2]
2(x_{3}-x_{1}) x_{0} + 2(y_{3}-y_{1}) y_{0} + 2(z_{3}-z_{1}) z_{0} = (x_{3}^{2} - x_{1}^{2}) + (y_{3}^{2} - y_{1}^{2}) + (z_{3}^{2} - z_{1}^{2}) ... [3]a_{1} x_{0} + b_{1} y_{0} + c_{1} z_{0} = d_{1} ... [4]
a_{2} x_{0} + b_{2} y_{0} + c_{2} z_{0} = d_{2} ... [5]
a_{3} x_{0} + b_{3} y_{0} + c_{3} z_{0} = d_{3} ... [6]
Note that the solution of this set of equations is trivial as there are infinite number of points which can be at equal distance from the three given points. The centre we are looking for is the point which is co-planar with the 3 given points. This requires some special mathematical treatment to find the coordinates of that unique point.For 2D plane when it is aligned to one of the Cartesian planes (XY or YZ or ZX):
First point | X11: | X12: | ||
Second point | X21: | X22: | ||
Third point | X31: | X32: | ||
This is opposite to the fan boundary conditions define above and like fan is also is considered to be infinitely thin membrane, and the discontinuous pressure drop across it is specified as a function of the velocity through the fan. The relationship may be a constant, a polynomial - of the form a + b*x^{2} + ... , or piecewise-linear, or piecewise-polynomial function, or a user-defined function. There is no similar boundary conditions available in CFX.
Examples of porous materials: note that a material which does not allow convection of fluid (e.g. brick, soil, cotton...) may allow diffusion of the fluid (brick, cotton...) such as water. In general, most of the applications of porous domain in CFD simulations refer to convection phenomena such as flow through the radiator core and honeycomb in a catalytic converter.
All the porous media formulation take the form: Δp = -L × (A.v + B.v^{2}) where v is the 'superficial' flow velocity and negative sign refers to the fact that pressure decreases along the flow direction. The 'superficial velocity' is calculated assuming there is no blockage of the flow. L is the thickness of the porous domain in the direction of the flow. Here, A and B are coefficients of viscous and inertial resistances.
Porosity (γ) provides information on the overall pore volume of a porous material and is defined as the ratio of the nonsolid volume (voids) to the total volume of the porous domain. The air permeability and porosity are inter-dependent, an increase in porosity (the free space) should lead to decrease in permeability (resistance to the air flow). v_{SUPERFICIAL} = γ v_{REAL} where v_{REAL} is the actual velocity that can be observed or measured in the pores of the porous media. Since the measurements of velocity in such micro-channels are difficult, a superficial velocity is used which is computed the mean velocity with respect to some reference area. The image below describes a simple arrangement of porous domain in a duct.Most of the time the pressure drop across a porous domain are reported in terms of v_{F} and not v_{REAL}. If one assumes that mathematically there is pressure drop (momentum sink) where porous domain is present though the physical boundaries (solid space) of the porous material is not present, v_{REAL} will become equal to v_{F} and in that case it is (v_{F}) known as superficial velocity.
In FLUENT and CFX, the equation used is: Δp/L = -(μ/α.v + C_{2}.ρ/2.v^{2}) where α is known as 'permeability' and μ is the dynamic viscosity of fluid flowing through the porous domain. Note that FLUENT takes 1/α as input for porous regions. This is a measure of flow resistance and the unit of permeability (α) is [m^{2}]. Other unit of measurement is the darcy [1 darcy = 0.987 μm^{2}], named after the French scientist who discovered the phenomenon. The unit of inertial loss coefficient (C_{2}) is [m^{-1}]. A good knowledge of soil permeability, also termed hydraulic conductivity is needed for estimating the quantity of seepage under dams and dewatering to facilitate underground construction. The degree of pore saturation with water is one of significant factors affecting the soil water permeability. The permeability coefficient has unit of [m/s] and is calculated by Q.L/A/H where Q is the leakage flow rate [m^{3}/s], L is the length of height of porous domain [m], A is cross-section area [m^{2}] and H is applied hydraulic head [m]. The value of soils varies from 10^{-4} to 10^{-10} [m/s]. Clay is the most porous sediment but is the least permeable which acts as an aquitard, impeding the flow of water.
Inertial loss coefficient is known as Quadratic Loss Coefficient in CFX. An alternate formulation in CFX is in terms of Linear Loss Coefficient [kg/m^{3}/s] and Quadratic Loss Coefficients [kg/m^{4}] which are defined as A and B respectively above or P_{v} and P_{i} below.
STAR-CCM+ uses the expression Δp/L = -(P_{v}.v + P_{i}.v^{2}) for a porous domain. Here P_{v} is Viscous Loss Coefficient and has unit [kg/m^{3}/s] and P_{i} is known as Inertial Loss Coefficient having unit [kg/m^{4}]. Note that it is the coefficient of v^{2} which is called Inertial which is contrary to the fact that viscous pressure drop is proportional to v^{2}.
The pressure drop is usually specified as Δp = ζ/2·ρ·v^{2} where ζ is 'equivalent loss coefficient' and is dimensionless. Darcy expressed the pressure gradient in the porous media as v = -[K/μ]·dP/dL where 'K' is the permeability and 'v' is the superficial velocity or the apparent velocity determined by dividing the flow rate by the cross-sectional area across which fluid is flowing. All the programs FLUENT, CFX and STAR-CCM+ uses superficial velocity.Steps to find out viscous and inertial resistances:
When pressure drop is linear function of velocity, it is known as Darcy's Law: Δp = - μ/α × v. This option is implemented explicitly in COMSOL.
In case porous domain is not aligned to any coordinate direction, the direction of unit vector along the flow and across the flow directions can be estimated from following Javascript program. Note that empty field is considered as 0.0. There is no check if a text value is specified in the input fields and the calculator will result in an error.
First point - X1: | |
First point - Y1: | |
First point - Z1: | |
Second point - X2: | |
Second point - Y2: | |
Second point - Z2: | |
Porous Loss Coefficient Calculator
Note that the loss coefficients are assumed to be (and most of the time they are so) parabalic. Hence, only 3 sets of points are required to derive the parabolic equation. It is recommended that the first point is close to no-flow conditions, second point close to the mid-point of operation and the last point towards the highest flow rate conditions. In the table below, P_{i} refers to pressure drop across the porous domain at V_{i}. Air is assumed to be the working fluid. Sutherland law is used to find out viscosity of air at specified temperature.
Thickness of HX | TH | [mm] | |
First point | V1 | [m/s] | |
First point | P1 | [Pa] | |
Second point | V2 | [m/s] | |
Second point | P2 | [Pa] | |
Third point | V3 | [m/s] | |
Third point | P3 | [Pa] | |
Gas Temperature | Tg | [°C] | |
Gas Pressure | Pg | [kPa] | |
Example calculation: x = {0.0 1.0 2.5}, Y = {0 50 75}, m_{1} = 50.0, m_{2} = 16.67, A = -13.33, B = 63.33, C = 0.0. The parabola generated by the curve fit coefficients looks like as shown below. Note that a parabola with local maxima towards higher value of X-axis may have a negative value of A or B.
The formula to derive A, B and C are as follows:
If the porous domain is not inclined the coordiate system of the entire model, a local coordinate system needs to be created to specify viscous and inertial loss coefficients. Most pre-processors such as FLUENT and STAR-CCM+ have option to create a coordinate system based on origin and two points which can be selected graphically using left mouse button. In case it is to be specified using a direction vector such as in macros or scripts, following transformation matrix should be used.
Note that the transformation matrix is for rotation about z-axis. Thus, a vector in Laboratory Reference Frame (1.0, 0.0, 0.0) shall become (cosθ, -sinθ, 0.0) when rotated by angle θ about z-axis.
Porosity provides information on the overall pore volume of a porous material and is defined as the ratio of the nonsolid volume (voids) to the total volume of the nonwoven fabric. The air permeability and porosity of fabric duct are inter-dependent. An increase in porosity (the free space) should lead to decrease in permeability (resistance to the air flow). From flow point of view, a volume based definition of porosity may not provide a direct relation to the the flow velocity which depends on flow area perpendicular to the flow direction.
There is decrease in atmospheric pressure and temperature with altitude as compared to height above sea level. Why sea level is considered as reference datum? This is because the lquids maintain uniform level and any point anywhere in the sea is expected to be same radial distance from the centre of the Earth.
p [Pa] = 101325 * (1 - 2.25577E-05 × H)^{5.2559} where altitude H is in [m].
T [K] = 288.15 - 0.0065 × H. You may use the following calculator to estimate ambient pressure and tempearture at higher altitudes. There is option to chose altitude in [m] or [ft]. However, the outputs are in SI units.
Altitude, H | |
Unit | |
In situations with multi-component flows (such as leakage of fuel or refrigerant) where diffusion dominates the correct specification of binary diffusion coefficient is very important. Following table specifies value at 1 [atm] and 300 [K]. Reference: Analytic Combustion by Anil W. Date.
Pair | D_{ab} [m^{2}/s] | Pair | D_{ab} [m^{2}/s] |
H2O - air | 24.0E-6 | C6H6 - air | 8.00E-6 |
CO - air | 19.0E-6 | C8H18 - air | 5.00E-6 |
CO2 - air | 14.0E-6 | C8H16 - air | 7.10E-6 |
H2 - air | 78.0E-6 | C10H22 - air | 6.00E-6 |
O2 - air | 19.0E-6 | O2 - H2 | 70.0E-6 |
SO2 - air | 13.0E-6 | CO2 - N2 | 11.0E-6 |
NH3 - air | 28.0E-6 | CO2 - H2 | 55.0E-6 |
CH3OH - air | 14.0E-6 | C6H14 - N2 | 8.00E-6 |
C2H5OH - air | 11.0E-6 | C8H18 - N2 | 7.00E-6 |
CH4 - air | 16.0E-6 | C10H22 - N2 | 6.40E-6 |
Here:
Dupuit-Thiem equation (based on Darcy Law in cylindrical coordinate system) is a widely used formula to estimate pressure drop across the wall for a known (oil extraction) flow rate in a circular reservoir that has a constant pressure at its outer boundary.
Dust Accumulation in Air Filters: There are many application of air filters such as automotive air cleaners. Dust Holding Capacity (DHC) is one of the key parameters of such filters. The filter are orthotropic porous media where the porous loss coefficients are different along the 3 directions. However, any CFD simulations to deal with dust accumulation will be a transient simulation where the behaviour of porous domain will change depending upon duct collection level and spatial distribution. This is because filter may not collect dust uniformly and hence permeability will change non-uniformly. For most practical applications, change in pressure drop can be assumed to be a linear function of duct loading (the amount of dust trapped in filters). How does one model the trapping of dust particles in the pores of the filter? Neither the filter pores nor the diameters of the particles are uniform in size and shape!
The rate of convergence slows a porous region is defined such that pressure drop is relatively large in the flow direction (e.g. the permeability is low or the inertial factor is large). This slow convergence can occur because the porous media pressure drop appears as a momentum source term yielding a loss of diagonal dominance in the matrix of equations solved. The best remedy for poor convergence of a problem involving a porous medium is to supply a good initial guess for the pressure drop across the medium. You can supply this guess by patching a value for the pressure in the fluid cells upstream and/or downstream of the medium. It is important to recall, when patching the pressure, that the pressures you input should be defined as the gauge pressures used by the solver (i.e. relative to the operating pressure defined in the simulation).
Another possible way to deal with poor convergence is to temporarily disable the porous media model and obtain an initial flow field without the effect of the porous region. Once an initial solution is obtained, or the calculation is proceeding steadily to convergence, enable the porous media model and continue the calculation with the porous region included. (This method is not recommended for porous media with high resistance.)
Simulations involving highly anisotropic porous media may, at times, pose convergence troubles. This issue can be addressed limiting the anisotropy of the porous media coefficients to two (10^{2}) or three (10^{3}) orders of magnitude. Even if the medium's resistance in one direction is infinite, it is not needed to set the resistance in that direction to be greater than 1000 times the resistance in the primary flow direction.
In the fluid mechanics of porous media, tortuosity is the ratio of the length of a streamline to that of the straight-line distance between those points. A measure of deviation from a straight line. It is the ratio of the actual distance traveled between two points, including any curves encountered, divided by the straight line distance. Tortuosity is used by drillers to describe wellbore trajectory, by log analysts to describe electrical current flow through rock and by geologists to describe pore systems in rock and the meander of rivers.
A related concept is fractal which is used to describe the effective length of rivers and used even for trading in stock markets.
[2]: A method for aerodynamic drop break-up that has been widely implemented is known as the Taylor Analogy Break-up model or TAB. The TAB model approximates a drop as a damped oscillator system and tracks a deformation parameter. When the deformation parameter is sufficiently large, break-up is assumed.
Surface tension and angle of contact with water and some of the non-metals are tabulated below as per www.accudynetest.com.
Surface Tension of Water with | Abbreviation | [N/m] | [°] |
Silicone Oxide | Glass | 0.0725 | ~ 0 |
Poly-Vinyl-Chloride | PVC | 0.0379 | 85.6 |
Poly-Tetra-Fluoro-Ethylene | PTFE | 0.0194 | 109 |
Poly-Amide-6-6 | Nylon-66 | 0.0422 | 68.3 |
Poly-Methyle-Meth-Acrylate | PMMA (Acrylic, Plexiglass) | 0.0375 | 70.9 |
Poly-Ethylene-Terephthalate | PET | 0.0390 | 72.5 |
Poly-Carbonate | PC | 0.0440 | 82.0 |
Acrylonitril Butadine Styrene | ABS | 0.0385 | 80.9 |
Poly-Ethylene | PE | 0.0316 | 96.0 |
Poly-Propylene | PP | 0.0305 | 102 |
One of the applications of capillary effect combined with fully-developed laminar (Poiseuille) flow is pipetting. Pipetting process is aspiration of a pre-determined volume of liquid by creating a vacuum above a tapered capillary tube (known as tip). The pressure in the pipette chamber during the process is in a dynamic equilibrium and is affected by the ambient pressure, viscosity, surface tension and density of the liquid, and the speed of the piston movement. The suction (aspiration) of liquid in pipette tips normally undergo following 4 phases:
Regime-I | Regime-II | Regime-III | Regime-IV |
Initial boundary effects important | Viscous effect negligible | Poiseuille flow: inertia effect negligible | Late viscous regime |
Surface tension forces dominant | Capillary rise resisted by fluid inertia | Lucas-Washburn law applies | Fluid rise approaches steady state height |
Capillary rise: z ∝ t^{2} | Capillary rise: z ∝ t | Capillary rise: z ∝ t^{0.5} | Capillary rise: z ∝ e^{-t} |
The format of the standard transient profile file:
( (inletvel transdata 3 0) (time 1 2 3 ) (invel 6 5 4 ) )Here '3' denotes the number of data points (number of rows in tabular format). '0' indicated that data is not time-periodic, use '1' for a time-periodic dataset. 'invel' shall appear as drop-down in boundary condition panel. Note:
The data set used above can also be specified as described below:
( (inletvel transdata 3 0) (time 1 2 3 ) (invel 6 5 4 ) )
Tabular Transient Profiles
The format of the tabular transient profile file is as follows. Use TUI command "file read-transient-table fileTr.dat" to upload the data into the solver. The options "inletvel time" and "inletvel invel" shall appear in the drop-down menu of boundary condition panels.inletvel transdata 3 0 time invel 1 6 2 5 3 4
Data Import into ANSYS FLUENT:
Interpolation file format: ANSYS FLUENT
2 --the interpolation file version 3 --dimension (2 or 3) 25000 --total number of points 3 --total number of fields (velocity, pressure...) included x-velocity --field variable to be interpolated pressure --field variable to be interpolated -0.100000 --x-coordinate of first node in the list -0.200000 --y-coordinate of first node in the list -0.300000 --z-coordinate of first node in the list ... ... ... 1.000000 --x-coordinate of last node in the list 2.000000 --y-coordinate of last node in the list 3.000000 --z-coordinate of last node in the list 0.500000 --x-velocity at first node in the list 0.100000 --x-velocity at second node in the list 0.200000 --x-velocity at third node in the list ... ... ... 5.000000 --x-velocity at last node in the list 100.0000 --pressure at first node in the list 200.0000 --pressure at second node in the list 300.0000 --pressure at third node in the list ... ... ... 500.0000 --pressure at last node in the list
At the end is a list of the field values at all the points in the same order as their names. The number of coordinate and field points (x-velocity and pressure in this case) should match the number given in line 3 (25000 in this case). Note that it cannot interpolate heat generation rate (W/m^{3}), though internal energy, enthalphy, heat transfer coefficient... can be interpolated. In other words, no volumetric source/sink term can be interpolated in the fluid region. On the boundaries, all field variables are available for interpolation.
Data Export into CSV Format: ANSYS FLUENT:
The header in CSV files for FLUENT and STAR-CCM+ uses different variable names. In STAR-CCM+ variables has to be specified with appropriate units: e.g. "Absolute Pressure (Pa)", "Velocity Magnitude (m/s)", "Velocity[i] (m/s)", "Velocity[j] (m/s)", "Velocity[k] (m/s)", "X (m)", "Y (m)", "Z (m)"... Note the space between variables and unit. The x-component of velocity is accessed by Velocity[i].
Data Mapper: STAR-CCM+:
Using a user field function, the table data which provides the variations of velocity with time can be imported and then defined as follow:
interpolateTable(@Table("VelData"), "Time", LINEAR, "VelMPS", $Time)
Here the field function type should be vector. For boundaries such as inlet, in the “Physics Values” option of the boundary in the “Method” field, select Field Function, and then right below this, select the field function defined earlier.
The name of the table is VelData, Time is the name for the first column in the .csv file, LINEAR (other option available is SPLINE) is the interpolation method creating values between the discrete time values specified in the table, VelMPS is the name of the second column in the table (corresponding to the velocity magnitude), $Time specifies that the velocity is function of time. In case velocity is function of spatial coordinates say y-axis: $${Position}[1] can be used in place of $Time to specify that the position coordinate is aligned with the y-axis (1 = y-direction, use '0' for x-direction and 2 for z-direction).
License Manager: LMTOOLS
Most of the CFD programs be it ANSYS FLUENT or STAR-CCM+ use LMTools to manage the licenses at their customers. The executable utility is installed in the same folder where main program such as ANSYS FLUENT is installed. The syntax to use LMSTAT is: lmstat [-a] [-c license_file_list] [-f [feature]] [-i [feature]] [-s[server]] [-S [vendor]] [-t timeout_value]. The purpose of each option is:Finding Names of Feature: to use -f option, one need to know the exact names of feasures. This can be found in license.dat file or check the output from -a options. lmutil.exe lmstat -s licServers -f "feature_name"
For scripting and partial automation of pre-processing, solver and post-processing activities, refer to this section of scripting page.
For dealing with file handling in Linux and Windows, refer to short summary of shell scripting.
No. | Checkpoint | Record [Y / N / NA] |
01 | Is the CAD data available in STEP format? | |
02 | Are the boundary condition values at inlets available: mass or volume flow rate, temperature, pressure? | |
03 | Are the boundary conditions at outlets available: mass flow rate, pressure? | |
04 | Are the boundary conditions such as surface heat sources at walls available? | |
05 | Are the volumetric heat generation rates available? | |
06 | If applicable, is the fan performance curve available? | |
07 | If applicable, is the pressure drop vs velocity curve available for porous domain? | |
08 | Are the material type and grade for all solids avaiable? | |
09 | Are the values of ambient conditions (temperature and pressure) available? | |
10 | Are the thermal conductivity values of non-standard materials available? |
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