Hi,
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
U=1.5
L=0.1
U_mean=2/3*U
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]
Draw(velocity,mesh,"u",sd=3)
Draw(gfu.components[1],mesh,"p",sd=3)
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"))
Redraw()
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')
plt.grid(True)
plt.show()
(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)
a.Assemble()
f = LinearForm(X)
f.Assemble()
inv_stokes = a.mat.Inverse(X.FreeDofs())
res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res
Redraw()
dt = 0.001
# matrix for implicit part of IMEX(1) scheme:
mstar = BilinearForm(X)
mstar += SymbolicBFI(InnerProduct(u,v) + dt*stokes)
mstar.Assemble()
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()
res.data = a.mat * gfu.vec + conv.vec
gfu.vec.data -= 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) )
#print(drag)
t = t + dt
Redraw()
# Plot drag force over time
plt.plot(time_vals, drag_x_vals)
plt.xlabel('time'); plt.ylabel('drag'); plt.title('drag'); plt.grid(True)
plt.show()
# Plot lift force over time
plt.plot(time_vals, drag_y_vals)
plt.xlabel('time'); plt.ylabel('lift'); plt.title('lift'); plt.grid(True)
plt.show()