Turbulent physics: Zero Equation Turbulence Model

Introduction

This tutorial walks you through the process of adding a algebraic (zero equation) turbulence model to the Modulus simulations. In this tutorial you will learn the following:

  1. How to implement the Zero equation turbulence model in Modulus.

  2. How to create nodes in the graph for arbitrary variables.

Note

This tutorial assumes that you have completed the Lid Driven Cavity Background tutorial on Lid Driven Cavity Flow and have familiarized yourself with the basics of the Modulus APIs.

Problem Description

In this tutorial you will add the zero equation turbulence for a Lid Driven Cavity flow. The problem setup is very similar to the one found in the Lid Driven Cavity Background. The Reynolds number is increased to 1000. The velocity profile is kept the same as the previous problem. To increase the Reynolds Number, the viscosity is reduced to 1 × 10−4 \(m^2/s\).

Case Setup

The case set up in this tutorial is very similar to the example in Lid Driven Cavity Background. It only describes the additions that are made to the previous code.

Note

The python script for this problem can be found at examples/ldc/ldc_2d_zeroEq.py

Importing the required packages

Import Modulus’ ZeroEquation and Node to help setup the problem. Other import are very similar to previous LDC.

from sympy import Symbol, Eq, Abs
import torch
import modulus
from modulus.hydra import to_absolute_path, to_yaml, instantiate_arch
from modulus.hydra.config import ModulusConfig
from modulus.csv_utils.csv_rw import csv_to_dict
from modulus.continuous.solvers.solver import Solver
from modulus.continuous.domain.domain import Domain
from modulus.geometry.csg.csg_2d import Rectangle

from modulus.continuous.constraints.constraint import (
    PointwiseBoundaryConstraint,
    PointwiseInteriorConstraint,
)
from modulus.continuous.monitor.monitor import PointwiseMonitor
from modulus.continuous.validator.validator import PointwiseValidator
from modulus.continuous.inferencer.inferencer import PointwiseInferencer
from modulus.key import Key
from modulus.PDES.navier_stokes import NavierStokes
from modulus.PDES.turbulence_zero_eq import ZeroEquation
from modulus.tensorboard_utils.plotter import ValidatorPlotter, InferencerPlotter
from modulus.architecture import layers
from modulus.node import Node

Defining the Equations, Networks and Nodes

In addition to the Navier Stokes equation, the Zero Equation turbulence model is included by instantiating the ZeroEquation equation class. The kinematic viscosity \(\nu\) is now a string variable in the Navier Stokes equation. The Zero equation turbulence model provides the effective viscosity \((\nu+\nu_t)\) to the Navier Stokes equations. The kinematic viscosity of the fluid calculated based on the Reynolds number is given as an input to the ZeroEquation class.

The Zero Equation turbulence model is defined in the equations below. Note, \(\mu_t = \rho \nu_t\).

(134)\[\mu_t=\rho l_m^2 \sqrt{G}\]
(135)\[G=2(u_x)^2 + 2(v_y)^2 + 2(w_z)^2 + (u_y + v_x)^2 + (u_z + w_x)^2 + (v_z + w_y)^2\]
(136)\[l_m=\min (0.419d, 0.09d_{max})\]

Where, \(l_m\) is the mixing length, \(d\) is the normal distance from wall, \(d_{max}\) is maximum normal distance and \(\sqrt{G}\) is the modulus of mean rate of strain tensor.

The zero equation turbulence model requires normal distance from no slip walls to compute the turbulent viscosity. This information is supplied in the form of variable 'normal_distance' by making a node from the SDF function. node_from_sympy can be used to create such nodes for arbitrary symbolic variables/functions.

The inputs and the outputs of the neural network are specified and the nodes of the architecture are made as before.

@modulus.main(config_path="conf_zeroEq", config_name="config")
def run(cfg: ModulusConfig) -> None:
    print(to_yaml(cfg))

    # add constraints to solver
    # make geometry
    height = 0.1
    width = 0.1
    x, y = Symbol("x"), Symbol("y")
    rec = Rectangle((-width / 2, -height / 2), (width / 2, height / 2))
    
    # make list of nodes to unroll graph on
    ns = NavierStokes(nu='nu', rho=1.0, dim=2, time=False)
    ze = ZeroEquation(nu=1e-4, dim=2, time=False, max_distance=height/2)
    flow_net = instantiate_arch(
        input_keys=[Key("x"), Key("y")],
        output_keys=[Key("u"), Key("v"), Key("p")],
        cfg=cfg.arch.fully_connected,
    )
    nodes = (ns.make_nodes() 
            + ze.make_nodes()
            + [Node.from_sympy(rec.sdf, "normal_distance")]
            + [flow_net.make_node(name="flow_network", jit=cfg.jit)])

Setting up domain, adding constraints and running the solver

This section of the code is very similar to LDC tutorial, so the code and final results is presented here.

    # make ldc domain
    ldc_domain = Domain()

    # top wall
    top_wall = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=rec,
        outvar={"u": 1.5, "v": 0},
        batch_size=cfg.batch_size.TopWall,
        lambda_weighting={"u": 1.0 - 20 * Abs(x), "v": 1.0},  # weight edges to be zero
        criteria=Eq(y, height / 2),
    )
    ldc_domain.add_constraint(top_wall, "top_wall")

    # no slip
    no_slip = PointwiseBoundaryConstraint(
        nodes=nodes,
        geometry=rec,
        outvar={"u": 0, "v": 0},
        batch_size=cfg.batch_size.NoSlip,
        criteria=y < height / 2,
    )
    ldc_domain.add_constraint(no_slip, "no_slip")

    # interior
    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=rec,
        outvar={"continuity": 0, "momentum_x": 0, "momentum_y": 0},
        batch_size=cfg.batch_size.Interior,
        bounds={x: (-width / 2, width / 2), y: (-height / 2, height / 2)},
        lambda_weighting={
            "continuity": rec.sdf,
            "momentum_x": rec.sdf,
            "momentum_y": rec.sdf,
        },
    )
    ldc_domain.add_constraint(interior, "interior")

    # add validator
    mapping = {"Points:0": "x", "Points:1": "y", "U:0": "u", "U:1": "v", "p": "p", "nuT": "nu"}
    openfoam_var = csv_to_dict(
        to_absolute_path("openfoam/cavity_uniformVel_zeroEqn_refined.csv"), mapping
    )
    openfoam_var["x"] += -width / 2  # center OpenFoam data
    openfoam_var["y"] += -height / 2  # center OpenFoam data
    openfoam_var["nu"] += 1e-4 # effective viscosity
    openfoam_invar_numpy = {
        key: value for key, value in openfoam_var.items() if key in ["x", "y"]
    }
    openfoam_outvar_numpy = {
        key: value for key, value in openfoam_var.items() if key in ["u", "v", "nu"]
    }
    openfoam_validator = PointwiseValidator(
        openfoam_invar_numpy,
        openfoam_outvar_numpy,
        nodes,
        batch_size=1024,
        plotter=ValidatorPlotter(),
    )
    ldc_domain.add_validator(openfoam_validator)

    # add inferencer data
    grid_inference = PointwiseInferencer(
        openfoam_invar_numpy,
        ["u", "v", "p", "nu"],
        nodes,
        batch_size=1024,
        plotter=InferencerPlotter(),
    )
    ldc_domain.add_inferencer(grid_inference, "inf_data")
    
    # add monitors
    global_monitor = PointwiseMonitor(
        rec.sample_interior(4000, bounds={x: (-width/2, width/2), y: (-height/2, height/2)}),
        output_names=["continuity", "momentum_x", "momentum_y"],
        metrics={
            "mass_imbalance": lambda var: torch.sum(
                var["area"] * torch.abs(var["continuity"])
            ),
            "momentum_imbalance": lambda var: torch.sum(
                var["area"]
                * (torch.abs(var["momentum_x"]) + torch.abs(var["momentum_y"]))
            ),
        },
        nodes=nodes,
    )
    ldc_domain.add_monitor(global_monitor)

    # make solver
    slv = Solver(cfg, ldc_domain)

    # start solver
    slv.solve()


if __name__ == "__main__":
    run()

The results of the turbulent lid driven cavity flow are shown below.

Visualizing variables from Inference domain

Fig. 47 Visualizing variables from Inference domain

Comparison with OpenFOAM data. Left: Modulus Prediction. Center:OpenFOAM, Right: Difference

Fig. 48 Comparison with OpenFOAM data. Left: Modulus Prediction. Center: OpenFOAM, Right: Difference