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	taw_Werner (double kinematic_viscosity, double velocity, double wall_distance)
	Wall shear stress from the cell-integrated Werner-Wengle relation.
double	u_Werner_explicit (double kinematic_viscosity, double wall_distance, double wall_shear_stress)
	Wall-parallel velocity at a distance, from the Werner-Wengle wall stress.
double	find_utau_Werner (double kinematic_viscosity, double velocity, double wall_distance, double initial_guess)
	Friction velocity from a point velocity, inverting u_Werner() exactly.
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]	user	Simulation context (unused but required for interface)
[in]	distance_reference	Wall-normal distance to reference point (sc)
[in]	distance_boundary	Wall-normal distance to boundary point (sb)
[in]	velocity_wall	Velocity at the wall (Ua), typically zero
[in]	velocity_reference	Velocity at reference point (Uc)
[out]	velocity_boundary	Computed velocity at boundary point (Ub)
[in]	normal_x	X-component of wall normal vector
[in]	normal_y	Y-component of wall normal vector
[in]	normal_z	Z-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.

{
    (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
	)

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]	user	Simulation context
[in]	distance_reference	Wall-normal distance to reference point
[in]	distance_boundary	Wall-normal distance to boundary point
[in]	velocity_wall	Velocity at the wall
[in]	velocity_reference	Velocity at reference point
[out]	velocity_boundary	Computed velocity at boundary point
[in]	normal_x	X-component of wall normal vector
[in]	normal_y	Y-component of wall normal vector
[in]	normal_z	Z-component of wall normal vector

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
	)

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_velocity	Friction velocity u_τ
[in]	roughness_height	Equivalent sand grain roughness height ks
[in]	kinematic_viscosity	Kinematic 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.

{
    // 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
	)

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_viscosity	Kinematic viscosity ν
[in]	wall_distance	Distance from wall y
[in]	friction_velocity	Friction velocity u_τ
[in]	roughness_height	Equivalent 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.

{
    // 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
	)

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_viscosity	Kinematic viscosity
[in]	known_velocity	Known velocity at distance y
[in]	wall_distance	Distance from wall
[in]	friction_velocity_guess	Current guess for u_τ
[in]	roughness_height	Wall 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.

{
    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
	)

Numerical derivative of residual function.

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

Parameters

[in]	kinematic_viscosity	Kinematic viscosity
[in]	known_velocity	Known velocity
[in]	wall_distance	Distance from wall
[in]	friction_velocity_guess	Current guess for u_τ
[in]	roughness_height	Wall 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.

{
    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
	)

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_viscosity	Kinematic viscosity
[in]	known_velocity	Velocity at reference point
[in]	wall_distance	Distance from wall to reference point
[in]	initial_guess	Initial guess for u_τ
[in]	roughness_height	Wall 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.

{
    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

)

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]

yplus

Normalized 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.

{
    // 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
	)

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_viscosity	Kinematic viscosity ν
[in]	wall_distance	Distance from wall y
[in]	friction_velocity	Friction velocity u_τ
[in]	integration_mode	0 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.

{
    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
	)

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_viscosity	Kinematic viscosity
[in]	friction_velocity	Friction velocity u_τ
[in]	wall_distance	Distance from wall
[in]	velocity	Velocity at wall_distance
[in]	pressure_gradient_tangent	Tangential 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.

{
    // 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
	)

Computes velocity using Cabot wall function.

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

Parameters

[in]	kinematic_viscosity	Kinematic viscosity
[in]	wall_distance	Distance from wall
[in]	friction_velocity	Friction velocity
[in]	pressure_gradient_tangent	Tangential pressure gradient
[in]	wall_shear_stress	Wall 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.

{
    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
	)

Residual function for Cabot wall function.

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

Parameters

[in]	kinematic_viscosity	Kinematic viscosity
[in]	velocity	Velocity
[in]	wall_distance	Distance from wall
[in]	friction_velocity_guess	Guess for u_τ
[in]	pressure_gradient_tangent	Tangential pressure gradient
[in]	pressure_gradient_normal	Normal 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.

{
    (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
	)

Numerical derivative for Cabot wall function.

Parameters

kinematic_viscosity	Fluid kinematic viscosity.
velocity	Resolved tangential velocity at the wall-model point.
wall_distance	Normal distance from the wall.
friction_velocity_guess	Current friction-velocity iterate.
pressure_gradient_tangent	Tangential pressure-gradient contribution.
pressure_gradient_normal	Normal 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.

{
    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
	)

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_viscosity	Kinematic viscosity
[in]	velocity	Velocity at reference point
[in]	wall_distance	Distance from wall
[in]	initial_guess	Initial guess for u_τ
[in]	pressure_gradient_tangent	Tangential pressure gradient
[in]	pressure_gradient_normal	Normal pressure gradient
[out]	friction_velocity	Converged friction velocity
[out]	wall_shear_velocity	Wall shear stress for velocity
[out]	wall_shear_normal	Wall 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.

{
    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
	)

Computes velocity using Werner-Wengle wall function.

The pointwise Werner-Wengle profile, and the relation the solver uses: the wall function evaluates it at the boundary cell with the friction velocity that find_utau_Werner() takes from the reference cell. 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_viscosity	Kinematic viscosity
[in]	wall_distance	Distance from wall
[in]	friction_velocity	Friction velocity

Returns

Tangential velocity

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
	)

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

Note

Not wired into the solver. The profile is assumed to hold across a cell of height wall_distance and integrated, so velocity must be the AVERAGE over [0, wall_distance] - a finite-volume first-cell value with its full height. The wall function's reference is instead a point value at the second cell's centre; read as an average it over-predicts u_tau by about (1 + B)^(1/(1+B)), 12-14% on a wall-modelled channel at Re_tau ~ 1000. The solver uses the pointwise pair find_utau_Werner() / u_Werner(), which is also closed form. Kept, with its inverse u_Werner_explicit(), for a caller that does supply a cell average.

|u| <= u_m : tau_w = 2 nu |u| / y |u| > u_m : tau_w = [ (1-B)/2 A^((1+B)/(1-B)) (nu/y)^(1+B)

• (1+B)/A (nu/y)^B |u| ]^(2/(1+B))

with A = 8.3, B = 1/7, and the switch at u_m = nu/(2y) A^(2/(1-B)). That threshold depends only on the geometry and the fluid, which is what makes the inversion explicit; a threshold written in terms of the friction velocity would reintroduce the iteration the model exists to avoid. The two branches meet exactly at u_m.

Pairs with u_Werner_explicit(), which inverts it. Do not pair it with u_Werner(), which is the pointwise profile rather than its cell integral.

Parameters

[in]	kinematic_viscosity	Kinematic viscosity.
[in]	velocity	Wall-parallel speed at wall_distance, sign ignored.
[in]	wall_distance	Normal distance the profile is integrated over.

Returns

Wall shear stress, equal to the square of the friction velocity.

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
	)

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

The exact inverse of taw_Werner(), branch for branch, so that a stress obtained from a velocity at one distance reproduces a consistent velocity at another. It is monotonically increasing in wall_distance, so a cell nearer the wall than the one the stress came from is always assigned the slower speed.

The branch is selected on the stress rather than the velocity it is about to produce, using the same threshold expressed as tau_m = (nu/y)^2 A^(2/(1-B)).

Parameters

[in]	kinematic_viscosity	Kinematic viscosity.
[in]	wall_distance	Normal distance to evaluate at.
[in]	wall_shear_stress	Wall shear stress from taw_Werner().

Returns

Wall-parallel speed at wall_distance, non-negative.

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
	)

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

The closed-form inverse of the pointwise Werner-Wengle profile, branch for branch, and the solver's path. Along that profile U y / nu = (y+)^2 in the sublayer, so the branch is chosen from the data alone - sublayer while U y / nu <= 11.81^2 - and each branch inverts explicitly:

sublayer : u_tau = sqrt(nu |U| / y) power law : u_tau = ( |U| / (A (y/nu)^B) )^(1/(1+B)), A = 8.3, B = 1/7

No iteration, so the residual it feeds stays smooth for a matrix-free Jacobian.

Parameters

[in]	kinematic_viscosity	Kinematic viscosity
[in]	velocity	Wall-parallel speed at the reference point, sign ignored
[in]	wall_distance	Distance from the wall to that point
[in]	initial_guess	Unused; kept so existing callers need no change

Returns

Friction velocity

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
	)

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_distance	Distance from wall
[in]	friction_velocity	Friction velocity
[in]	roughness_length	Roughness 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 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
	)

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

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

Parameters

[in]	velocity	Known velocity
[in]	wall_distance	Distance from wall
[in]	roughness_length	Roughness length scale

Returns

Friction velocity

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
	)

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]	user	Simulation context
[in]	distance_reference	Distance to reference point
[in]	distance_boundary	Distance to boundary point
[in]	velocity_wall	Wall velocity
[in]	velocity_reference	Reference velocity
[out]	velocity_boundary	Output boundary velocity
[out]	friction_velocity	Output friction velocity
[in]	normal_x	X-component of wall normal
[in]	normal_y	Y-component of wall normal
[in]	normal_z	Z-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 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
	)

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]	user	Simulation context
[in]	roughness_height	Wall roughness height ks
[in]	distance_reference	Distance to reference point
[in]	distance_boundary	Distance to boundary point
[in]	velocity_wall	Wall velocity (typically zero)
[in]	velocity_reference	Reference velocity
[out]	velocity_boundary	Output boundary velocity
[out]	friction_velocity	Output friction velocity u_τ
[in]	normal_x	X-component of wall normal
[in]	normal_y	Y-component of wall normal
[in]	normal_z	Z-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 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
	)

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]	user	Simulation context
[in]	roughness_height	Wall roughness (currently unused in this version)
[in]	distance_reference	Distance to reference point
[in]	distance_boundary	Distance to boundary point
[in]	velocity_wall	Wall velocity
[in]	velocity_reference	Reference velocity
[out]	velocity_boundary	Output boundary velocity
[out]	friction_velocity	Output friction velocity
[in]	normal_x	X-component of wall normal
[in]	normal_y	Y-component of wall normal
[in]	normal_z	Z-component of wall normal
[in]	pressure_gradient_x	X-component of pressure gradient
[in]	pressure_gradient_y	Y-component of pressure gradient
[in]	pressure_gradient_z	Z-component of pressure gradient
[in]	iteration_count	Current 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 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: