wallfunction.c File Reference

Wall function implementations for near-wall turbulence modeling. More...

#include "wallfunction.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

Include dependency graph for wallfunction.c:

Go to the source code of this file.

Macros
#define	KAPPA   0.41
	von Karman constant (universal turbulence constant)
#define	LOGLAW_B   5.5
	Log-law intercept constant B for smooth walls.
#define	VISCOUS_SUBLAYER_YPLUS   11.81
	Viscous sublayer thickness y+ threshold.
#define	ROUGHNESS_TRANSITION_YPLUS   2.25
	Smooth-to-rough transition y+ threshold.
#define	FULLY_ROUGH_YPLUS   90.0
	Fully rough regime y+ threshold.
#define	DAMPING_COEFFICIENT   19.0
	Eddy viscosity damping coefficient (van Driest damping)
#define	__FUNCT__   "taw_Werner"
#define	__FUNCT__   "u_Werner_explicit"

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)
	Internal helper implementation: noslip().
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)
	Internal helper implementation: freeslip().
double	E_coeff (double friction_velocity, double roughness_height, double kinematic_viscosity)
	Internal helper implementation: E_coeff().
double	u_hydset_roughness (double kinematic_viscosity, double wall_distance, double friction_velocity, double roughness_height)
	Internal helper implementation: u_hydset_roughness().
double	f_hydset (double kinematic_viscosity, double known_velocity, double wall_distance, double friction_velocity_guess, double roughness_height)
	Internal helper implementation: f_hydset().
double	df_hydset (double kinematic_viscosity, double known_velocity, double wall_distance, double friction_velocity_guess, double roughness_height)
	Implementation of df_hydset().
double	find_utau_hydset (double kinematic_viscosity, double known_velocity, double wall_distance, double initial_guess, double roughness_height)
	Implementation of find_utau_hydset().
double	nu_t (double yplus)
	Internal helper implementation: nu_t().
double	integrate_1 (double kinematic_viscosity, double wall_distance, double friction_velocity, int integration_mode)
	Implementation of integrate_1().
double	taw (double kinematic_viscosity, double friction_velocity, double wall_distance, double velocity, double pressure_gradient_tangent)
	Implementation of taw().
double	u_Cabot (double kinematic_viscosity, double wall_distance, double friction_velocity, double pressure_gradient_tangent, double wall_shear_stress)
	Implementation of u_Cabot().
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().
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().
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().
double	u_Werner (double kinematic_viscosity, double wall_distance, double friction_velocity)
	Internal helper implementation: u_Werner().
double	taw_Werner (double kinematic_viscosity, double velocity, double wall_distance)
	Implementation of taw_Werner().
double	u_Werner_explicit (double kinematic_viscosity, double wall_distance, double wall_shear_stress)
	Implementation of u_Werner_explicit().
double	find_utau_Werner (double kinematic_viscosity, double velocity, double wall_distance, double initial_guess)
	Implementation of find_utau_Werner().
double	u_loglaw (double wall_distance, double friction_velocity, double roughness_length)
	Implementation of u_loglaw().
double	find_utau_loglaw (double velocity, double wall_distance, double roughness_length)
	Internal helper implementation: find_utau_loglaw().
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)
	Implementation of wall_function().
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)
	Internal helper implementation: wall_function_loglaw().
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)
	Internal helper implementation: wall_function_Cabot().

Detailed Description

Wall function implementations for near-wall turbulence modeling.

This file contains various wall function models that bridge the gap between the wall and the first computational grid point. Wall functions allow simulations to avoid resolving the viscous sublayer, reducing computational cost while maintaining reasonable accuracy for turbulent wall-bounded flows.

PHYSICAL BACKGROUND: The near-wall region in turbulent flows is characterized by three layers:

1. Viscous sublayer (y+ < 5): Dominated by viscous effects, u+ = y+
2. Buffer layer (5 < y+ < 30): Transition region
3. Log-law region (y+ > 30): Inertial layer, u+ = (1/κ) ln(y+) + B where: u+ = u / u_τ (normalized velocity) y+ = y * u_τ / ν (normalized wall distance) u_τ = sqrt(τ_w / ρ) (friction velocity) κ ≈ 0.41 (von Karman constant) B ≈ 5.5 (log-law intercept for smooth walls)

IMPLEMENTED MODELS:

• No-slip: Linear interpolation for stationary walls
• Free-slip: Normal interpolation, tangential extrapolation
• Werner-Wengle: Algebraic wall function
• Log-law: Standard equilibrium wall function with roughness
• Cabot: Non-equilibrium wall function with pressure gradient effects

Author

Original implementation adapted from legacy code

Date

Enhanced with comprehensive documentation

Definition in file wallfunction.c.

Macro Definition Documentation

KAPPA

#define KAPPA   0.41

von Karman constant (universal turbulence constant)

Definition at line 43 of file wallfunction.c.

LOGLAW_B

#define LOGLAW_B   5.5

Log-law intercept constant B for smooth walls.

Definition at line 46 of file wallfunction.c.

VISCOUS_SUBLAYER_YPLUS

#define VISCOUS_SUBLAYER_YPLUS   11.81

Viscous sublayer thickness y+ threshold.

Definition at line 49 of file wallfunction.c.

ROUGHNESS_TRANSITION_YPLUS

#define ROUGHNESS_TRANSITION_YPLUS   2.25

Smooth-to-rough transition y+ threshold.

Definition at line 52 of file wallfunction.c.

FULLY_ROUGH_YPLUS

#define FULLY_ROUGH_YPLUS   90.0

Fully rough regime y+ threshold.

Definition at line 55 of file wallfunction.c.

DAMPING_COEFFICIENT

#define DAMPING_COEFFICIENT   19.0

Eddy viscosity damping coefficient (van Driest damping)

Definition at line 58 of file wallfunction.c.

__FUNCT__ [1/2]

#define __FUNCT__   "taw_Werner"

Definition at line 517 of file wallfunction.c.

__FUNCT__ [2/2]

#define __FUNCT__   "u_Werner_explicit"

Definition at line 517 of file wallfunction.c.

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
	)

Internal helper implementation: noslip().

Applies no-slip wall boundary condition with linear interpolation.

Local to this translation unit.

Definition at line 68 of file wallfunction.c.

{
    (void)user;
    (void)normal_x;
    (void)normal_y;
    (void)normal_z;
    // Compute velocity difference between reference point and wall
    double delta_u = velocity_reference.x - velocity_wall.x;
    double delta_v = velocity_reference.y - velocity_wall.y;
    double delta_w = velocity_reference.z - velocity_wall.z;
    
    // Linear interpolation factor
    double interpolation_factor = distance_boundary / distance_reference;
    
    // Apply linear interpolation
    (*velocity_boundary).x = interpolation_factor * delta_u;
    (*velocity_boundary).y = interpolation_factor * delta_v;
    (*velocity_boundary).z = interpolation_factor * delta_w;
    
    // Add wall velocity (shift to absolute frame)
    (*velocity_boundary).x += velocity_wall.x;
    (*velocity_boundary).y += velocity_wall.y;
    (*velocity_boundary).z += velocity_wall.z;
}

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
	)

Internal helper implementation: freeslip().

Applies free-slip wall boundary condition.

Local to this translation unit.

Definition at line 101 of file wallfunction.c.

{
    (void)user;
    // Extract normal components of velocities
    double wall_normal_velocity = velocity_wall.x * normal_x + 
                                   velocity_wall.y * normal_y + 
                                   velocity_wall.z * normal_z;
    
    double reference_normal_velocity = velocity_reference.x * normal_x + 
                                        velocity_reference.y * normal_y + 
                                        velocity_reference.z * normal_z;
    
    // Interpolate normal velocity component
    double boundary_normal_velocity = wall_normal_velocity + 
        (reference_normal_velocity - wall_normal_velocity) * (distance_boundary / distance_reference);
    
    // Extract tangential velocity components (extrapolated from reference point)
    double tangential_u = velocity_reference.x - reference_normal_velocity * normal_x;
    double tangential_v = velocity_reference.y - reference_normal_velocity * normal_y;
    double tangential_w = velocity_reference.z - reference_normal_velocity * normal_z;
    
    // Reconstruct total velocity: U = U_t + U_n * n
    (*velocity_boundary).x = tangential_u + boundary_normal_velocity * normal_x;
    (*velocity_boundary).y = tangential_v + boundary_normal_velocity * normal_y;
    (*velocity_boundary).z = tangential_w + boundary_normal_velocity * normal_z;
}

Here is the caller graph for this function:

E_coeff()

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

Internal helper implementation: E_coeff().

Computes roughness-modified log-law coefficient E.

Local to this translation unit.

Definition at line 140 of file wallfunction.c.

{
    // Compute roughness Reynolds number
    double roughness_reynolds = friction_velocity * roughness_height / kinematic_viscosity;
    
    double roughness_correction;
    
    if (roughness_reynolds <= ROUGHNESS_TRANSITION_YPLUS) {
        // Hydraulically smooth regime: no correction needed
        roughness_correction = 0.0;
    }
    else if (roughness_reynolds < FULLY_ROUGH_YPLUS) {
        // Transitional regime: smooth interpolation using sine function
        // This provides a gradual transition between smooth and rough wall behavior
        double max_correction = LOGLAW_B - 8.5 + (1.0 / KAPPA) * log(roughness_reynolds);
        
        roughness_correction = max_correction * 
            sin(0.4258 * (log(roughness_reynolds) - 0.811));
    }
    else {
        // Fully rough regime: maximum correction
        roughness_correction = LOGLAW_B - 8.5 + (1.0 / KAPPA) * log(roughness_reynolds);
    }
    
    // Return E = exp[κ(B - ΔB)]
    return exp(KAPPA * (LOGLAW_B - roughness_correction));
}

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
	)

Internal helper implementation: u_hydset_roughness().

Computes velocity from log-law for rough walls.

Local to this translation unit.

Definition at line 172 of file wallfunction.c.

{
    // Normalized wall distance
    double yplus = friction_velocity * wall_distance / kinematic_viscosity;
    
    // Threshold values for regime transitions
    const double viscous_limit = VISCOUS_SUBLAYER_YPLUS;
    const double loglayer_limit = 300.0;
    
    double tangential_velocity;
    
    if (yplus <= viscous_limit) {
        // Viscous sublayer: linear velocity profile
        tangential_velocity = friction_velocity * yplus;
    }
    else if (yplus <= loglayer_limit) {
        // Log-layer: logarithmic velocity profile with roughness correction
        double E = E_coeff(friction_velocity, roughness_height, kinematic_viscosity);
        tangential_velocity = (friction_velocity / KAPPA) * log(E * yplus);
    }
    else {
        // Outside valid range (wake region or beyond)
        tangential_velocity = -1.0;
    }
    
    return tangential_velocity;
}

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
	)

Internal helper implementation: f_hydset().

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

Local to this translation unit.

Definition at line 209 of file wallfunction.c.

{
    double yplus = friction_velocity_guess * wall_distance / kinematic_viscosity;
    double residual;
    
    if (yplus <= VISCOUS_SUBLAYER_YPLUS) {
        // Viscous sublayer: u = u_τ * y+
        residual = friction_velocity_guess * yplus - known_velocity;
    }
    else {
        // Log-layer: u = (u_τ / κ) * ln(E * y+)
        double E = E_coeff(friction_velocity_guess, roughness_height, kinematic_viscosity);
        residual = friction_velocity_guess * (1.0 / KAPPA * log(E * yplus)) - known_velocity;
    }
    
    return residual;
}

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
	)

Implementation of df_hydset().

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.

{
    const double perturbation = 1.e-7;
    
    double f_plus = f_hydset(kinematic_viscosity, known_velocity, wall_distance,
                            friction_velocity_guess + perturbation, roughness_height);
    
    double f_minus = f_hydset(kinematic_viscosity, known_velocity, wall_distance,
                             friction_velocity_guess - perturbation, roughness_height);
    
    return (f_plus - f_minus) / (2.0 * perturbation);
}

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
	)

Implementation of find_utau_hydset().

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.

{
    double friction_velocity = initial_guess;
    double friction_velocity_old = initial_guess;
    
    const int max_iterations = 30;
    const double convergence_tolerance = 1.e-7;
    
    int iteration;
    for (iteration = 0; iteration < max_iterations; iteration++) {
        // Newton-Raphson update
        double residual = f_hydset(kinematic_viscosity, known_velocity, wall_distance,
                                  friction_velocity_old, roughness_height);
        double derivative = df_hydset(kinematic_viscosity, known_velocity, wall_distance,
                                     friction_velocity_old, roughness_height);
        
        friction_velocity = friction_velocity_old - residual / derivative;
        
        // Check convergence
        if (fabs(friction_velocity - friction_velocity_old) < convergence_tolerance) {
            break;
        }
        
        friction_velocity_old = friction_velocity;
    }
    
    // Warn if convergence not achieved
    if (iteration == max_iterations) {
        LOG_ALLOW(LOCAL, LOG_DEBUG, 
                 "WARNING: u_tau iteration did not converge after %d iterations. "
                 "Final difference: %le\n", 
                 iteration, fabs(friction_velocity - friction_velocity_old));
    }
    
    return friction_velocity;
}

Here is the call graph for this function:

Here is the caller graph for this function:

nu_t()

double nu_t

(

double

yplus

)

Internal helper implementation: nu_t().

Computes turbulent eddy viscosity ratio (ν_t / ν)

Local to this translation unit.

Definition at line 303 of file wallfunction.c.

{
    // van Driest damping function: [1 - exp(-y+ / A+)]²
    double damping_function = 1.0 - exp(-yplus / DAMPING_COEFFICIENT);
    
    // Mixing length model: ν_t / ν = κ * y+ * D²
    return KAPPA * yplus * damping_function * damping_function;
}

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
	)

Implementation of integrate_1().

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.

{
    const int num_points = 30;
    double step_size = wall_distance / (num_points - 1);
    double integral_values[num_points];
    
    // Compute integrand at each point
    for (int i = 0; i < num_points; i++) {
        double current_distance = i * step_size;
        double yplus = current_distance * friction_velocity / kinematic_viscosity;
        double eddy_viscosity_ratio = nu_t(yplus);
        
        if (integration_mode == 0) {
            // Mode 0: integrand = 1 / [(1 + ν_t/ν) * ν]
            integral_values[i] = 1.0 / ((1.0 + eddy_viscosity_ratio) * kinematic_viscosity);
        }
        else {
            // Mode 1: integrand = y / [(1 + ν_t/ν) * ν]
            integral_values[i] = current_distance / 
                                ((1.0 + eddy_viscosity_ratio) * kinematic_viscosity);
        }
    }
    
    // Trapezoidal rule integration with corrected endpoints
    double integral_sum = 0.0;
    
    // Interior points: weight = (f[i-1] + 2*f[i] + f[i+1]) * dy/4
    for (int i = 1; i < num_points - 1; i++) {
        integral_sum += (integral_values[i-1] + 2.0 * integral_values[i] + 
                        integral_values[i+1]) * 0.25 * step_size;
    }
    
    // Endpoint correction for higher accuracy
    integral_sum += 0.1667 * step_size * 
                   (2.0 * integral_values[0] + integral_values[1] +
                    2.0 * integral_values[num_points-1] + integral_values[num_points-2]);
    
    return integral_sum;
}

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
	)

Implementation of taw().

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.

{
    // Compute integrated viscosity profiles
    double F1 = integrate_1(kinematic_viscosity, wall_distance, friction_velocity, 0);
    double Fy = integrate_1(kinematic_viscosity, wall_distance, friction_velocity, 1);
    
    // Apply momentum balance: τ_w * F1 = u - (dp/dx) * Fy
    return (1.0 / F1) * (velocity - pressure_gradient_tangent * Fy);
}

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
	)

Implementation of u_Cabot().

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.

{
    double F1 = integrate_1(kinematic_viscosity, wall_distance, friction_velocity, 0);
    double Fy = integrate_1(kinematic_viscosity, wall_distance, friction_velocity, 1);
    
    return wall_shear_stress * F1 + pressure_gradient_tangent * Fy;
}

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
	)

Internal helper implementation: f_Cabot().

Residual function for Cabot wall function.

Local to this translation unit.

Definition at line 400 of file wallfunction.c.

{
    (void)pressure_gradient_normal;
    double wall_shear = taw(kinematic_viscosity, friction_velocity_guess,
                           wall_distance, velocity, pressure_gradient_tangent);
    
    return friction_velocity_guess - sqrt(fabs(wall_shear));
}

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
	)

Implementation of df_Cabot().

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.

{
    const double perturbation = 1.e-7;
    
    double f_plus = f_Cabot(kinematic_viscosity, velocity, wall_distance,
                           friction_velocity_guess + perturbation,
                           pressure_gradient_tangent, pressure_gradient_normal);
    
    double f_minus = f_Cabot(kinematic_viscosity, velocity, wall_distance,
                            friction_velocity_guess - perturbation,
                            pressure_gradient_tangent, pressure_gradient_normal);
    
    return (f_plus - f_minus) / (2.0 * perturbation);
}

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
	)

Implementation of find_utau_Cabot().

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.

{
    double current_guess = initial_guess;
    double new_guess;
    
    const int max_iterations = 30;
    const double convergence_tolerance = 1.e-10;
    
    int iteration;
    for (iteration = 0; iteration < max_iterations; iteration++) {
        double residual = f_Cabot(kinematic_viscosity, velocity, wall_distance,
                                 current_guess, pressure_gradient_tangent,
                                 pressure_gradient_normal);
        
        double derivative = df_Cabot(kinematic_viscosity, velocity, wall_distance,
                                    current_guess, pressure_gradient_tangent,
                                    pressure_gradient_normal);
        
        new_guess = fabs(current_guess - residual / derivative);
        
        if (fabs(current_guess - new_guess) < convergence_tolerance) {
            break;
        }
        
        current_guess = new_guess;
    }
    
    if (fabs(current_guess - new_guess) > 1.e-5 && iteration >= 29) {
        LOG_ALLOW(GLOBAL, LOG_WARNING,
                  "Cabot wall-function iteration reached its limit without the requested tolerance.\n");
    }
    
    *friction_velocity = new_guess;
    *wall_shear_velocity = taw(kinematic_viscosity, new_guess, wall_distance,
                              velocity, pressure_gradient_tangent);
    *wall_shear_normal = taw(kinematic_viscosity, new_guess, wall_distance,
                            0.0, pressure_gradient_normal);
}

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
	)

Internal helper implementation: u_Werner().

Computes velocity using Werner-Wengle wall function.

Local to this translation unit.

Definition at line 490 of file wallfunction.c.

{
    /* The pointwise profile, and the production path: `wall_function()` evaluates it at
       the boundary cell with the friction velocity `find_utau_Werner()` takes from the
       reference cell, so one relation is used in both directions. */
    double yplus = friction_velocity * wall_distance / kinematic_viscosity;
    
    // Werner-Wengle constants
    const double power_law_coefficient = 8.3;
    const double power_law_exponent = 1.0 / 7.0;
    
    double velocity;
    
    if (yplus <= VISCOUS_SUBLAYER_YPLUS) {
        // Viscous sublayer: u+ = y+
        velocity = yplus * friction_velocity;
    }
    else {
        // Power-law region: u+ = A * (y+)^B
        velocity = power_law_coefficient * pow(yplus, power_law_exponent) * friction_velocity;
    }
    
    return velocity;
}

Here is the caller graph for this function:

taw_Werner()

double taw_Werner	(	double	kinematic_viscosity,
		double	velocity,
		double	wall_distance
	)

Implementation of taw_Werner().

Wall shear stress from the cell-integrated Werner-Wengle relation.

Full API contract is documented with the header declaration in include/wallfunction.h.

Definition at line 523 of file wallfunction.c.

{
    const double A = 8.3;
    const double B = 1.0 / 7.0;
    const double speed = fabs(velocity);
    const double nu_over_y = kinematic_viscosity / wall_distance;
 
    /* The switch is written in terms of the geometry alone. Expressing it through the
       friction velocity, as the iterative residual does, makes the branch depend on the
       answer and forces a Newton solve the closed form does not need. */
    const double transition_speed = 0.5 * nu_over_y * pow(A, 2.0 / (1.0 - B));
 
    if (speed <= transition_speed) {
        /* Integrating u+ = y+ across the cell leaves the factor of two that separates
           this from the pointwise relation tau_w = nu u / y. */
        return 2.0 * kinematic_viscosity * speed / wall_distance;
    }
 
    {
        const double term1 = 0.5 * (1.0 - B) * pow(A, (1.0 + B) / (1.0 - B)) *
                             pow(nu_over_y, 1.0 + B);
        const double term2 = (1.0 + B) / A * pow(nu_over_y, B) * speed;
 
        return pow(term1 + term2, 2.0 / (1.0 + B));
    }
}

Here is the caller graph for this function:

u_Werner_explicit()

double u_Werner_explicit	(	double	kinematic_viscosity,
		double	wall_distance,
		double	wall_shear_stress
	)

Implementation of u_Werner_explicit().

Wall-parallel velocity at a distance, from the Werner-Wengle wall stress.

Full API contract is documented with the header declaration in include/wallfunction.h.

Definition at line 557 of file wallfunction.c.

{
    const double A = 8.3;
    const double B = 1.0 / 7.0;
    const double stress = fabs(wall_shear_stress);
    const double nu_over_y = kinematic_viscosity / wall_distance;
 
    /* The same threshold as taw_Werner(), carried over to the stress side. The two
       branches meet here exactly, so selecting on the stress cannot land on the wrong
       side of the velocity switch. */
    const double transition_stress = nu_over_y * nu_over_y * pow(A, 2.0 / (1.0 - B));
 
    if (stress <= transition_stress) {
        return stress * wall_distance / (2.0 * kinematic_viscosity);
    }
 
    {
        const double offset = 0.5 * (1.0 - B) * pow(A, (1.0 + B) / (1.0 - B)) *
                              pow(nu_over_y, 1.0 + B);
        /* Non-negative by construction above the threshold - it is zero only in the
           limit of the switch - so the max guards roundoff, not a modelling case. */
        const double bracket = PetscMax(pow(stress, 0.5 * (1.0 + B)) - offset, 0.0);
 
        return A / (1.0 + B) * pow(1.0 / nu_over_y, B) * bracket;
    }
}

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
	)

Implementation of find_utau_Werner().

Friction velocity from a point velocity, inverting u_Werner() exactly.

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 591 of file wallfunction.c.

{
    const double A = 8.3;
    const double B = 1.0 / 7.0;
    const double speed = fabs(velocity);
    const double y_over_nu = wall_distance / kinematic_viscosity;
    (void)initial_guess;
 
    /* The exact inverse of `u_Werner()`, branch for branch. Along that profile
       U y / nu = (y+)^2 in the sublayer, so the branch it was evaluated on is decided by
       the data alone: the sublayer holds while U y / nu <= VISCOUS_SUBLAYER_YPLUS^2, the
       same switch `u_Werner()` applies to y+. Both branches then invert explicitly. */
    if (speed * y_over_nu <= VISCOUS_SUBLAYER_YPLUS * VISCOUS_SUBLAYER_YPLUS) {
        return sqrt(speed / y_over_nu);                           /* U = u_tau^2 y / nu */
    }
    return pow(speed / (A * pow(y_over_nu, B)), 1.0 / (1.0 + B)); /* U = A u_tau^(1+B) (y/nu)^B */
}

Here is the caller graph for this function:

u_loglaw()

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

Implementation of u_loglaw().

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 620 of file wallfunction.c.

{
    return friction_velocity * (1.0 / KAPPA) * log((roughness_length + wall_distance) / roughness_length);
}

Here is the caller graph for this function:

find_utau_loglaw()

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

Internal helper implementation: find_utau_loglaw().

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

Local to this translation unit.

Definition at line 629 of file wallfunction.c.

{
    return KAPPA * velocity / log((wall_distance + roughness_length) / roughness_length);
}

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
	)

Implementation of wall_function().

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 648 of file wallfunction.c.

{
    SimCtx *simCtx = user->simCtx;
    double kinematic_viscosity = 1.0 / simCtx->ren;
    
    // Decompose velocity into components
    double delta_u = velocity_reference.x - velocity_wall.x;
    double delta_v = velocity_reference.y - velocity_wall.y;
    double delta_w = velocity_reference.z - velocity_wall.z;
    
    double normal_velocity = delta_u * normal_x + delta_v * normal_y + delta_w * normal_z;
    
    double tangential_u = delta_u - normal_velocity * normal_x;
    double tangential_v = delta_v - normal_velocity * normal_y;
    double tangential_w = delta_w - normal_velocity * normal_z;
    
    double tangential_magnitude = sqrt(tangential_u * tangential_u +
                                      tangential_v * tangential_v +
                                      tangential_w * tangential_w);
    
    /* The reference speed is a point value - the velocity at the reference cell's centre,
       `distance_reference` from the wall - so the wall stress comes from the pointwise
       Werner-Wengle profile inverted there, and the boundary cell's speed from that same
       profile at its own centre. One relation in both directions, as the log-law path
       does, and both directions are closed form, so there is still no inner iteration.
       The cell-integral relation (`taw_Werner()`) reads its input as the average over
       [0, distance_reference]; given this point value it over-predicted u_tau by about
       (1 + 1/7)^(7/8), 12-14%, on a wall-modelled channel at Re_tau ~ 1000. */
    const double friction = find_utau_Werner(kinematic_viscosity, tangential_magnitude,
                                             distance_reference, 0.0);
    double tangential_modeled = u_Werner(kinematic_viscosity, distance_boundary, friction);
    if (friction_velocity) {
        *friction_velocity = (PetscReal)friction;
    }
    
    // Scale tangential components
    if (tangential_magnitude > 1.e-10) {
        tangential_u *= tangential_modeled / tangential_magnitude;
        tangential_v *= tangential_modeled / tangential_magnitude;
        tangential_w *= tangential_modeled / tangential_magnitude;
    }
    else {
        tangential_u = tangential_v = tangential_w = 0.0;
    }
    
    // Reconstruct velocity: U = U_t + U_n * n
    (*velocity_boundary).x = tangential_u + (distance_boundary / distance_reference) * normal_velocity * normal_x;
    (*velocity_boundary).y = tangential_v + (distance_boundary / distance_reference) * normal_velocity * normal_y;
    (*velocity_boundary).z = tangential_w + (distance_boundary / distance_reference) * normal_velocity * normal_z;
    
    // Add wall velocity
    (*velocity_boundary).x += velocity_wall.x;
    (*velocity_boundary).y += velocity_wall.y;
    (*velocity_boundary).z += velocity_wall.z;
}

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
	)

Internal helper implementation: wall_function_loglaw().

Applies log-law wall function with roughness correction.

Local to this translation unit.

Definition at line 711 of file wallfunction.c.

{
    SimCtx *simCtx = user->simCtx;
    double kinematic_viscosity = 1.0 / simCtx->ren;
    
    // Decompose velocity
    double delta_u = velocity_reference.x - velocity_wall.x;
    double delta_v = velocity_reference.y - velocity_wall.y;
    double delta_w = velocity_reference.z - velocity_wall.z;
    
    double normal_velocity = delta_u * normal_x + delta_v * normal_y + delta_w * normal_z;
    
    double tangential_u = delta_u - normal_velocity * normal_x;
    double tangential_v = delta_v - normal_velocity * normal_y;
    double tangential_w = delta_w - normal_velocity * normal_z;
    
    double tangential_magnitude = sqrt(tangential_u * tangential_u +
                                      tangential_v * tangential_v +
                                      tangential_w * tangential_w);
    
    // Solve for friction velocity
    *friction_velocity = find_utau_hydset(kinematic_viscosity, tangential_magnitude,
                                         distance_reference, 0.001, roughness_height);
    
    // Compute wall function velocity
    double tangential_modeled = u_hydset_roughness(kinematic_viscosity, distance_boundary,
                                                   *friction_velocity, roughness_height);
    
    // Check if outside valid range (y+ > 300)
    if (tangential_modeled < 0.0) {
        tangential_modeled = tangential_magnitude;
    }
    
    // Scale tangential components
    if (tangential_magnitude > 1.e-10) {
        tangential_u *= tangential_modeled / tangential_magnitude;
        tangential_v *= tangential_modeled / tangential_magnitude;
        tangential_w *= tangential_modeled / tangential_magnitude;
    }
    else {
        tangential_u = tangential_v = tangential_w = 0.0;
    }
    
    // Reconstruct total velocity
    (*velocity_boundary).x = tangential_u + (distance_boundary / distance_reference) * normal_velocity * normal_x;
    (*velocity_boundary).y = tangential_v + (distance_boundary / distance_reference) * normal_velocity * normal_y;
    (*velocity_boundary).z = tangential_w + (distance_boundary / distance_reference) * normal_velocity * normal_z;
    
    // Add wall velocity
    (*velocity_boundary).x += velocity_wall.x;
    (*velocity_boundary).y += velocity_wall.y;
    (*velocity_boundary).z += velocity_wall.z;
}

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
	)

Internal helper implementation: wall_function_Cabot().

Applies Cabot non-equilibrium wall function with pressure gradients.

Local to this translation unit.

Definition at line 773 of file wallfunction.c.

{
    (void)roughness_height;
    SimCtx *simCtx = user->simCtx;
    double kinematic_viscosity = 1.0 / simCtx->ren;
    
    // Decompose velocity
    double delta_u = velocity_reference.x - velocity_wall.x;
    double delta_v = velocity_reference.y - velocity_wall.y;
    double delta_w = velocity_reference.z - velocity_wall.z;
    
    double normal_velocity = delta_u * normal_x + delta_v * normal_y + delta_w * normal_z;
    
    double tangential_u = delta_u - normal_velocity * normal_x;
    double tangential_v = delta_v - normal_velocity * normal_y;
    double tangential_w = delta_w - normal_velocity * normal_z;
    
    double tangential_magnitude = sqrt(tangential_u * tangential_u +
                                      tangential_v * tangential_v +
                                      tangential_w * tangential_w);
    
    // Compute tangential pressure gradient
    double pressure_gradient_tangent = 0.0;
    if (tangential_magnitude > 1.e-10) {
        pressure_gradient_tangent = (pressure_gradient_x * tangential_u +
                                    pressure_gradient_y * tangential_v +
                                    pressure_gradient_z * tangential_w) / tangential_magnitude;
    }
    
    // Normal pressure gradient (currently set to zero)
    double pressure_gradient_normal = 0.0;
    
    // Solve for friction velocity (only on first iteration or periodically)
    if (iteration_count == 0 || iteration_count > 4) {
        double utau, wall_shear1, wall_shear2;
        find_utau_Cabot(kinematic_viscosity, tangential_magnitude, distance_reference,
                       0.01, pressure_gradient_tangent, pressure_gradient_normal,
                       &utau, &wall_shear1, &wall_shear2);
        *friction_velocity = utau;
    }
    
    // Compute wall function velocity with pressure gradient correction
    double tangential_modeled = u_Cabot(kinematic_viscosity, distance_boundary,
                                       *friction_velocity, pressure_gradient_tangent,
                                       taw(kinematic_viscosity, *friction_velocity,
                                           distance_reference, tangential_magnitude,
                                           pressure_gradient_tangent));
    
    // Scale tangential components
    if (tangential_magnitude > 1.e-10) {
        tangential_u *= tangential_modeled / tangential_magnitude;
        tangential_v *= tangential_modeled / tangential_magnitude;
        tangential_w *= tangential_modeled / tangential_magnitude;
    }
    else {
        tangential_u = tangential_v = tangential_w = 0.0;
    }
    
    // Reconstruct total velocity
    (*velocity_boundary).x = tangential_u + (distance_boundary / distance_reference) * normal_velocity * normal_x;
    (*velocity_boundary).y = tangential_v + (distance_boundary / distance_reference) * normal_velocity * normal_y;
    (*velocity_boundary).z = tangential_w + (distance_boundary / distance_reference) * normal_velocity * normal_z;
    
    // Add wall velocity
    (*velocity_boundary).x += velocity_wall.x;
    (*velocity_boundary).y += velocity_wall.y;
    (*velocity_boundary).z += velocity_wall.z;
}

Here is the call graph for this function:

Here is the caller graph for this function: