Skip to content

3-Element Windkessel Boundary Conditions: Interactive Tutorial

Reference

This tutorial is based on the comprehensive review by Westerhof, Lankhaar & Westerhof (2009) "The arterial Windkessel" in Medical & Biological Engineering & Computing, 47:131-141. DOI: 10.1007/s11517-008-0359-2

Overview

The 3-element Windkessel model (RCR) represents the optimal balance between accuracy and simplicity for cardiovascular CFD simulations. This model uses three lumped electrical circuit elements that represent physiological properties with clear physical meaning, providing accurate systolic pressure prediction while maintaining computational efficiency.

Historical Context

Historical Development

  • Hales (1735): First measured blood pressure and noted elasticity's role in pressure variations
  • Weber (1827): Compared arterial elasticity to fire engine windkessels
  • Frank (1899): Quantitatively formulated the original two-element model
  • Modern Development: The 3-element RCR model evolved to address systolic pressure accuracy limitations

Key Concepts

Physiological Basis

Large arteries (aorta, carotids) provide compliance - storing blood during systole and releasing during diastole

Resistance vessels (arterioles) create peripheral resistance following Poiseuille's law (R ∝ 1/r⁴)

3-Element RCR Model Components

  • R: Total peripheral resistance = (P̄ₐₒ - P̄ᵥₑₙ)/CO
  • C: Total arterial compliance = ΔV/ΔP
  • Zc: Characteristic impedance = ρ·PWV/A (proximal resistance)

Frequency Response

  • Low freq (< 0.2 HR): Peripheral resistance dominates
  • Mid freq (0.2-3 HR): Total arterial compliance controls
  • High freq (> 3 HR): Characteristic impedance determines response

Clinical Relevance

Pulse pressure is a major predictor of cardiovascular mortality. The 3-element RCR model provides optimal accuracy for afterload modeling by combining resistance, compliance, and wave propagation effects.

The 3-Element RCR Model

The 3-element Windkessel model represents the gold standard for cardiovascular CFD boundary conditions, providing excellent accuracy while maintaining computational efficiency.

Mathematical Foundation

Mathematical Model

Governing Equation:

\[P(t) = Z_c \cdot Q(t) + P_d(t)\]

Distal Pressure (Windkessel compartment):

\[P_d(t) = R \cdot Q_d(t) + \frac{1}{C}\int_0^t Q_d(\tau) d\tau\]

Flow Conservation:

\[Q(t) = Q_d(t) + C \cdot \frac{dP_d}{dt}\]

Characteristic Impedance:

\[Z_c = \frac{\rho \cdot c}{A} = \sqrt{\frac{\rho \cdot E \cdot h}{2 \cdot r \cdot A}}\]

Where: - ρ = blood density (1060 kg/m³) - c = pulse wave velocity (4-8 m/s) - A = cross-sectional area - E = elastic modulus - h = wall thickness - r = vessel radius

Components: - Zc: Characteristic impedance (ρ·PWV/A, ~5-7% of R) - C: Total arterial compliance (aorta + large arteries) - R: Total peripheral resistance (arterioles)

Key Advantages

Why Choose the 3-Element RCR Model?

  • Accurate systolic pressure prediction: Characteristic impedance captures wave propagation effects
  • Physiologically meaningful: Each parameter represents distinct vascular properties
  • Computationally efficient: Simple enough for routine CFD simulations
  • Well-validated: Extensively tested in cardiovascular research
  • Optimal balance: Best compromise between accuracy and simplicity

Numerical Implementation

The 3-element RCR model requires proper discretization for CFD applications.

Discretization for CFD

1st Order (Euler)

\[P^{n+1} = P^n + \Delta t \left[ \frac{Q^n - Q^{n-1}}{C \Delta t} + Z \frac{Q^n - Q^{n-1}}{\Delta t} \right] + R Q^n\]

2nd Order (Adams-Bashforth)

\[P^{n+1} = P^n + \frac{\Delta t}{2C} \left[ 3Q^n - Q^{n-1} \right] + Z \frac{3Q^n - Q^{n-1}}{2} + R Q^n\]

3rd Order (Higher Accuracy)

\[P^{n+1} = P^n + \frac{\Delta t}{12C} \left[ 23Q^n - 16Q^{n-1} + 5Q^{n-2} \right] + Z \frac{23Q^n - 16Q^{n-1} + 5Q^{n-2}}{12} + R Q^n\]

RCR Parameter Relationships

Parameter Physical Meaning Typical Range Units
R Peripheral resistance 10⁶ - 10⁹ Pa·s/m³
C Arterial compliance 10⁻¹⁰ - 10⁻⁸ m³/Pa
Zc Characteristic impedance 10⁵ - 10⁷ Pa·s/m³

Clinical Insight from Westerhof et al.

"At mean flow (0 Hz), peripheral resistance determines input impedance. For intermediate frequencies (0.2-3 HR), total arterial compliance dominates. At high frequencies (>3 HR), only characteristic impedance matters."

RCR Parameter Estimation

Multiple validated methods exist for estimating 3-element Windkessel parameters, as comprehensively reviewed by Westerhof et al.

Clinical Methods

Clinical Data Approach for RCR Model

  1. Decay Time Method: RC from diastolic pressure exponential decay
  2. Pulse Pressure Method: C ≈ SV/PP (overestimates by ~60%)
  3. Area Method: RC = ∫P dt / (P₁ - P₂) - reduces noise
  4. Complete RCR Fitting: Simultaneous estimation of R, C, and Zc from pressure/flow data
  5. Characteristic Impedance: Zc ≈ 5-7% of total peripheral resistance R

Literature-Based Values

Typical RCR Values for Aortic Outlets (from published studies):

Outlet Location R [×10⁸ Pa·s/m³] C [×10⁻⁹ m³/Pa] Zc [×10⁶ Pa·s/m³]
Descending Aorta 0.5 - 1.0 1.0 - 2.0 0.5 - 1.0
Brachiocephalic 1.0 - 2.0 0.5 - 1.0 1.0 - 2.0
Carotid Arteries 2.0 - 4.0 0.3 - 0.8 2.0 - 4.0
Subclavian 3.0 - 6.0 0.2 - 0.5 3.0 - 6.0 
// Adjust for patient age
scalar age_factor = 1.0 + 0.02 * (age - 50); // 2% per year after 50
R_adjusted = R_baseline * age_factor;
C_adjusted = C_baseline / age_factor;  // Decreased compliance

Murray's Law Scaling

// Based on outlet diameter relative to reference for RCR model
scalar diameter = 0.025;  // Outlet diameter [m]
scalar d_ref = 0.030;     // Reference diameter [m]
scalar scaling = pow(d_ref/diameter, 4);  // Poiseuille scaling

// Scale RCR parameters
scalar R = R_reference * scaling;
scalar C = C_reference / scaling;
scalar Zc = Zc_reference * pow(d_ref/diameter, 2);

// Flow split consideration
scalar flow_fraction = 0.65;  // 65% to descending aorta
R_effective = R / flow_fraction;
C_effective = C * flow_fraction;

OpenFOAM Implementation

The modularWKPressure boundary condition implements the 3-element RCR Windkessel model in OpenFOAM with backflow stabilization.

Installation Steps

# 1. Clone the OpenFOAM-WK repository
git clone https://github.com/JieWangnk/OpenFOAM-WK.git
cd OpenFOAM-WK

# 2. Source OpenFOAM environment
source /opt/openfoam12/etc/bashrc

# 3. Compile the boundary condition
cd src/modularWKPressure
wmake

Pressure Boundary Configuration (0/p)

outlet
{
    type            modularWKPressure;
    phi             phi;           // Flux field name
    order           2;             // Time discretization order (1-3)

    // 3-Element RCR Windkessel parameters
    R               8.0e7;         // Peripheral resistance [Pa·s/m³]
    C               1.5e-9;        // Compliance [m³/Pa]
    Z               1.2e7;         // Characteristic impedance [Pa·s/m³]

    // Initial conditions
    p0              12000;         // Initial pressure [Pa]
    value           uniform 12000;

    // Historical values for higher-order accuracy
    p_1             12000;         // Pressure at t-dt
    q_1             0;             // Flow at t-dt
    q_2             0;             // Flow at t-2dt
    q_3             0;             // Flow at t-3dt
}

Velocity Boundary Configuration (0/U)

outlet
{
    type                stabilizedWindkesselVelocity;
    beta                1.0;       // Stabilization coefficient
    enableStabilization true;      // Enable backflow stabilization
}

Stabilization Parameters

Stabilization Coefficients

  • β = 0.5: Conservative stabilization for mild backflow
  • β = 1.0: Standard stabilization (recommended)
  • β = 1.5: Strong stabilization for severe backflow issues

Backflow Stabilization

Complex 3D cardiovascular geometries often experience flow reversal, leading to numerical instabilities. The stabilized implementation prevents these issues.

The Stabilization Method

Based on Esmaily Moghadam et al. (2011), the stabilization adds a velocity penalty term:

\[\mathbf{u} = \mathbf{u}_{WK} - \beta \frac{\max(0, -\mathbf{u} \cdot \mathbf{n})}{\|\mathbf{u}\|} \mathbf{u}\]

Where: - \(\mathbf{u}_{WK}\): Standard Windkessel velocity - \(\beta\): Stabilization coefficient - \(\mathbf{n}\): Outward normal vector

Performance Impact

CoA Test Case Results:

Configuration Max Δt Stability Simulation Time
No Stabilization 1e-8 s ❌ Crashes N/A
β = 0.5 1e-5 s ⚠️ Marginal 10x slower
β = 1.0 1e-4 s ✅ Stable Normal
β = 1.5 1e-4 s ✅ Very Stable Normal

When to Use Stabilization

Recommended for:

  • ✅ 3D patient-specific geometries
  • ✅ Complex branching patterns
  • ✅ Cases with flow reversal
  • ✅ High Reynolds number flows

Not needed for:

  • 2D simplified geometries
  • Steady-state simulations
  • Cases without backflow
  • Well-conditioned meshes

Practical Implementation Guide

Step 1: Case Setup

// system/controlDict
application     pimpleFoam;  // or pisoFoam
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         1.0;         // One cardiac cycle
deltaT          0.0001;      // Small timestep for stability
writeControl    adjustableRunTime;
writeInterval   0.01;

libs
(
    "libmodularWKPressure.so"  // Load WK library
);

Step 2: Multiple Outlets

// Different parameters for different outlets
outlet_desc
{
    type            modularWKPressure;
    phi             phi;
    order           2;
    R               8.0e7;    // Descending aorta
    C               1.5e-9;
    Z               1.2e7;
    p0              12000;    // 90 mmHg
    value           uniform 12000;
    q_1             0;
}

outlet_bca
{
    type            modularWKPressure;
    phi             phi;
    order           2;
    R               1.5e8;    // Brachiocephalic (higher resistance)
    C               8.0e-10;  // Lower compliance
    Z               2.0e7;
    p0              12000;
    value           uniform 12000;
    q_1             0;
}

Step 3: Solver Settings

// system/fvSolution
PIMPLE
{
    nOuterCorrectors    3;
    nCorrectors         2;
    nNonOrthogonalCorrectors 1;

    // Important for stability with WK
    momentumPredictor   true;
    consistent          true;

    residualControl
    {
        U               1e-6;
        p               1e-6;
    }
}

Validation and Verification

Pressure Validation

Target Ranges (Healthy Adult) (clinical literature):

Location Systolic Diastolic Mean Pulse Pressure
Ascending Aorta 120±10 mmHg 80±5 mmHg 90±8 mmHg 40±8 mmHg
Arch Vessels 115±10 mmHg 75±5 mmHg 85±8 mmHg 40±8 mmHg
Descending Aorta 110±10 mmHg 70±5 mmHg 80±8 mmHg 40±8 mmHg

Flow Validation

Expected Flow Patterns: - Peak flow: 400-600 mL/s (systole) - Mean flow: 80-100 mL/s (cardiac output) - Flow reversal: < 10% of forward flow - Flow distribution: Murray's law scaling

Expected Flow Distribution

// Expected flow distribution (%)
scalar flow_desc = 65.0;      // Descending aorta
scalar flow_bca = 15.0;       // Brachiocephalic
scalar flow_lcca = 12.0;      // Left carotid
scalar flow_lsca = 8.0;       // Left subclavian

Troubleshooting Common Issues

Problem 1: Simulation Crashes/Instability

Symptoms

  • Timestep drops below 1e-8
  • Continuity errors explode (> 1e10)
  • Solution diverges rapidly
  • "Maximum number of iterations exceeded"

Solutions

// 1. Enable stabilization
outlet
{
    type                stabilizedWindkesselVelocity;
    beta                1.5;  // Increase if still unstable
    enableStabilization true;
}

// 2. Adjust solver settings
PIMPLE
{
    nOuterCorrectors        5;     // Increase iterations
    nCorrectors             3;
    nNonOrthogonalCorrectors 2;

    residualControl
    {
        U               1e-7;  // Tighter convergence
        p               1e-7;
    }
}

// 3. Conservative startup
order           1;      // Start with 1st order
// Gradually increase to 2nd/3rd order after stabilization

Problem 2: Unphysiological Pressure Values

High Pressures (> 20 kPa)

// Reduce resistance
R               5.0e7;     // Lower than 8.0e7

// Increase compliance
C               2.5e-9;    // Higher than 1.5e-9

// Check inlet conditions
inlet
{
    type            flowRateInletVelocity;
    volumetricFlowRate  4.0e-5;  // Reduce if too high (5 L/min)
}

Low Pressures (< 8 kPa)

// Increase resistance
R               1.2e8;     // Higher than 8.0e7

// Decrease compliance
C               8.0e-10;   // Lower than 1.5e-9

// Verify inlet flow rate
// Should be ~5 L/min (8.33e-5 m³/s) for resting conditions

Problem 3: Flow Reversal Issues

Excessive Backflow

// Strengthen stabilization
beta                2.0;   // Increase from 1.0

// Alternative: Use outflow boundary for severe cases
outlet
{
    type            outletPhaseMeanVelocity;
    Umean           0.1;   // Prevent strong backflow
    alpha           alpha.water;
}

No Backflow (Unphysical)

// Reduce stabilization
beta                0.5;   // Decrease from 1.0

// Or disable for testing
enableStabilization false;

// Check mesh quality at outlets
checkMesh -allGeometry -constant

Clinical Applications & Use Cases

Based on Westerhof et al. (2009), Windkessel models have extensive clinical applications:

Cardiovascular Risk Assessment

  • Pulse pressure prediction: Major predictor of cardiovascular mortality
  • Arterial stiffness quantification: Essential for aging and hypertension studies
  • Ventriculo-arterial coupling: Understanding heart-vessel interactions

Therapeutic Device Development

  • Ventricular assist devices: Windkessel load testing for artificial hearts
  • Artificial valve testing: Physiological pressure-flow conditions
  • Cardiac output monitoring: Non-invasive estimation from pressure waveforms

Disease Modeling

  • Hypertension studies: Separating resistance vs. compliance contributions
  • Aging effects: Age-related compliance changes (factor 2-4 decrease)
  • Pulmonary circulation: Right heart afterload assessment

CFD Boundary Condition Applications

  • Patient-specific simulations: Realistic outlet boundary conditions
  • Surgical planning: Pre/post-operative hemodynamic assessment
  • Device optimization: Stent, graft, and bypass performance evaluation

Key Clinical Insight

"For similar percentage changes in total arterial compliance and peripheral resistance, the resistance contributes a factor 4 more to pressure." - Westerhof et al. (2009)

However, with aging, resistance may increase 10-20% while compliance decreases 2-4 fold, making both parameters clinically crucial.

Quick Reference

Key Implementation Steps

  1. Install OpenFOAM-WK boundary condition library
  2. Estimate parameters from clinical data or literature
  3. Configure pressure BC with modularWKPressure
  4. Add velocity stabilization for complex geometries
  5. Validate results against physiological ranges
  6. Monitor stability and adjust β coefficient as needed

References

Primary Reference: Westerhof, N., Lankhaar, J. W., & Westerhof, B. E. (2009). The arterial Windkessel. Medical & Biological Engineering & Computing, 47(2), 131-141. DOI: 10.1007/s11517-008-0359-2

Additional References: ¹ Esmaily Moghadam, M., Bazilevs, Y., Hsia, T. Y., Vignon-Clementel, I. E., & Marsden, A. L. (2011). A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations. Computational Mechanics, 48(3), 277-291.

Resources

Found an issue or have a suggestion for this page?

Open a GitHub issue