Reproducing Tutorial 3.2


I am trying to re-produce lift-drag graphs obtained in Tutorial 3.2, that will be a baseline for further lift/drag calculations. I simply used the code given in the Tutorial (without dx`s but with SymbolicBFI, I am using an old version) but the lift graph I obtain is slightly different than the one in the tutorial.

Since this will be my baseline, I am confused, which graph I should follow as a reference.Why are they different ?

I have attached the code I use here, (it is ready to run) and the picture of lift graph I obtain, the difference can be seen clearly.

Thanks in advance.

[attachment=undefined]Fi.png[/attachment] [attachment=undefined]Fi.png[/attachment]

from netgen import gui
from ngsolve import *
from netgen.geom2d import SplineGeometry

geo = SplineGeometry()
geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc="cyl")
mesh = Mesh( geo.GenerateMesh(maxh=0.08))
mesh.Curve(3); Draw(mesh)
# viscosity
nu = 0.001

k = 3
V = VectorH1(mesh,order=k, dirichlet="wall|cyl|inlet")
Q = H1(mesh,order=k-1)
X = FESpace([V,Q])

gfu = GridFunction(X)
velocity = gfu.components[0]
from ngsolve.internal import visoptions
visoptions.scalfunction = "u:0"

# parabolic inflow at bc=1:
uin = CoefficientFunction((U*4*y*(0.41-y)/(0.41*0.41),0))
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

import matplotlib 
import numpy as np 
import matplotlib.pyplot as plt
#matplotlib inline

s = np.arange(0.0, 0.42, 0.01)
bvs = U*4*s*(0.41-s)/(0.41)**2
plt.plot(s, bvs)
plt.xlabel('y'); plt.ylabel('$u_x(0,y)$'); plt.title('inflow profile')

(u,p), (v,q) = X.TnT()

a = BilinearForm(X)
stokes = nu*InnerProduct(grad(u),grad(v))-div(u)*q-div(v)*p

a += SymbolicBFI(stokes)

f = LinearForm(X)

inv_stokes = a.mat.Inverse(X.FreeDofs())

res = f.vec.CreateVector() = f.vec - a.mat*gfu.vec += inv_stokes * res


dt = 0.001

# matrix for implicit part of IMEX(1) scheme:
mstar = BilinearForm(X)
mstar += SymbolicBFI(InnerProduct(u,v) + dt*stokes)

inv = mstar.mat.Inverse(X.FreeDofs())

conv = LinearForm(X)
conv += SymbolicLFI(InnerProduct(grad(velocity)*velocity,v))

t = 0
tend = 0

#### For Lift Drag Part ####
drag_x_test = GridFunction(X)
drag_x_test.components[0].Set(CoefficientFunction((-2/(U_mean**2*L),0)), definedon=mesh.Boundaries("cyl"))
drag_y_test = GridFunction(X)
drag_y_test.components[0].Set(CoefficientFunction((0,-2/(U_mean**2*L))), definedon=mesh.Boundaries("cyl"))

time_vals = []
drag_x_vals = []
drag_y_vals = []

# implicit Euler/explicit Euler splitting method:
tend += 1.0
while t < tend-0.5*dt:
    print ("\rt=", t, end="")

    conv.Assemble() = a.mat * gfu.vec + conv.vec -= dt * inv * res
    time_vals.append( t )
    drag_x_vals.append(InnerProduct(res, drag_x_test.vec) )
    drag_y_vals.append(InnerProduct(res, drag_y_test.vec) )

    t = t + dt

# Plot drag force over time
plt.plot(time_vals, drag_x_vals)
plt.xlabel('time'); plt.ylabel('drag'); plt.title('drag'); plt.grid(True)

# Plot lift force over time
plt.plot(time_vals, drag_y_vals)
plt.xlabel('time'); plt.ylabel('lift'); plt.title('lift'); plt.grid(True)


when I run your script, the drag/lift pictures look very similar to the ones in the iTutorial. The small differences are probably due to slightly different meshes, created my newer (or older) versions of netgen. Since this is a relatively coarse mesh and the solution is not fully periodic after t=1, neither could be seen as a “reference” .

Best wishes, Henry