Help calculating the torque of a H-VAWT

[IsTrivial]

I am trying to replicate the results of a real-life turbine to validate the results of the simulation made in FluidX3D. However, I am running into a problem: The torque calculated is way higher than it should be. I am currently using the approach recommended in another closed issue by @ProjectPhysX: "For thrust measurement, it's better to not measure forces on the propeller, but average flux (from velocity field) through an area behind the propeller". I am very lost because the real-life turbine had an average torque of 52 Nm, and I have an average torque of 300 Nm in my simulation. Below, I attach all the information used to build the setup:

The 3-blade H-VAWT should have a diameter of 0.9 m, and each airfoil should have a chord length of 0.13275 m. I generate the .stl for the turbine using the following Python script:

import numpy as np
from stl import mesh
import math

def get_airfoil_points(profile='naca0025', num_points=50):
    """Generate airfoil coordinates (x, y) normalized to chord length 1."""
    if profile.lower() == 'naca0025':
        t = 0.25  # Max thickness as fraction of chord
        x = np.linspace(0, 1, num_points // 2)
        yt = 5 * t * (0.2969 * np.sqrt(x) - 0.1260 * x - 0.3516 * x**2 + 0.2843 * x**3 - 0.1015 * x**4)
        upper = np.column_stack((x, yt))
        lower = np.column_stack((x, -yt))
        points = np.vstack((upper[::-1], lower[1:-1]))
    else:
        points = np.array(profile)
        points[:, 0] /= points[:, 0].max()  # Normalize to chord length 1
    return points

def create_hvawt_stl(
    airfoil_profile='naca0025',
    num_blades=3,
    chord_length=0.13275,
    diameter=0.9,
    height=1.0,
    angle_of_attack=0.0,
    twist_angle=0.0,
    num_sections=50,
    num_airfoil_points=50,
    filename='hvawt.stl'
):
    """
    Generate a binary STL file for a Darrieus VAWT with 3 blades.
    
    Parameters:
    - airfoil_profile: 'naca0025' or list of (x, y) coordinates
    - num_blades: Number of blades (int)
    - chord_length: Chord length of airfoil (m)
    - diameter: Turbine diameter (m)
    - height: Turbine height (m)
    - angle_of_attack: Angle of attack in degrees
    - twist_angle: Total helical twist angle in degrees (0 for Darrieus)
    - num_sections: Number of vertical sections
    - num_airfoil_points: Number of points per airfoil profile
    - filename: Output STL filename
    """
    radius = diameter / 2  # 0.45 m
    aoa_rad = np.radians(angle_of_attack)
    twist_rad = np.radians(twist_angle)
    airfoil = get_airfoil_points(airfoil_profile, num_airfoil_points)
    
    # Scale airfoil to chord length
    airfoil *= chord_length
    print(f"Airfoil points (scaled to chord length {chord_length}): {airfoil[:5]}")  # Debug
    
    # Initialize lists for vertices and faces
    all_vertices = []
    all_faces = []
    
    # Generate each blade independently
    points_per_section = len(airfoil)
    for blade in range(num_blades):
        blade_angle = 2 * np.pi * blade / num_blades  # 120° apart
        blade_vertices = []
        vertex_offset = len(all_vertices)  # Track starting index for this blade
        
        # Generate sections along the height
        for i in range(num_sections + 1):
            z = i * height / num_sections
            theta = twist_rad * i / num_sections  # 0 for Darrieus
            rot_angle = blade_angle + theta
            
            # Airfoil center at radius
            center_x = radius * np.cos(rot_angle)
            center_y = radius * np.sin(rot_angle)
            
            # Rotate airfoil to be tangential to the circular path
            # Chord should be perpendicular to the radius (rotate by rot_angle + 90°)
            tangential_angle = rot_angle + np.pi / 2  # +90° to align chord tangentially
            
            # Rotate and translate airfoil points
            section_vertices = []
            for x, y in airfoil:
                # Step 1: Rotate to align chord tangentially
                x_tan = x * np.cos(tangential_angle) - y * np.sin(tangential_angle)
                y_tan = x * np.sin(tangential_angle) + y * np.cos(tangential_angle)
                
                # Step 2: Apply angle of attack (relative to tangential direction)
                x_rot = x_tan * np.cos(aoa_rad) - y_tan * np.sin(aoa_rad)
                y_rot = x_tan * np.sin(aoa_rad) + y_tan * np.cos(aoa_rad)
                
                # Step 3: Translate to center
                x_pos = center_x + x_rot
                y_pos = center_y + y_rot
                section_vertices.append([x_pos, y_pos, z])
            
            blade_vertices.extend(section_vertices)
            
            # Debug: Print position of first point in first and last section
            if i == 0 or i == num_sections:
                print(f"Blade {blade}, Section {i}, First point: {section_vertices[0]}")
        
        all_vertices.extend(blade_vertices)
        
        # Triangulate the blade surface
        for i in range(num_sections):
            section_start = i * points_per_section
            next_section_start = (i + 1) * points_per_section
            for j in range(points_per_section):
                j_next = (j + 1) % points_per_section
                # Two triangles per quad
                all_faces.append([
                    vertex_offset + section_start + j,
                    vertex_offset + next_section_start + j,
                    vertex_offset + next_section_start + j_next
                ])
                all_faces.append([
                    vertex_offset + section_start + j,
                    vertex_offset + next_section_start + j_next,
                    vertex_offset + section_start + j_next
                ])
    
    # Add top and bottom caps for each blade
    for blade in range(num_blades):
        blade_start = blade * (num_sections + 1) * points_per_section
        
        # Bottom cap (z=0)
        bottom_center = [
            radius * np.cos(2 * np.pi * blade / num_blades),
            radius * np.sin(2 * np.pi * blade / num_blades),
            0
        ]
        bottom_center_idx = len(all_vertices)
        all_vertices.append(bottom_center)
        
        for j in range(points_per_section):
            j_next = (j + 1) % points_per_section
            all_faces.append([
                bottom_center_idx,
                blade_start + j,
                blade_start + j_next
            ])
        
        # Top cap (z=height)
        top_section_start = blade_start + num_sections * points_per_section
        top_angle = 2 * np.pi * blade / num_blades + twist_rad
        top_center = [
            radius * np.cos(top_angle),
            radius * np.sin(top_angle),
            height
        ]
        top_center_idx = len(all_vertices)
        all_vertices.append(top_center)
        
        for j in range(points_per_section):
            j_next = (j + 1) % points_per_section
            all_faces.append([
                top_center_idx,
                top_section_start + j_next,
                top_section_start + j
            ])
    
    # Create STL mesh
    vertices = np.array(all_vertices, dtype=np.float32)
    faces = np.array(all_faces, dtype=np.uint32)
    turbine_mesh = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
    for i, f in enumerate(faces):
        for j in range(3):
            turbine_mesh.vectors[i][j] = vertices[f[j]]
    
    # Save as binary STL
    turbine_mesh.save(filename)
    print(f"STL file saved as {filename}")

# Example usage
if __name__ == "__main__":
    create_hvawt_stl(
        airfoil_profile='naca0025',
        num_blades=3,
        chord_length=0.13275,
        diameter=0.9,
        height=0.7,
        angle_of_attack=0.0,
        twist_angle=0.0,  # Darrieus (no twist)
        num_sections=200,
        num_airfoil_points=200,
        filename='hvawt.stl'
    )

For the setup.cpp, I set a velocity inlet of 1.62 m/s for water and 60 rpm for the turbine angular speed:

#include "defines.hpp"
#include "lbm.hpp"
#include "shapes.hpp"
#include <cmath>
#include <cstddef>
#include <algorithm>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

void main_setup() {
    // ==================== SIMULATION SETUP ====================
    const uint3 lbm_N = resolution(float3(7.2f, 7.2f, 1.8f), 500u); // 8x8x2 diameters
    const float si_length = 7.2f;
    const float si_velocity = 1.62f;
    const float si_density = 1000.0f;
    const float si_nu = 1.0e-6f;
    const float lbm_u = 0.1f;

    units.set_m_kg_s((float)lbm_N.x, lbm_u, 1.0f, si_length, si_velocity, si_density);
    const float lbm_nu = units.nu(si_nu);
    const float dt = units.si_t(1.0f);

    const float si_rpm = 60.0f;
    const float si_omega = 2.0f * M_PI * (si_rpm / 60.0f); // rad/s
    const float radius_si = 0.45f; // 0.9 m diameter
    const float lbm_radius = radius_si / units.si_x(1.0f);
    const float lbm_omega = si_omega * dt;

    const ulong lbm_T = 48000ull;
    const ulong lbm_dt = 10ull;

    printf("=== Conversion Factors ===\n");
    printf("SI length per cell: %.6f m\n", units.si_x(1.0f));
    printf("SI time per step: %.6f s\n", dt);
    printf("SI velocity per LBM u: %.6f m/s\n", units.si_u(1.0f));
    printf("LBM viscosity: %.6e, SI viscosity: %.6e m²/s\n", lbm_nu, units.si_nu(lbm_nu));
    printf("Turbine diameter in LBM units: %.4f, SI units: %.4f m\n", 2.0f * lbm_radius, units.si_x(2.0f * lbm_radius));

    LBM lbm(lbm_N, 1u, 1u, 1u, lbm_nu);

    // ==================== GEOMETRY SETUP ====================
    const float3 center = float3(3.0f * lbm_radius, lbm.center().y, lbm.center().z);
    Mesh* mesh = read_stl(get_exe_path() + "../stl/hvawt.stl", lbm.size(), center, 2.0f * lbm_radius);
    if (mesh == nullptr) {
        printf("Error: Failed to load mesh\n");
        exit(1);
    }

    const uint Nx = lbm.get_Nx(), Ny = lbm.get_Ny(), Nz = lbm.get_Nz();
    parallel_for(lbm.get_N(), [&](ulong n) {
        uint x = 0u, y = 0u, z = 0u; lbm.coordinates(n, x, y, z);
        if (x == 0u) {
            lbm.flags[n] = TYPE_E;
            lbm.u.x[n] = lbm_u;
        } else if (x == Nx - 1u) {
            lbm.flags[n] = TYPE_E;
            lbm.rho[n] = 1.0f;
        } else if (y == 0u || y == Ny - 1u || z == 0u) {
            lbm.flags[n] = TYPE_S;
        } else if (z == Nz - 1u) {
            lbm.flags[n] = TYPE_E;
            lbm.rho[n] = 1.0f;
        }
    });

    // ==================== MAIN SIMULATION LOOP ====================
    lbm.graphics.visualization_modes = VIS_FLAG_LATTICE | VIS_FLAG_SURFACE | VIS_Q_CRITERION;
    lbm.run(0u);

    // Medición del wake a 3 radios aguas abajo
    const uint x_measure = (uint)(center.x + 3.0f * lbm_radius);
    float torque_sum = 0.0f;
    ulong measurement_count = 0ull;

    const float dx = units.si_x(1.0f);
    const float cell_area = dx * dx;
    const float swept_area = 0.9f * 0.7f; // Turbine swept area en m² (diámetro x altura real)

    while (lbm.get_t() < lbm_T) {
        if (mesh == nullptr) {
            printf("Error: Mesh is null at time %lu\n", lbm.get_t());
            exit(1);
        }

        mesh->rotate(float3x3(float3(0.0f, 0.0f, 1.0f), lbm_omega * (float)lbm_dt));
        lbm.voxelize_mesh_on_device(mesh, TYPE_S, center, float3(0.0f), float3(0.0f, 0.0f, lbm_omega));
        lbm.run(lbm_dt);

        if (lbm.get_t() % 100 == 0 && lbm.get_t() >= 2000) {
            lbm.u.read_from_device();
            float u_wake_sum = 0.0f;
            int valid_cells = 0;

            for (uint y = 0; y < Ny; y++) {
                for (uint z = 0; z < Nz; z++) {
                    ulong n = lbm.index(x_measure, y, z);
                    float dy = (y - center.y) * dx;
                    float dz = (z - center.z) * dx;
                    float r = sqrtf(dy * dy + dz * dz);

                    if (r <= 2.0f * radius_si) {
                        float3 u_si = units.si_u(lbm.u[n]);
                        u_wake_sum += u_si.x;
                        valid_cells++;
                    }
                }
            }

            float u_wake_avg = valid_cells > 0 ? u_wake_sum / valid_cells : 0.0f;
            const float momentum_in = si_density * si_velocity * si_velocity * swept_area;
            const float momentum_flux = si_density * u_wake_avg * u_wake_avg * swept_area;
            const float thrust = momentum_in - momentum_flux;
            const float power_extracted = thrust * si_velocity;
            const float torque = power_extracted / si_omega;

            torque_sum += torque;
            measurement_count++;

            printf("Time %lu: Thrust = %.2f N, Power extracted = %.2f W, Torque = %.4f N·m, Avg Torque = %.4f N·m, Wake Velocity = %.4f m/s (Cells: %d)\n",
                   lbm.get_t(), thrust, power_extracted, torque, torque_sum / measurement_count, u_wake_avg, valid_cells);
        }

        #if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)
        if (lbm.graphics.next_frame(lbm_T, 30.0f)) {
            lbm.graphics.set_camera_free(float3(0.5f * Nx, 0.5f * Ny, 2.0f * Nz), 0.0f, 90.0f, 100.0f);
            lbm.graphics.write_frame();
        }
        #endif
    }

    if (mesh != nullptr) {
        delete mesh;
        mesh = nullptr;
    }
}

This is an example of the values I am obtaining by running the simulation with this setup:

| 2348 | 261 GB/s | 246 | 47999 90% | 0s |Time 48000: Thrust = 1158.82 N, Power extracted = 1877.29 W, Torque = 298.7801 N·m, Avg Torque = 298.0444 N·m, Wake Velocity = 0.8860 m/s (Cells: 5580)

This is how the setup looks like:

This is how the velocity field looks like:

I hope you can help me with this problem because I feel really lost trying to solve this problem.

[cTatu]

I think this is a FluidX3D limitation. The Ahmed Body simulation also returns wrong force values and it's because:
"(...) expect drag force to be too large by about a factor 2, due to lack of grid refinement and no wall model. (...) I'm stuck with the literal corner case in 3D where specular reflection of DDFs reflect them from space diagonals to axis-aligned directions, which causes trouble with the different collision operator weights for those directions. There is no solution to this in literature either." - ProjectPhysX