Fluid Dynamics Predictions Within Helical Geometry (Inlet / Outlet Problem) Are Inaccurate

Hello Everyone,

I have designed a 3D helical geometry which has some girth, using NVIDIA Modulus. Defined the inlet and outlet as two very thin cylinders at the top and bottom of the helix, added integral continuity throughout the right side of the helix which consists of several thin cylinders, then finally, the helix has walls (no slip) and an interior space plus some sample discs throughout the geometry. I am trying to simulate water flowing through this geometry from inlet to outlet; top to bottom.

You can see the images above:

  • Left: Shows the helix walls, interior, inlet (top disc), and outlet (bottom disc).

  • Middle: Shows the same as left with the added integral continuity discs that occur at the RHS of the helix.

  • Right: Shows the same as the middle but with the added sample discs at both the LHS and RHS of the helix.

I’ve given the PINN constraints at this different boundaries such as:

Inlet

u, v, w = circular_parabola(
    Symbol("x"),
    Symbol("y"),
    Symbol("z"),
    center=(r + R, 0, n_turns*pitch),
    normal=(0, -1, 0),
    radius=r,
    max_vel=inlet_vel,
)

inlet = PointwiseBoundaryConstraint(
    nodes=nodes,
    geometry=inlet,
    outvar={"u": u, "v": v, "w": w},
    batch_size=cfg.batch_size.inlet,
)
domain.add_constraint(inlet, "inlet")

Which makes the inlet velocity follow a parabolic motion with higher velocity at the center of the cylinder and 0 at the walls.

Outlet

outlet = PointwiseBoundaryConstraint(
    nodes=nodes,
    geometry=outlet,
    outvar={"p": 0},
    batch_size=cfg.batch_size.outlet,
)
domain.add_constraint(outlet, "outlet")

Specified pressure will be 0 at the outlet as is the case with all inlet / outlet flow problems of this caliber.

No Slip

no_slip = PointwiseBoundaryConstraint(
    nodes=nodes,
    geometry=geo,
    outvar={"u": 0, "v": 0, "w": 0},
    batch_size=cfg.batch_size.no_slip,
)
domain.add_constraint(no_slip, "no_slip")

Specified velocity will be 0 at the walls of the helix.

Interior

interior = PointwiseInteriorConstraint(
    nodes=nodes,
    geometry=geo,
    outvar={"continuity": 0, "momentum_x": 0, "momentum_y": 0, "momentum_z": 0},
    batch_size=cfg.batch_size.interior,
    lambda_weighting={
        "continuity": 2 * Symbol("sdf"),
        "momentum_x": 2 * Symbol("sdf"),
        "momentum_y": 2 * Symbol("sdf"),
        "momentum_z": 2 * Symbol("sdf")
    },
)
domain.add_constraint(interior, "interior")

Typical interior constraints, not sure if the lambda_weighting is correct though.

Sample Interior

sample_interior = PointwiseInteriorConstraint(
    nodes=nodes,
    geometry=sample_discs,
    outvar={'continuity': 0, 'momentum_x': 0, 'momentum_y': 0, 'momentum_z': 0},
    lambda_weighting={
        'continuity': 5 * Symbol('sdf'),
        'momentum_x': 2 * Symbol('sdf'),
        'momentum_y': 2 * Symbol('sdf'),
        'momentum_z': 2 * Symbol('sdf'),
    },
    batch_size=cfg.batch_size.interior_sample,
)
domain.add_constraint(sample_interior, 'sample_discs')

I’ve added some sample discs so I can later process velocity and pressure predictions, given the same constraints as the interior of the helix.

Integral Continuity Discs

def int_cont(n, r, R, pitch, inlet_vel):

    integral_disc = Cylinder(center=(r + R, 0, 0), radius=r, height=0.00001).rotate(angle=np.pi/2, axis='x', center=(r + R, 0, 0)).translate((0, 0, n*pitch))

    integral_continuity = IntegralBoundaryConstraint(
        nodes=nodes,
        geometry=integral_disc,
        outvar={"normal_dot_vel": inlet_vel},
        batch_size=1,
        integral_batch_size=cfg.batch_size.integral_continuity,
        lambda_weighting={"normal_dot_vel": 1},
        # parameterization=z_pos_range,
    )
    domain.add_constraint(integral_continuity, "integral_continuity_" + str(n))

    print("Made Intergral Continuity Disc " + str(n))

for n in range(1, n_turns):
    int_cont(n, r, R, pitch, inlet_vel)

Made a function which places integral continuity discs throughout one side of the helix and specifies that the ‘normal_dot_vel’ will be equal to the inlet velocity at each disc.

Geometric Setup

config = [{"R": 0.5, "r": 0.055, "pitch": 0.2, "n_turns": 12, "inlet_vel": 17.538}]

num = 0
R = config[num].get("R")
r = config[num].get("r")
pitch = config[num].get("pitch")
n_turns = config[num].get("n_turns")
inlet_vel = config[num].get("inlet_vel")

# define sympy variables to parametrize domain curves
x, y, z = Symbol("x"), Symbol("y"), Symbol("z")

# define geometry

def Helix(R, r, pitch, n_turns):

    n = 50
    epsilon = 0.00001

    theta = np.linspace(0, 2*n_turns*np.pi, n_turns*n)
    height = np.linspace(0, n_turns*pitch, n_turns*n)

    geo = Sphere(center=((r + R)*np.cos(theta[0]), (r + R)*np.sin(theta[0]), height[0]), radius=r)
    for i, j in zip(theta, height):
        geo += Sphere(center=((r + R)*np.cos(i), (r + R)*np.sin(i), j), radius=r)
        
    sample_discs = Cylinder(center=(r + R, 0, 0), radius=r, height=epsilon).rotate(angle=np.pi/2, axis='x', center=(r + R, 0, 0)).translate((0, 0, pitch))
    for k in range(2, n_turns):
        sample_discs += Cylinder(center=(r + R, 0, 0), radius=r, height=epsilon).rotate(angle=np.pi/2, axis='x', center=(r + R, 0, 0)).translate((0, 0, k*pitch))
    for l in range(1, n_turns + 1):
        sample_discs += Cylinder(center=(r + R, 0, 0), radius=r, height=epsilon).rotate(angle=np.pi/2, axis='x', center=(r + R, 0, 0)).translate((-2*(R + r), 0, l*pitch - 0.5*pitch))

    b0 = Box(((r + R)*np.cos(theta[0]) - r, (r + R)*np.sin(theta[0]) - r, height[0] - r), ((r + R)*np.cos(theta[0]) + r, (r + R)*np.sin(theta[0]), height[0] + r))
    b1 = Box(((r + R)*np.cos(i) - r, (r + R)*np.sin(i), j - r), ((r + R)*np.cos(i) + r, (r + R)*np.sin(i) + r, j + r))
    geo -= (b0 + b1)

    outlet = geo & b0.translate((0, epsilon, 0))
    inlet = geo & b1.translate((0, -epsilon, 0))

    return inlet, outlet, geo, sample_discs

inlet, outlet, geo, sample_discs = Helix(R, r, pitch, n_turns)

sample_plane = Box((-1.1*r, -0.00001, -1.1*r), (1.1*r, 0.00001, n_turns*pitch + 1.1*r)) 

This code basically defines some parameters like:

  • R: Helix external radius
  • r: Helix internal radius
  • pitch: Vertical distance (z-axis) between each turn of the helix
  • n_turns: Number of turns the helix has
  • inlet_vel: Inlet velocity of the helix

Within the Helix() function, I’ve written a script that creates a helix out of (n=50 x n_turns=12) spheres with radius r and places them along the equation of a helix, which you can see being used in the geo = … line.

The code also adds the sample discs here for sampling velocity and pressure after the model has trained.

b0 and b1 are small boxes which are used to clip the ends of the helix after it has been made out of hundreds of spheres so that the ends of the helix are flat and the inlet and outlet can be placed there.

The inlet and outlet are defined as small discs using the cylinder geometry and they’re placed at the top and bottom of the helix respectively.

PDE and Model Info

ns = NavierStokes(nu=0.01, rho=1.0, dim=3, time=False)
    normal_dot_vel = NormalDotVec(["u", "v", "w"])
    flow_net = instantiate_arch(
        input_keys=[Key("x"), Key("y"), Key("z")],
        output_keys=[Key("u"), Key("v"), Key("w"), Key("p")],
        cfg=cfg.arch.fully_connected,
    )
    nodes = (
        ns.make_nodes()
        + normal_dot_vel.make_nodes()
        + [flow_net.make_node(name="flow_network")]
    )

This is how I have initialized the model. Used the Navier Stokes PDEs and given nu=0.01, rho=1.0, dim=3, time=False.

Inputs are: x, y, x.

Outputs are: u, v, w, p.

Specified this is a fully connected PINN.

Problem

The velocity and pressure predictions are miles off and I have tried many things to fix this issue like scaling all of the parameters, my nu and rho values used to be 1e-6 and 1000.0 for example, now they have been altered so that they’re less orders of magnitude apart. I did this with the geometric parameters too. Results were different but they’re far away from the validation still.

Does anyone have any clue how to fix this? Help would be massively appreciated because I have officially ran out of ideas. Also please let me know if you’d like additional information about the problem.

Many thanks,
Thomas