PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
Loading...
Searching...
No Matches
Functions
wallfunction.h File Reference
#include <petscpf.h>
#include <petscdmswarm.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <petsctime.h>
#include <petscsys.h>
#include <petscdmcomposite.h>
#include <petscsystypes.h>
#include "variables.h"
#include "ParticleSwarm.h"
#include "walkingsearch.h"
#include "grid.h"
#include "logging.h"
#include "io.h"
#include "interpolation.h"
#include "ParticleMotion.h"
#include "Boundaries.h"
Include dependency graph for wallfunction.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

void noslip (UserCtx *user, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, double normal_x, double normal_y, double normal_z)
 Applies no-slip wall boundary condition with linear interpolation.
 
void freeslip (UserCtx *user, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, double normal_x, double normal_y, double normal_z)
 Applies free-slip wall boundary condition.
 
double E_coeff (double friction_velocity, double roughness_height, double kinematic_viscosity)
 Computes roughness-modified log-law coefficient E.
 
double u_hydset_roughness (double kinematic_viscosity, double wall_distance, double friction_velocity, double roughness_height)
 Computes velocity from log-law for rough walls.
 
double f_hydset (double kinematic_viscosity, double known_velocity, double wall_distance, double friction_velocity_guess, double roughness_height)
 Residual function for friction velocity equation (log-law with roughness)
 
double df_hydset (double kinematic_viscosity, double known_velocity, double wall_distance, double friction_velocity_guess, double roughness_height)
 Numerical derivative of residual function.
 
double find_utau_hydset (double kinematic_viscosity, double known_velocity, double wall_distance, double initial_guess, double roughness_height)
 Solves for friction velocity using Newton-Raphson iteration.
 
double nu_t (double yplus)
 Computes turbulent eddy viscosity ratio (ν_t / ν)
 
double integrate_1 (double kinematic_viscosity, double wall_distance, double friction_velocity, int integration_mode)
 Integrates eddy viscosity profile from wall to distance y.
 
double taw (double kinematic_viscosity, double friction_velocity, double wall_distance, double velocity, double pressure_gradient_tangent)
 Computes wall shear stress with pressure gradient effects.
 
double u_Cabot (double kinematic_viscosity, double wall_distance, double friction_velocity, double pressure_gradient_tangent, double wall_shear_stress)
 Computes velocity using Cabot wall function.
 
double f_Cabot (double kinematic_viscosity, double velocity, double wall_distance, double friction_velocity_guess, double pressure_gradient_tangent, double pressure_gradient_normal)
 Residual function for Cabot wall function.
 
double df_Cabot (double kinematic_viscosity, double velocity, double wall_distance, double friction_velocity_guess, double pressure_gradient_tangent, double pressure_gradient_normal)
 Numerical derivative for Cabot wall function.
 
void find_utau_Cabot (double kinematic_viscosity, double velocity, double wall_distance, double initial_guess, double pressure_gradient_tangent, double pressure_gradient_normal, double *friction_velocity, double *wall_shear_velocity, double *wall_shear_normal)
 Solves for friction velocity using Cabot wall function.
 
double u_Werner (double kinematic_viscosity, double wall_distance, double friction_velocity)
 Computes velocity using Werner-Wengle wall function.
 
double f_Werner (double kinematic_viscosity, double velocity, double wall_distance, double friction_velocity)
 Residual function for Werner-Wengle iteration.
 
double df_Werner (double kinematic_viscosity, double velocity, double wall_distance, double friction_velocity)
 Numerical derivative for Werner-Wengle iteration.
 
double find_utau_Werner (double kinematic_viscosity, double velocity, double wall_distance, double initial_guess)
 Solves for friction velocity using Werner-Wengle wall function.
 
double u_loglaw (double wall_distance, double friction_velocity, double roughness_length)
 Computes velocity using simple log-law (smooth wall with roughness offset)
 
double find_utau_loglaw (double velocity, double wall_distance, double roughness_length)
 Solves for friction velocity using simple log-law (explicit formula)
 
void wall_function (UserCtx *user, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, PetscReal *friction_velocity, double normal_x, double normal_y, double normal_z)
 Applies standard wall function with Werner-Wengle model.
 
void wall_function_loglaw (UserCtx *user, double roughness_height, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, PetscReal *friction_velocity, double normal_x, double normal_y, double normal_z)
 Applies log-law wall function with roughness correction.
 
void wall_function_Cabot (UserCtx *user, double roughness_height, double distance_reference, double distance_boundary, Cmpnts velocity_wall, Cmpnts velocity_reference, Cmpnts *velocity_boundary, PetscReal *friction_velocity, double normal_x, double normal_y, double normal_z, double pressure_gradient_x, double pressure_gradient_y, double pressure_gradient_z, int iteration_count)
 Applies Cabot non-equilibrium wall function with pressure gradients.
 

Function Documentation

◆ noslip()

void noslip ( UserCtx user,
double  distance_reference,
double  distance_boundary,
Cmpnts  velocity_wall,
Cmpnts  velocity_reference,
Cmpnts velocity_boundary,
double  normal_x,
double  normal_y,
double  normal_z 
)

Applies no-slip wall boundary condition with linear interpolation.

This function enforces a no-slip boundary condition (zero velocity at the wall) by linearly interpolating between the wall velocity (typically zero) and the velocity at a reference point in the flow.

MATHEMATICAL FORMULATION: For a point at distance sb from the wall, with a reference velocity Uc at distance sc from the wall: U_boundary = U_wall + (U_reference - U_wall) * (sb / sc)

PHYSICAL INTERPRETATION: This provides a first-order approximation assuming linear velocity variation in the near-wall region, which is valid in the viscous sublayer.

Parameters
[in]userSimulation context (unused but required for interface)
[in]distance_referenceWall-normal distance to reference point (sc)
[in]distance_boundaryWall-normal distance to boundary point (sb)
[in]velocity_wallVelocity at the wall (Ua), typically zero
[in]velocity_referenceVelocity at reference point (Uc)
[out]velocity_boundaryComputed velocity at boundary point (Ub)
[in]normal_xX-component of wall normal vector
[in]normal_yY-component of wall normal vector
[in]normal_zZ-component of wall normal vector
Note
For moving walls, velocity_wall would be non-zero
The normal vector should point INTO the fluid domain

Applies no-slip wall boundary condition with linear interpolation.

Local to this translation unit.

Definition at line 68 of file wallfunction.c.

73{
74 (void)user;
75 (void)normal_x;
76 (void)normal_y;
77 (void)normal_z;
78 // Compute velocity difference between reference point and wall
79 double delta_u = velocity_reference.x - velocity_wall.x;
80 double delta_v = velocity_reference.y - velocity_wall.y;
81 double delta_w = velocity_reference.z - velocity_wall.z;
82
83 // Linear interpolation factor
84 double interpolation_factor = distance_boundary / distance_reference;
85
86 // Apply linear interpolation
87 (*velocity_boundary).x = interpolation_factor * delta_u;
88 (*velocity_boundary).y = interpolation_factor * delta_v;
89 (*velocity_boundary).z = interpolation_factor * delta_w;
90
91 // Add wall velocity (shift to absolute frame)
92 (*velocity_boundary).x += velocity_wall.x;
93 (*velocity_boundary).y += velocity_wall.y;
94 (*velocity_boundary).z += velocity_wall.z;
95}
PetscScalar x
Definition variables.h:103
PetscScalar z
Definition variables.h:103
PetscScalar y
Definition variables.h:103
Here is the caller graph for this function:

◆ freeslip()

void freeslip ( UserCtx user,
double  distance_reference,
double  distance_boundary,
Cmpnts  velocity_wall,
Cmpnts  velocity_reference,
Cmpnts velocity_boundary,
double  normal_x,
double  normal_y,
double  normal_z 
)

Applies free-slip wall boundary condition.

Free-slip conditions allow tangential flow but enforce zero normal velocity. This is appropriate for inviscid walls or symmetry planes where there is no shear stress but flow cannot penetrate the boundary.

DECOMPOSITION: Velocity is decomposed into normal and tangential components: U = U_n * n + U_t where U_n = U · n (normal component) U_t = U - U_n * n (tangential component)

BOUNDARY CONDITIONS:

  • Normal component: Interpolated from interior (∂U_n/∂n ≠ 0)
  • Tangential component: Extrapolated from interior (∂U_t/∂n = 0)
Parameters
[in]userSimulation context
[in]distance_referenceWall-normal distance to reference point
[in]distance_boundaryWall-normal distance to boundary point
[in]velocity_wallVelocity at the wall
[in]velocity_referenceVelocity at reference point
[out]velocity_boundaryComputed velocity at boundary point
[in]normal_xX-component of wall normal vector
[in]normal_yY-component of wall normal vector
[in]normal_zZ-component of wall normal vector

Applies free-slip wall boundary condition.

Local to this translation unit.

Definition at line 101 of file wallfunction.c.

106{
107 (void)user;
108 // Extract normal components of velocities
109 double wall_normal_velocity = velocity_wall.x * normal_x +
110 velocity_wall.y * normal_y +
111 velocity_wall.z * normal_z;
112
113 double reference_normal_velocity = velocity_reference.x * normal_x +
114 velocity_reference.y * normal_y +
115 velocity_reference.z * normal_z;
116
117 // Interpolate normal velocity component
118 double boundary_normal_velocity = wall_normal_velocity +
119 (reference_normal_velocity - wall_normal_velocity) * (distance_boundary / distance_reference);
120
121 // Extract tangential velocity components (extrapolated from reference point)
122 double tangential_u = velocity_reference.x - reference_normal_velocity * normal_x;
123 double tangential_v = velocity_reference.y - reference_normal_velocity * normal_y;
124 double tangential_w = velocity_reference.z - reference_normal_velocity * normal_z;
125
126 // Reconstruct total velocity: U = U_t + U_n * n
127 (*velocity_boundary).x = tangential_u + boundary_normal_velocity * normal_x;
128 (*velocity_boundary).y = tangential_v + boundary_normal_velocity * normal_y;
129 (*velocity_boundary).z = tangential_w + boundary_normal_velocity * normal_z;
130}
Here is the caller graph for this function:

◆ E_coeff()

double E_coeff ( double  friction_velocity,
double  roughness_height,
double  kinematic_viscosity 
)

Computes roughness-modified log-law coefficient E.

The coefficient E accounts for wall roughness effects on the log-law: u+ = (1/κ) ln(E * y+)

ROUGHNESS REGIMES:

  1. Hydraulically smooth (ks+ < 2.25): E = exp(κ*B), no roughness effect
  2. Transitional (2.25 < ks+ < 90): Smooth interpolation between regimes
  3. Fully rough (ks+ > 90): E modified by roughness height where ks+ = u_τ * ks / ν (roughness Reynolds number)
Parameters
[in]friction_velocityFriction velocity u_τ
[in]roughness_heightEquivalent sand grain roughness height ks
[in]kinematic_viscosityKinematic viscosity ν
Returns
Roughness-modified log-law coefficient E
Note
This follows the Nikuradse sand grain roughness correlation

Computes roughness-modified log-law coefficient E.

Local to this translation unit.

Definition at line 140 of file wallfunction.c.

141{
142 // Compute roughness Reynolds number
143 double roughness_reynolds = friction_velocity * roughness_height / kinematic_viscosity;
144
145 double roughness_correction;
146
147 if (roughness_reynolds <= ROUGHNESS_TRANSITION_YPLUS) {
148 // Hydraulically smooth regime: no correction needed
149 roughness_correction = 0.0;
150 }
151 else if (roughness_reynolds < FULLY_ROUGH_YPLUS) {
152 // Transitional regime: smooth interpolation using sine function
153 // This provides a gradual transition between smooth and rough wall behavior
154 double max_correction = LOGLAW_B - 8.5 + (1.0 / KAPPA) * log(roughness_reynolds);
155
156 roughness_correction = max_correction *
157 sin(0.4258 * (log(roughness_reynolds) - 0.811));
158 }
159 else {
160 // Fully rough regime: maximum correction
161 roughness_correction = LOGLAW_B - 8.5 + (1.0 / KAPPA) * log(roughness_reynolds);
162 }
163
164 // Return E = exp[κ(B - ΔB)]
165 return exp(KAPPA * (LOGLAW_B - roughness_correction));
166}
#define ROUGHNESS_TRANSITION_YPLUS
Smooth-to-rough transition y+ threshold.
#define LOGLAW_B
Log-law intercept constant B for smooth walls.
#define FULLY_ROUGH_YPLUS
Fully rough regime y+ threshold.
#define KAPPA
von Karman constant (universal turbulence constant)
Here is the caller graph for this function:

◆ u_hydset_roughness()

double u_hydset_roughness ( double  kinematic_viscosity,
double  wall_distance,
double  friction_velocity,
double  roughness_height 
)

Computes velocity from log-law for rough walls.

Calculates the tangential velocity at a given wall distance using the roughness-modified log-law of the wall.

APPLICABLE REGIMES:

  • Viscous sublayer (y+ < 11.53): u+ = y+
  • Log-layer (11.53 < y+ < 300): u+ = (1/κ) ln(E * y+)
  • Outer layer (y+ > 300): Returns -1 (invalid, wake region)
Parameters
[in]kinematic_viscosityKinematic viscosity ν
[in]wall_distanceDistance from wall y
[in]friction_velocityFriction velocity u_τ
[in]roughness_heightEquivalent sand grain roughness ks
Returns
Tangential velocity u_t, or -1 if outside valid range
Note
The upper limit y+ < 300 ensures we stay within the log-layer
For y+ > 300, wake effects become important and log-law is invalid

Computes velocity from log-law for rough walls.

Local to this translation unit.

Definition at line 172 of file wallfunction.c.

174{
175 // Normalized wall distance
176 double yplus = friction_velocity * wall_distance / kinematic_viscosity;
177
178 // Threshold values for regime transitions
179 const double viscous_limit = VISCOUS_SUBLAYER_YPLUS;
180 const double loglayer_limit = 300.0;
181
182 double tangential_velocity;
183
184 if (yplus <= viscous_limit) {
185 // Viscous sublayer: linear velocity profile
186 tangential_velocity = friction_velocity * yplus;
187 }
188 else if (yplus <= loglayer_limit) {
189 // Log-layer: logarithmic velocity profile with roughness correction
190 double E = E_coeff(friction_velocity, roughness_height, kinematic_viscosity);
191 tangential_velocity = (friction_velocity / KAPPA) * log(E * yplus);
192 }
193 else {
194 // Outside valid range (wake region or beyond)
195 tangential_velocity = -1.0;
196 }
197
198 return tangential_velocity;
199}
#define VISCOUS_SUBLAYER_YPLUS
Viscous sublayer thickness y+ threshold.
double E_coeff(double friction_velocity, double roughness_height, double kinematic_viscosity)
Internal helper implementation: E_coeff().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ f_hydset()

double f_hydset ( double  kinematic_viscosity,
double  known_velocity,
double  wall_distance,
double  friction_velocity_guess,
double  roughness_height 
)

Residual function for friction velocity equation (log-law with roughness)

This function computes the residual for Newton-Raphson iteration: f(u_τ) = u_predicted(u_τ) - u_known where u_predicted comes from the log-law or linear law depending on y+.

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]known_velocityKnown velocity at distance y
[in]wall_distanceDistance from wall
[in]friction_velocity_guessCurrent guess for u_τ
[in]roughness_heightWall roughness height
Returns
Residual value (should be zero at solution)

Residual function for friction velocity equation (log-law with roughness)

Local to this translation unit.

Definition at line 209 of file wallfunction.c.

212{
213 double yplus = friction_velocity_guess * wall_distance / kinematic_viscosity;
214 double residual;
215
216 if (yplus <= VISCOUS_SUBLAYER_YPLUS) {
217 // Viscous sublayer: u = u_τ * y+
218 residual = friction_velocity_guess * yplus - known_velocity;
219 }
220 else {
221 // Log-layer: u = (u_τ / κ) * ln(E * y+)
222 double E = E_coeff(friction_velocity_guess, roughness_height, kinematic_viscosity);
223 residual = friction_velocity_guess * (1.0 / KAPPA * log(E * yplus)) - known_velocity;
224 }
225
226 return residual;
227}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ df_hydset()

double df_hydset ( double  kinematic_viscosity,
double  known_velocity,
double  wall_distance,
double  friction_velocity_guess,
double  roughness_height 
)

Numerical derivative of residual function.

Computes df/du_τ using finite differences for Newton-Raphson iteration.

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]known_velocityKnown velocity
[in]wall_distanceDistance from wall
[in]friction_velocity_guessCurrent guess for u_τ
[in]roughness_heightWall roughness
Returns
Derivative df/du_τ

Numerical derivative of residual function.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
df_hydset()

Definition at line 235 of file wallfunction.c.

238{
239 const double perturbation = 1.e-7;
240
241 double f_plus = f_hydset(kinematic_viscosity, known_velocity, wall_distance,
242 friction_velocity_guess + perturbation, roughness_height);
243
244 double f_minus = f_hydset(kinematic_viscosity, known_velocity, wall_distance,
245 friction_velocity_guess - perturbation, roughness_height);
246
247 return (f_plus - f_minus) / (2.0 * perturbation);
248}
double f_hydset(double kinematic_viscosity, double known_velocity, double wall_distance, double friction_velocity_guess, double roughness_height)
Internal helper implementation: f_hydset().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_utau_hydset()

double find_utau_hydset ( double  kinematic_viscosity,
double  known_velocity,
double  wall_distance,
double  initial_guess,
double  roughness_height 
)

Solves for friction velocity using Newton-Raphson iteration.

Given a known velocity at a known distance from the wall, this function iteratively solves for the friction velocity u_τ that satisfies the roughness-modified log-law or linear law.

ALGORITHM: Newton-Raphson: u_τ^(n+1) = u_τ^n - f(u_τ^n) / f'(u_τ^n) Convergence criterion: |u_τ^(n+1) - u_τ^n| < 1e-7

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]known_velocityVelocity at reference point
[in]wall_distanceDistance from wall to reference point
[in]initial_guessInitial guess for u_τ
[in]roughness_heightWall roughness height
Returns
Converged friction velocity u_τ
Note
Maximum 30 iterations; warns if convergence not achieved
Typical convergence in 3-5 iterations for good initial guess

Solves for friction velocity using Newton-Raphson iteration.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
find_utau_hydset()

Definition at line 256 of file wallfunction.c.

259{
260 double friction_velocity = initial_guess;
261 double friction_velocity_old = initial_guess;
262
263 const int max_iterations = 30;
264 const double convergence_tolerance = 1.e-7;
265
266 int iteration;
267 for (iteration = 0; iteration < max_iterations; iteration++) {
268 // Newton-Raphson update
269 double residual = f_hydset(kinematic_viscosity, known_velocity, wall_distance,
270 friction_velocity_old, roughness_height);
271 double derivative = df_hydset(kinematic_viscosity, known_velocity, wall_distance,
272 friction_velocity_old, roughness_height);
273
274 friction_velocity = friction_velocity_old - residual / derivative;
275
276 // Check convergence
277 if (fabs(friction_velocity - friction_velocity_old) < convergence_tolerance) {
278 break;
279 }
280
281 friction_velocity_old = friction_velocity;
282 }
283
284 // Warn if convergence not achieved
285 if (iteration == max_iterations) {
287 "WARNING: u_tau iteration did not converge after %d iterations. "
288 "Final difference: %le\n",
289 iteration, fabs(friction_velocity - friction_velocity_old));
290 }
291
292 return friction_velocity;
293}
#define LOCAL
Logging scope definitions for controlling message output.
Definition logging.h:45
#define LOG_ALLOW(scope, level, fmt,...)
Logging macro that checks both the log level and whether the calling function is in the allowed-funct...
Definition logging.h:200
@ LOG_DEBUG
Detailed debugging information.
Definition logging.h:32
double df_hydset(double kinematic_viscosity, double known_velocity, double wall_distance, double friction_velocity_guess, double roughness_height)
Implementation of df_hydset().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ nu_t()

double nu_t ( double  yplus)

Computes turbulent eddy viscosity ratio (ν_t / ν)

Uses the mixing length model with van Driest damping: ν_t / ν = κ * y+ * [1 - exp(-y+ / A+)]² where A+ ≈ 19 is the damping coefficient.

PHYSICAL INTERPRETATION:

  • Near wall (y+ → 0): ν_t → 0 (viscous effects dominate)
  • Far from wall (y+ >> A+): ν_t ~ κ * y+ (mixing length theory)
  • Damping function smoothly transitions between these limits
Parameters
[in]yplusNormalized wall distance y+ = y * u_τ / ν
Returns
Eddy viscosity ratio ν_t / ν
Note
This is the standard mixing length model with van Driest damping
Valid in the inner layer (y+ < 50-100)

Computes turbulent eddy viscosity ratio (ν_t / ν)

Local to this translation unit.

Definition at line 303 of file wallfunction.c.

304{
305 // van Driest damping function: [1 - exp(-y+ / A+)]²
306 double damping_function = 1.0 - exp(-yplus / DAMPING_COEFFICIENT);
307
308 // Mixing length model: ν_t / ν = κ * y+ * D²
309 return KAPPA * yplus * damping_function * damping_function;
310}
#define DAMPING_COEFFICIENT
Eddy viscosity damping coefficient (van Driest damping)
Here is the caller graph for this function:

◆ integrate_1()

double integrate_1 ( double  kinematic_viscosity,
double  wall_distance,
double  friction_velocity,
int  integration_mode 
)

Integrates eddy viscosity profile from wall to distance y.

Computes integrals needed for non-equilibrium wall functions: If mode=0: ∫[0 to y] dy / (ν + ν_t) If mode=1: ∫[0 to y] y dy / (ν + ν_t)

These integrals appear in the solution of the momentum equation with pressure gradients in the near-wall region.

NUMERICAL METHOD:

  • Trapezoidal rule with 30 integration points
  • Higher-order correction at endpoints for accuracy
Parameters
[in]kinematic_viscosityKinematic viscosity ν
[in]wall_distanceDistance from wall y
[in]friction_velocityFriction velocity u_τ
[in]integration_mode0 for F integral, 1 for Fy integral
Returns
Value of the integral
Note
Used in Cabot wall function for pressure gradient effects

Integrates eddy viscosity profile from wall to distance y.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
integrate_1()

Definition at line 318 of file wallfunction.c.

320{
321 const int num_points = 30;
322 double step_size = wall_distance / (num_points - 1);
323 double integral_values[num_points];
324
325 // Compute integrand at each point
326 for (int i = 0; i < num_points; i++) {
327 double current_distance = i * step_size;
328 double yplus = current_distance * friction_velocity / kinematic_viscosity;
329 double eddy_viscosity_ratio = nu_t(yplus);
330
331 if (integration_mode == 0) {
332 // Mode 0: integrand = 1 / [(1 + ν_t/ν) * ν]
333 integral_values[i] = 1.0 / ((1.0 + eddy_viscosity_ratio) * kinematic_viscosity);
334 }
335 else {
336 // Mode 1: integrand = y / [(1 + ν_t/ν) * ν]
337 integral_values[i] = current_distance /
338 ((1.0 + eddy_viscosity_ratio) * kinematic_viscosity);
339 }
340 }
341
342 // Trapezoidal rule integration with corrected endpoints
343 double integral_sum = 0.0;
344
345 // Interior points: weight = (f[i-1] + 2*f[i] + f[i+1]) * dy/4
346 for (int i = 1; i < num_points - 1; i++) {
347 integral_sum += (integral_values[i-1] + 2.0 * integral_values[i] +
348 integral_values[i+1]) * 0.25 * step_size;
349 }
350
351 // Endpoint correction for higher accuracy
352 integral_sum += 0.1667 * step_size *
353 (2.0 * integral_values[0] + integral_values[1] +
354 2.0 * integral_values[num_points-1] + integral_values[num_points-2]);
355
356 return integral_sum;
357}
double nu_t(double yplus)
Internal helper implementation: nu_t().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ taw()

double taw ( double  kinematic_viscosity,
double  friction_velocity,
double  wall_distance,
double  velocity,
double  pressure_gradient_tangent 
)

Computes wall shear stress with pressure gradient effects.

Solves the integrated momentum equation in the near-wall region: τ_w = (u - dp/dx * F_y) / F_1 where F_1 and F_y are integrals of the effective viscosity profile.

This accounts for non-equilibrium effects due to streamwise pressure gradients.

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]friction_velocityFriction velocity u_τ
[in]wall_distanceDistance from wall
[in]velocityVelocity at wall_distance
[in]pressure_gradient_tangentTangential pressure gradient dp/ds
Returns
Wall shear stress τ_w
Note
This is more accurate than equilibrium wall functions in separated flows

Computes wall shear stress with pressure gradient effects.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
taw()

Definition at line 369 of file wallfunction.c.

371{
372 // Compute integrated viscosity profiles
373 double F1 = integrate_1(kinematic_viscosity, wall_distance, friction_velocity, 0);
374 double Fy = integrate_1(kinematic_viscosity, wall_distance, friction_velocity, 1);
375
376 // Apply momentum balance: τ_w * F1 = u - (dp/dx) * Fy
377 return (1.0 / F1) * (velocity - pressure_gradient_tangent * Fy);
378}
double integrate_1(double kinematic_viscosity, double wall_distance, double friction_velocity, int integration_mode)
Implementation of integrate_1().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ u_Cabot()

double u_Cabot ( double  kinematic_viscosity,
double  wall_distance,
double  friction_velocity,
double  pressure_gradient_tangent,
double  wall_shear_stress 
)

Computes velocity using Cabot wall function.

Reconstructs velocity from wall shear stress and pressure gradient: u = τ_w * F1 + (dp/dx) * Fy

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]wall_distanceDistance from wall
[in]friction_velocityFriction velocity
[in]pressure_gradient_tangentTangential pressure gradient
[in]wall_shear_stressWall shear stress
Returns
Velocity at wall_distance

Computes velocity using Cabot wall function.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
u_Cabot()

Definition at line 386 of file wallfunction.c.

389{
390 double F1 = integrate_1(kinematic_viscosity, wall_distance, friction_velocity, 0);
391 double Fy = integrate_1(kinematic_viscosity, wall_distance, friction_velocity, 1);
392
393 return wall_shear_stress * F1 + pressure_gradient_tangent * Fy;
394}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ f_Cabot()

double f_Cabot ( double  kinematic_viscosity,
double  velocity,
double  wall_distance,
double  friction_velocity_guess,
double  pressure_gradient_tangent,
double  pressure_gradient_normal 
)

Residual function for Cabot wall function.

Computes: f(u_τ) = u_τ - sqrt(|τ_w(u_τ)|) This enforces consistency between u_τ and τ_w.

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]velocityVelocity
[in]wall_distanceDistance from wall
[in]friction_velocity_guessGuess for u_τ
[in]pressure_gradient_tangentTangential pressure gradient
[in]pressure_gradient_normalNormal pressure gradient (currently unused)
Returns
Residual value

Residual function for Cabot wall function.

Local to this translation unit.

Definition at line 400 of file wallfunction.c.

403{
404 (void)pressure_gradient_normal;
405 double wall_shear = taw(kinematic_viscosity, friction_velocity_guess,
406 wall_distance, velocity, pressure_gradient_tangent);
407
408 return friction_velocity_guess - sqrt(fabs(wall_shear));
409}
double taw(double kinematic_viscosity, double friction_velocity, double wall_distance, double velocity, double pressure_gradient_tangent)
Implementation of taw().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ df_Cabot()

double df_Cabot ( double  kinematic_viscosity,
double  velocity,
double  wall_distance,
double  friction_velocity_guess,
double  pressure_gradient_tangent,
double  pressure_gradient_normal 
)

Numerical derivative for Cabot wall function.

Parameters
kinematic_viscosityFluid kinematic viscosity.
velocityResolved tangential velocity at the wall-model point.
wall_distanceNormal distance from the wall.
friction_velocity_guessCurrent friction-velocity iterate.
pressure_gradient_tangentTangential pressure-gradient contribution.
pressure_gradient_normalNormal pressure-gradient contribution.
Returns
Derivative of the Cabot residual with respect to friction velocity.

Numerical derivative for Cabot wall function.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
df_Cabot()

Definition at line 417 of file wallfunction.c.

420{
421 const double perturbation = 1.e-7;
422
423 double f_plus = f_Cabot(kinematic_viscosity, velocity, wall_distance,
424 friction_velocity_guess + perturbation,
425 pressure_gradient_tangent, pressure_gradient_normal);
426
427 double f_minus = f_Cabot(kinematic_viscosity, velocity, wall_distance,
428 friction_velocity_guess - perturbation,
429 pressure_gradient_tangent, pressure_gradient_normal);
430
431 return (f_plus - f_minus) / (2.0 * perturbation);
432}
double f_Cabot(double kinematic_viscosity, double velocity, double wall_distance, double friction_velocity_guess, double pressure_gradient_tangent, double pressure_gradient_normal)
Internal helper implementation: f_Cabot().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_utau_Cabot()

void find_utau_Cabot ( double  kinematic_viscosity,
double  velocity,
double  wall_distance,
double  initial_guess,
double  pressure_gradient_tangent,
double  pressure_gradient_normal,
double *  friction_velocity,
double *  wall_shear_velocity,
double *  wall_shear_normal 
)

Solves for friction velocity using Cabot wall function.

This non-equilibrium wall function accounts for pressure gradient effects, making it more accurate in separated or strongly accelerating/decelerating flows.

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]velocityVelocity at reference point
[in]wall_distanceDistance from wall
[in]initial_guessInitial guess for u_τ
[in]pressure_gradient_tangentTangential pressure gradient
[in]pressure_gradient_normalNormal pressure gradient
[out]friction_velocityConverged friction velocity
[out]wall_shear_velocityWall shear stress for velocity
[out]wall_shear_normalWall shear stress for normal pressure gradient

Solves for friction velocity using Cabot wall function.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
find_utau_Cabot()

Definition at line 440 of file wallfunction.c.

444{
445 double current_guess = initial_guess;
446 double new_guess;
447
448 const int max_iterations = 30;
449 const double convergence_tolerance = 1.e-10;
450
451 int iteration;
452 for (iteration = 0; iteration < max_iterations; iteration++) {
453 double residual = f_Cabot(kinematic_viscosity, velocity, wall_distance,
454 current_guess, pressure_gradient_tangent,
455 pressure_gradient_normal);
456
457 double derivative = df_Cabot(kinematic_viscosity, velocity, wall_distance,
458 current_guess, pressure_gradient_tangent,
459 pressure_gradient_normal);
460
461 new_guess = fabs(current_guess - residual / derivative);
462
463 if (fabs(current_guess - new_guess) < convergence_tolerance) {
464 break;
465 }
466
467 current_guess = new_guess;
468 }
469
470 if (fabs(current_guess - new_guess) > 1.e-5 && iteration >= 29) {
472 "Cabot wall-function iteration reached its limit without the requested tolerance.\n");
473 }
474
475 *friction_velocity = new_guess;
476 *wall_shear_velocity = taw(kinematic_viscosity, new_guess, wall_distance,
477 velocity, pressure_gradient_tangent);
478 *wall_shear_normal = taw(kinematic_viscosity, new_guess, wall_distance,
479 0.0, pressure_gradient_normal);
480}
#define GLOBAL
Scope for global logging across all processes.
Definition logging.h:46
@ LOG_WARNING
Non-critical issues that warrant attention.
Definition logging.h:30
double df_Cabot(double kinematic_viscosity, double velocity, double wall_distance, double friction_velocity_guess, double pressure_gradient_tangent, double pressure_gradient_normal)
Implementation of df_Cabot().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ u_Werner()

double u_Werner ( double  kinematic_viscosity,
double  wall_distance,
double  friction_velocity 
)

Computes velocity using Werner-Wengle wall function.

Algebraic wall function that provides explicit relation: u+ = y+ for y+ < 11.81 (viscous sublayer) u+ = A * (y+)^B for y+ ≥ 11.81 (power law) where A = 8.3, B = 1/7 are empirical constants.

ADVANTAGES:

  • Explicit (no iteration required)
  • Computationally cheap
  • Reasonable accuracy for equilibrium boundary layers

LIMITATIONS:

  • Less accurate than log-law for high Reynolds numbers
  • No pressure gradient effects
Parameters
[in]kinematic_viscosityKinematic viscosity
[in]wall_distanceDistance from wall
[in]friction_velocityFriction velocity
Returns
Tangential velocity

Computes velocity using Werner-Wengle wall function.

Local to this translation unit.

Definition at line 490 of file wallfunction.c.

492{
493 double yplus = friction_velocity * wall_distance / kinematic_viscosity;
494
495 // Werner-Wengle constants
496 const double power_law_coefficient = 8.3;
497 const double power_law_exponent = 1.0 / 7.0;
498
499 double velocity;
500
501 if (yplus <= VISCOUS_SUBLAYER_YPLUS) {
502 // Viscous sublayer: u+ = y+
503 velocity = yplus * friction_velocity;
504 }
505 else {
506 // Power-law region: u+ = A * (y+)^B
507 velocity = power_law_coefficient * pow(yplus, power_law_exponent) * friction_velocity;
508 }
509
510 return velocity;
511}
Here is the caller graph for this function:

◆ f_Werner()

double f_Werner ( double  kinematic_viscosity,
double  velocity,
double  wall_distance,
double  friction_velocity 
)

Residual function for Werner-Wengle iteration.

Computes residual: f(u_τ) = u_τ² - g(u, y, ν) where g is derived from the velocity profile inversion.

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]velocityKnown velocity
[in]wall_distanceDistance from wall
[in]friction_velocityGuess for friction velocity
Returns
Residual value

Residual function for Werner-Wengle iteration.

Local to this translation unit.

Definition at line 517 of file wallfunction.c.

519{
520 const double power_law_coefficient = 8.3;
521 const double power_law_exponent = 1.0 / 7.0;
522
523 // Transition point between viscous and power-law regions
524 double transition_distance = VISCOUS_SUBLAYER_YPLUS * kinematic_viscosity / friction_velocity;
525
526 // Transition velocity
527 double transition_velocity = kinematic_viscosity / (2.0 * transition_distance) *
528 pow(power_law_coefficient, 2.0 / (1.0 - power_law_exponent));
529
530 double residual;
531
532 if (fabs(velocity) <= transition_velocity) {
533 // Viscous sublayer regime
534 residual = friction_velocity * friction_velocity -
535 velocity / wall_distance * kinematic_viscosity;
536 }
537 else {
538 // Power-law regime (more complex inversion formula)
539 double term1 = 0.5 * (1.0 - power_law_exponent) *
540 pow(power_law_coefficient, (1.0 + power_law_exponent) /
541 (1.0 - power_law_exponent)) *
542 pow(kinematic_viscosity / wall_distance, 1.0 + power_law_exponent);
543
544 double term2 = (1.0 + power_law_exponent) / power_law_coefficient *
545 pow(kinematic_viscosity / wall_distance, power_law_exponent) *
546 fabs(velocity);
547
548 residual = friction_velocity * friction_velocity -
549 velocity / fabs(velocity) * pow(term1 + term2, 2.0 / (1.0 + power_law_exponent));
550 }
551
552 return residual;
553}
Here is the caller graph for this function:

◆ df_Werner()

double df_Werner ( double  kinematic_viscosity,
double  velocity,
double  wall_distance,
double  friction_velocity 
)

Numerical derivative for Werner-Wengle iteration.

Parameters
kinematic_viscosityFluid kinematic viscosity.
velocityResolved velocity at the wall-model point.
wall_distanceNormal distance from the wall.
friction_velocityCurrent friction-velocity iterate.
Returns
Derivative of the Werner–Wengle residual with respect to friction velocity.

Numerical derivative for Werner-Wengle iteration.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
df_Werner()

Definition at line 561 of file wallfunction.c.

563{
564 const double perturbation = 1.e-7;
565
566 double f_plus = f_Werner(kinematic_viscosity, velocity, wall_distance,
567 friction_velocity + perturbation);
568
569 double f_minus = f_Werner(kinematic_viscosity, velocity, wall_distance,
570 friction_velocity - perturbation);
571
572 return (f_plus - f_minus) / (2.0 * perturbation);
573}
double f_Werner(double kinematic_viscosity, double velocity, double wall_distance, double friction_velocity)
Internal helper implementation: f_Werner().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_utau_Werner()

double find_utau_Werner ( double  kinematic_viscosity,
double  velocity,
double  wall_distance,
double  initial_guess 
)

Solves for friction velocity using Werner-Wengle wall function.

Parameters
[in]kinematic_viscosityKinematic viscosity
[in]velocityVelocity at reference point
[in]wall_distanceDistance from wall
[in]initial_guessInitial guess for u_τ
Returns
Converged friction velocity

Solves for friction velocity using Werner-Wengle wall function.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
find_utau_Werner()

Definition at line 581 of file wallfunction.c.

583{
584 double current_guess = initial_guess;
585 double new_guess;
586
587 const int max_iterations = 20;
588 const double convergence_tolerance = 1.e-7;
589
590 int iteration;
591 for (iteration = 0; iteration < max_iterations; iteration++) {
592 double residual = f_Werner(kinematic_viscosity, velocity, wall_distance, current_guess);
593 double derivative = df_Werner(kinematic_viscosity, velocity, wall_distance, current_guess);
594
595 new_guess = current_guess - residual / derivative;
596
597 if (fabs(current_guess - new_guess) < convergence_tolerance) {
598 break;
599 }
600
601 current_guess = new_guess;
602 }
603
604 if (fabs(current_guess - new_guess) > 1.e-5 && iteration >= 19) {
606 "Werner-Wengle wall-function iteration reached its limit without the requested tolerance.\n");
607 }
608
609 return new_guess;
610}
double df_Werner(double kinematic_viscosity, double velocity, double wall_distance, double friction_velocity)
Implementation of df_Werner().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ u_loglaw()

double u_loglaw ( double  wall_distance,
double  friction_velocity,
double  roughness_length 
)

Computes velocity using simple log-law (smooth wall with roughness offset)

Simple logarithmic profile: u = (u_τ / κ) * ln((y + y0) / y0) where y0 is a roughness length scale.

Parameters
[in]wall_distanceDistance from wall
[in]friction_velocityFriction velocity
[in]roughness_lengthRoughness length scale y0
Returns
Velocity
Note
This is a simplified model; use u_hydset_roughness for better accuracy

Computes velocity using simple log-law (smooth wall with roughness offset)

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
u_loglaw()

Definition at line 622 of file wallfunction.c.

623{
624 return friction_velocity * (1.0 / KAPPA) * log((roughness_length + wall_distance) / roughness_length);
625}
Here is the caller graph for this function:

◆ find_utau_loglaw()

double find_utau_loglaw ( double  velocity,
double  wall_distance,
double  roughness_length 
)

Solves for friction velocity using simple log-law (explicit formula)

Explicit inversion: u_τ = κ * u / ln((y + y0) / y0)

Parameters
[in]velocityKnown velocity
[in]wall_distanceDistance from wall
[in]roughness_lengthRoughness length scale
Returns
Friction velocity

Solves for friction velocity using simple log-law (explicit formula)

Local to this translation unit.

Definition at line 631 of file wallfunction.c.

632{
633 return KAPPA * velocity / log((wall_distance + roughness_length) / roughness_length);
634}
Here is the caller graph for this function:

◆ wall_function()

void wall_function ( UserCtx user,
double  distance_reference,
double  distance_boundary,
Cmpnts  velocity_wall,
Cmpnts  velocity_reference,
Cmpnts velocity_boundary,
PetscReal *  friction_velocity,
double  normal_x,
double  normal_y,
double  normal_z 
)

Applies standard wall function with Werner-Wengle model.

This is a high-level interface that:

  1. Decomposes velocity into normal/tangential components
  2. Applies wall function to tangential velocity
  3. Interpolates normal velocity
  4. Reconstructs total velocity
Parameters
[in]userSimulation context
[in]distance_referenceDistance to reference point
[in]distance_boundaryDistance to boundary point
[in]velocity_wallWall velocity
[in]velocity_referenceReference velocity
[out]velocity_boundaryOutput boundary velocity
[out]friction_velocityOutput friction velocity
[in]normal_xX-component of wall normal
[in]normal_yY-component of wall normal
[in]normal_zZ-component of wall normal

Applies standard wall function with Werner-Wengle model.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/wallfunction.h.

See also
wall_function()

Definition at line 650 of file wallfunction.c.

654{
655 SimCtx *simCtx = user->simCtx;
656 double kinematic_viscosity = 1.0 / simCtx->ren;
657
658 // Decompose velocity into components
659 double delta_u = velocity_reference.x - velocity_wall.x;
660 double delta_v = velocity_reference.y - velocity_wall.y;
661 double delta_w = velocity_reference.z - velocity_wall.z;
662
663 double normal_velocity = delta_u * normal_x + delta_v * normal_y + delta_w * normal_z;
664
665 double tangential_u = delta_u - normal_velocity * normal_x;
666 double tangential_v = delta_v - normal_velocity * normal_y;
667 double tangential_w = delta_w - normal_velocity * normal_z;
668
669 double tangential_magnitude = sqrt(tangential_u * tangential_u +
670 tangential_v * tangential_v +
671 tangential_w * tangential_w);
672
673 // Apply Werner-Wengle wall function
674 double utau_guess = 0.05;
675 double tangential_modeled = u_Werner(kinematic_viscosity, distance_boundary, utau_guess);
676 if (friction_velocity) {
677 // Werner-Wengle path currently uses a fixed u_tau estimate.
678 *friction_velocity = (PetscReal)utau_guess;
679 }
680
681 // Scale tangential components
682 if (tangential_magnitude > 1.e-10) {
683 tangential_u *= tangential_modeled / tangential_magnitude;
684 tangential_v *= tangential_modeled / tangential_magnitude;
685 tangential_w *= tangential_modeled / tangential_magnitude;
686 }
687 else {
688 tangential_u = tangential_v = tangential_w = 0.0;
689 }
690
691 // Reconstruct velocity: U = U_t + U_n * n
692 (*velocity_boundary).x = tangential_u + (distance_boundary / distance_reference) * normal_velocity * normal_x;
693 (*velocity_boundary).y = tangential_v + (distance_boundary / distance_reference) * normal_velocity * normal_y;
694 (*velocity_boundary).z = tangential_w + (distance_boundary / distance_reference) * normal_velocity * normal_z;
695
696 // Add wall velocity
697 (*velocity_boundary).x += velocity_wall.x;
698 (*velocity_boundary).y += velocity_wall.y;
699 (*velocity_boundary).z += velocity_wall.z;
700}
SimCtx * simCtx
Back-pointer to the master simulation context.
Definition variables.h:909
PetscReal ren
Definition variables.h:744
The master context for the entire simulation.
Definition variables.h:695
double u_Werner(double kinematic_viscosity, double wall_distance, double friction_velocity)
Internal helper implementation: u_Werner().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ wall_function_loglaw()

void wall_function_loglaw ( UserCtx user,
double  roughness_height,
double  distance_reference,
double  distance_boundary,
Cmpnts  velocity_wall,
Cmpnts  velocity_reference,
Cmpnts velocity_boundary,
PetscReal *  friction_velocity,
double  normal_x,
double  normal_y,
double  normal_z 
)

Applies log-law wall function with roughness correction.

This is the recommended wall function interface for most applications. It uses the roughness-corrected log-law which is accurate for both smooth and rough walls.

Parameters
[in]userSimulation context
[in]roughness_heightWall roughness height ks
[in]distance_referenceDistance to reference point
[in]distance_boundaryDistance to boundary point
[in]velocity_wallWall velocity (typically zero)
[in]velocity_referenceReference velocity
[out]velocity_boundaryOutput boundary velocity
[out]friction_velocityOutput friction velocity u_τ
[in]normal_xX-component of wall normal
[in]normal_yY-component of wall normal
[in]normal_zZ-component of wall normal
Note
Sets velocity to reference value if y+ > 300 (outside log-layer)

Applies log-law wall function with roughness correction.

Local to this translation unit.

Definition at line 706 of file wallfunction.c.

711{
712 SimCtx *simCtx = user->simCtx;
713 double kinematic_viscosity = 1.0 / simCtx->ren;
714
715 // Decompose velocity
716 double delta_u = velocity_reference.x - velocity_wall.x;
717 double delta_v = velocity_reference.y - velocity_wall.y;
718 double delta_w = velocity_reference.z - velocity_wall.z;
719
720 double normal_velocity = delta_u * normal_x + delta_v * normal_y + delta_w * normal_z;
721
722 double tangential_u = delta_u - normal_velocity * normal_x;
723 double tangential_v = delta_v - normal_velocity * normal_y;
724 double tangential_w = delta_w - normal_velocity * normal_z;
725
726 double tangential_magnitude = sqrt(tangential_u * tangential_u +
727 tangential_v * tangential_v +
728 tangential_w * tangential_w);
729
730 // Solve for friction velocity
731 *friction_velocity = find_utau_hydset(kinematic_viscosity, tangential_magnitude,
732 distance_reference, 0.001, roughness_height);
733
734 // Compute wall function velocity
735 double tangential_modeled = u_hydset_roughness(kinematic_viscosity, distance_boundary,
736 *friction_velocity, roughness_height);
737
738 // Check if outside valid range (y+ > 300)
739 if (tangential_modeled < 0.0) {
740 tangential_modeled = tangential_magnitude;
741 }
742
743 // Scale tangential components
744 if (tangential_magnitude > 1.e-10) {
745 tangential_u *= tangential_modeled / tangential_magnitude;
746 tangential_v *= tangential_modeled / tangential_magnitude;
747 tangential_w *= tangential_modeled / tangential_magnitude;
748 }
749 else {
750 tangential_u = tangential_v = tangential_w = 0.0;
751 }
752
753 // Reconstruct total velocity
754 (*velocity_boundary).x = tangential_u + (distance_boundary / distance_reference) * normal_velocity * normal_x;
755 (*velocity_boundary).y = tangential_v + (distance_boundary / distance_reference) * normal_velocity * normal_y;
756 (*velocity_boundary).z = tangential_w + (distance_boundary / distance_reference) * normal_velocity * normal_z;
757
758 // Add wall velocity
759 (*velocity_boundary).x += velocity_wall.x;
760 (*velocity_boundary).y += velocity_wall.y;
761 (*velocity_boundary).z += velocity_wall.z;
762}
double u_hydset_roughness(double kinematic_viscosity, double wall_distance, double friction_velocity, double roughness_height)
Internal helper implementation: u_hydset_roughness().
double find_utau_hydset(double kinematic_viscosity, double known_velocity, double wall_distance, double initial_guess, double roughness_height)
Implementation of find_utau_hydset().
Here is the call graph for this function:
Here is the caller graph for this function:

◆ wall_function_Cabot()

void wall_function_Cabot ( UserCtx user,
double  roughness_height,
double  distance_reference,
double  distance_boundary,
Cmpnts  velocity_wall,
Cmpnts  velocity_reference,
Cmpnts velocity_boundary,
PetscReal *  friction_velocity,
double  normal_x,
double  normal_y,
double  normal_z,
double  pressure_gradient_x,
double  pressure_gradient_y,
double  pressure_gradient_z,
int  iteration_count 
)

Applies Cabot non-equilibrium wall function with pressure gradients.

This wall function accounts for pressure gradient effects and is most accurate in separated or strongly accelerating/decelerating flows.

Parameters
[in]userSimulation context
[in]roughness_heightWall roughness (currently unused in this version)
[in]distance_referenceDistance to reference point
[in]distance_boundaryDistance to boundary point
[in]velocity_wallWall velocity
[in]velocity_referenceReference velocity
[out]velocity_boundaryOutput boundary velocity
[out]friction_velocityOutput friction velocity
[in]normal_xX-component of wall normal
[in]normal_yY-component of wall normal
[in]normal_zZ-component of wall normal
[in]pressure_gradient_xX-component of pressure gradient
[in]pressure_gradient_yY-component of pressure gradient
[in]pressure_gradient_zZ-component of pressure gradient
[in]iteration_countCurrent iteration count
Note
Only updates friction velocity every 5 iterations for stability
Normal pressure gradient currently set to zero

Applies Cabot non-equilibrium wall function with pressure gradients.

Local to this translation unit.

Definition at line 768 of file wallfunction.c.

775{
776 (void)roughness_height;
777 SimCtx *simCtx = user->simCtx;
778 double kinematic_viscosity = 1.0 / simCtx->ren;
779
780 // Decompose velocity
781 double delta_u = velocity_reference.x - velocity_wall.x;
782 double delta_v = velocity_reference.y - velocity_wall.y;
783 double delta_w = velocity_reference.z - velocity_wall.z;
784
785 double normal_velocity = delta_u * normal_x + delta_v * normal_y + delta_w * normal_z;
786
787 double tangential_u = delta_u - normal_velocity * normal_x;
788 double tangential_v = delta_v - normal_velocity * normal_y;
789 double tangential_w = delta_w - normal_velocity * normal_z;
790
791 double tangential_magnitude = sqrt(tangential_u * tangential_u +
792 tangential_v * tangential_v +
793 tangential_w * tangential_w);
794
795 // Compute tangential pressure gradient
796 double pressure_gradient_tangent = 0.0;
797 if (tangential_magnitude > 1.e-10) {
798 pressure_gradient_tangent = (pressure_gradient_x * tangential_u +
799 pressure_gradient_y * tangential_v +
800 pressure_gradient_z * tangential_w) / tangential_magnitude;
801 }
802
803 // Normal pressure gradient (currently set to zero)
804 double pressure_gradient_normal = 0.0;
805
806 // Solve for friction velocity (only on first iteration or periodically)
807 if (iteration_count == 0 || iteration_count > 4) {
808 double utau, wall_shear1, wall_shear2;
809 find_utau_Cabot(kinematic_viscosity, tangential_magnitude, distance_reference,
810 0.01, pressure_gradient_tangent, pressure_gradient_normal,
811 &utau, &wall_shear1, &wall_shear2);
812 *friction_velocity = utau;
813 }
814
815 // Compute wall function velocity with pressure gradient correction
816 double tangential_modeled = u_Cabot(kinematic_viscosity, distance_boundary,
817 *friction_velocity, pressure_gradient_tangent,
818 taw(kinematic_viscosity, *friction_velocity,
819 distance_reference, tangential_magnitude,
820 pressure_gradient_tangent));
821
822 // Scale tangential components
823 if (tangential_magnitude > 1.e-10) {
824 tangential_u *= tangential_modeled / tangential_magnitude;
825 tangential_v *= tangential_modeled / tangential_magnitude;
826 tangential_w *= tangential_modeled / tangential_magnitude;
827 }
828 else {
829 tangential_u = tangential_v = tangential_w = 0.0;
830 }
831
832 // Reconstruct total velocity
833 (*velocity_boundary).x = tangential_u + (distance_boundary / distance_reference) * normal_velocity * normal_x;
834 (*velocity_boundary).y = tangential_v + (distance_boundary / distance_reference) * normal_velocity * normal_y;
835 (*velocity_boundary).z = tangential_w + (distance_boundary / distance_reference) * normal_velocity * normal_z;
836
837 // Add wall velocity
838 (*velocity_boundary).x += velocity_wall.x;
839 (*velocity_boundary).y += velocity_wall.y;
840 (*velocity_boundary).z += velocity_wall.z;
841}
double u_Cabot(double kinematic_viscosity, double wall_distance, double friction_velocity, double pressure_gradient_tangent, double wall_shear_stress)
Implementation of u_Cabot().
void find_utau_Cabot(double kinematic_viscosity, double velocity, double wall_distance, double initial_guess, double pressure_gradient_tangent, double pressure_gradient_normal, double *friction_velocity, double *wall_shear_velocity, double *wall_shear_normal)
Implementation of find_utau_Cabot().
Here is the call graph for this function:
Here is the caller graph for this function: