Reproducing Tutorial 3.2

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()