I am running a 2D airfoil simulation in Modulus 24.01 and I am observing the flow inside my airfoil.
I have created a rectangle and subtracted a line (created by naca_airfoil.py from the examples) from it.
height = 5.0
width = 5.0
x, y = Symbol("x"), Symbol("y")
recOrg = Rectangle((-width / 2, -height / 2), (width / 2, height / 2))
# naca profile
geoNaca = NacaGeometry()
# geoNaca = geoNaca.geo.sample_boundary(nr_points=100000)
rec = recOrg - geoNaca
When I sample the inside the samples look correct but in the results the flow appears inside the airfoil.
What am I doing wrong?
BCs:
sampled interior:
flow:
my code:
# SPDX-FileCopyrightText: Copyright (c) 2023 - 2024 NVIDIA CORPORATION & AFFILIATES.
# SPDX-FileCopyrightText: All rights reserved.
# SPDX-License-Identifier: Apache-2.0
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import os
import warnings
import torch
import numpy as np
from sympy import Symbol, sqrt, Max, Eq, Abs
import modulus.sym
from modulus.sym.hydra import to_absolute_path, instantiate_arch, ModulusConfig
from modulus.sym.solver import Solver
from modulus.sym.domain import Domain
from modulus.sym.geometry.primitives_2d import Rectangle
from modulus.sym.domain.constraint import (
PointwiseBoundaryConstraint,
PointwiseInteriorConstraint,
IntegralBoundaryConstraint,
)
from modulus.sym.domain.inferencer import PointwiseInferencer, PointVTKInferencer
from modulus.sym.domain.validator import PointwiseValidator
from modulus.sym.domain.monitor import PointwiseMonitor
from modulus.sym.key import Key
from modulus.sym.eq.pdes.navier_stokes import NavierStokes
from modulus.sym.eq.pdes.basic import NormalDotVec
from modulus.sym.utils.io import (
csv_to_dict,
ValidatorPlotter,
InferencerPlotter,
)
from modulus.sym.geometry.tessellation import Tessellation
from modulus.sym.utils.io.vtk import var_to_polyvtk, VTKUniformGrid
from modulus.sym.hydra import instantiate_arch
from modulus.sym.eq.pdes.turbulence_zero_eq import ZeroEquation
from naca_airfoil import *
@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# make geometry
height = 5.0
width = 5.0
x, y = Symbol("x"), Symbol("y")
recOrg = Rectangle((-width / 2, -height / 2), (width / 2, height / 2))
# naca profile
geoNaca = NacaGeometry()
# geoNaca = geoNaca.geo.sample_boundary(nr_points=100000)
rec = recOrg - geoNaca
s = rec.sample_interior(nr_points=100000)
var_to_polyvtk(s, "rec_geo")
# params
nu = 0.00001488
# nu=0.002
# inlet_vel = 115.129
u = 10
v = 0
w = 0
volumetric_flow = 0.05415
# make aneurysm domain
domain = Domain()
# make list of nodes to unroll graph on
# ns = NavierStokes(nu, rho=1.5, dim=3, time=False)
# model turbulencji + NavierStokes z turbulencjami
ze = ZeroEquation(nu=nu, dim=2, time=False, max_distance=0.1)
ns = NavierStokes(nu=ze.equations["nu"], rho=1.5, dim=2, time=False)
normal_dot_vel = NormalDotVec(["u", "v"])
# normal_dot_vel = NormalDotVec()
flow_net = instantiate_arch(
input_keys=[Key("x"), Key("y")],
output_keys=[Key("u"), Key("v"), Key("p")],
cfg=cfg.arch.fully_connected,
# cfg=cfg.arch.siren,
)
nodes = (
ns.make_nodes()
+ normal_dot_vel.make_nodes()
+ ze.make_nodes()
+ [flow_net.make_node(name="flow_network")]
# + ze.make_nodes()
)
# inlet
inlet = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
criteria=Eq(x, -width / 2),
outvar={"u": u, "v": v},
batch_size=cfg.batch_size.inlet,
lambda_weighting={
"u": 1.0,
"v": 1.0,
},
# batch_per_epoch=5,
)
domain.add_constraint(inlet, "inlet")
# outlet
outlet = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=rec,
criteria=Eq(x, width / 2),
outvar={"u": u, "v": v},
# outvar={},
batch_size=cfg.batch_size.outlet,
lambda_weighting={
"u": 1.0,
"v": 1.0,
},
# batch_per_epoch=5,
)
domain.add_constraint(outlet, "outlet")
# bottomWall - hithub
bottomWall = PointwiseBoundaryConstraint(
nodes=nodes,
outvar={"u": u, "v": v},
geometry=rec,
batch_size=cfg.batch_size.BottomWall,
lambda_weighting={"u": 1.0, "v": 1.0}, # weight edges to be zero
criteria=Eq(y, -height / 2),
)
domain.add_constraint(bottomWall, name="BottomWall")
# topWall - hithub
topWall = PointwiseBoundaryConstraint(
nodes=nodes,
outvar={"u": u, "v": v},
geometry=rec,
batch_size=cfg.batch_size.TopWall,
lambda_weighting={"u": 1.0, "v": 1.0}, # weight edges to be zero
criteria=Eq(y, height / 2),
)
domain.add_constraint(topWall, name="topWall")
# no slip
no_slip = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geoNaca,
outvar={"u": 0, "v": 0},
batch_size=cfg.batch_size.no_slip,
lambda_weighting={
"u": 1.0,
"v": 1.0,
},
# batch_per_epoch=5,
)
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,
lambda_weighting={
"continuity": Symbol("sdf"),
"momentum_x": Symbol("sdf"),
"momentum_y": Symbol("sdf"),
},
compute_sdf_derivatives=True
# batch_per_epoch=5,
)
domain.add_constraint(interior, "interior")
# add monitors
# metric for mass and momentum imbalance
global_monitor = PointwiseMonitor(
rec.sample_interior(1024),
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,
requires_grad=True,
)
domain.add_monitor(global_monitor)
# metric for force on inner sphere
force_monitor = PointwiseMonitor(
geoNaca.sample_boundary(1024),
output_names=["p"],
metrics={
"force_x": lambda var: torch.sum(var["normal_x"] * var["area"] * var["p"]),
"force_y": lambda var: torch.sum(var["normal_y"] * var["area"] * var["p"]),
},
nodes=nodes,
)
domain.add_monitor(force_monitor)
vtk_obj = VTKUniformGrid(
bounds=[[-width / 2, width / 2], [-height / 2, height / 2]],
npoints=[1024, 1024],
export_map={"U": ["u", "v", None], "p": ["p"]},
)
grid_inference = PointVTKInferencer(
vtk_obj=vtk_obj,
nodes=nodes,
input_vtk_map={"x": "x", "y": "y"},
output_names=["u", "v", "p"],
requires_grad=False,
batch_size=1024,
)
domain.add_inferencer(grid_inference, "vtk_inf")
# make solver
slv = Solver(cfg, domain)
# start solver
slv.solve()
if __name__ == "__main__":
run()
naca_airfoil.py
# SPDX-FileCopyrightText: Copyright (c) 2023 - 2024 NVIDIA CORPORATION & AFFILIATES.
# SPDX-FileCopyrightText: All rights reserved.
# SPDX-License-Identifier: Apache-2.0
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import math
import matplotlib.pyplot as plt
import numpy as np
from sympy import Number, Symbol, Heaviside, atan, sin, cos, sqrt
import os
from modulus.sym.geometry.primitives_2d import Polygon
from modulus.sym.geometry.parameterization import Parameterization, Parameter
from modulus.sym.utils.io.vtk import var_to_polyvtk
# Naca implementation modified from https://stackoverflow.com/questions/31815041/plotting-a-naca-4-series-airfoil
# https://en.wikipedia.org/wiki/NACA_airfoil#Equation_for_a_cambered_4-digit_NACA_airfoil
def camber_line(x, m, p, c):
cl = []
for xi in x:
cond_1 = Heaviside(xi, 0) * Heaviside((c * p) - xi, 0)
cond_2 = Heaviside(-xi, 0) + Heaviside(xi - (c * p), 0)
v_1 = m * (xi / p**2) * (2.0 * p - (xi / c))
v_2 = m * ((c - xi) / (1 - p) ** 2) * (1.0 + (xi / c) - 2.0 * p)
cl.append(cond_1 * v_1 + cond_2 * v_2)
return cl
def dyc_over_dx(x, m, p, c):
dd = []
for xi in x:
cond_1 = Heaviside(xi) * Heaviside((c * p) - xi)
cond_2 = Heaviside(-xi) + Heaviside(xi - (c * p))
v_1 = ((2.0 * m) / p**2) * (p - xi / c)
v_2 = (2.0 * m) / (1 - p**2) * (p - xi / c)
dd.append(atan(cond_1 * v_1 + cond_2 * v_2))
return dd
def thickness(x, t, c):
th = []
for xi in x:
term1 = 0.2969 * (sqrt(xi / c))
term2 = -0.1260 * (xi / c)
term3 = -0.3516 * (xi / c) ** 2
term4 = 0.2843 * (xi / c) ** 3
term5 = -0.1015 * (xi / c) ** 4
th.append(5 * t * c * (term1 + term2 + term3 + term4 + term5))
return th
def naca4(x, m, p, t, c=1):
th = dyc_over_dx(x, m, p, c)
yt = thickness(x, t, c)
yc = camber_line(x, m, p, c)
line = []
for xi, thi, yti, yci in zip(x, th, yt, yc):
line.append((xi - yti * sin(thi), yci + yti * cos(thi)))
x.reverse()
th.reverse()
yt.reverse()
yc.reverse()
for xi, thi, yti, yci in zip(x, th, yt, yc):
line.append((xi + yti * sin(thi), yci - yti * cos(thi)))
return line
def NacaGeometry():
# make parameters for naca airfoil
m = 0.02
p = 0.4
t = 0.12
c = 1.0
# make naca geometry
x = [x for x in np.linspace(0, 0.2, 10)] + [x for x in np.linspace(0.2, 1.0, 10)][
1:
] # higher res in front
line = naca4(x, m, p, t, c)[:-1]
geo = Polygon(line)
# sample different parameters
s = geo.sample_boundary(nr_points=100000)
var_to_polyvtk(s, "naca_boundary")
s = geo.sample_interior(nr_points=100000)
var_to_polyvtk(s, "naca_interior")
return geo