# cylinderflow.py¶

In this script we solve the Navier-Stokes equations around a cylinder, using the same Raviart-Thomas discretization as in drivencavity-compatible.py but in curvilinear coordinates. The mesh is constructed such that all elements are shape similar, growing exponentially with radius such that the artificial exterior boundary is placed at large (configurable) distance.

 10 import nutils, numpy 

The main function defines the parameter space for the script. Configurable parameters are the mesh density (in number of elements along the cylinder wall), polynomial degree, Reynolds number, rotational velocity of the cylinder, time step, exterior radius, and the end time of the simulation (by default it will run forever).

 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 def main(nelems: 'number of elements' = 24, degree: 'polynomial degree' = 3, reynolds: 'reynolds number' = 1000., rotation: 'cylinder rotation speed' = 0., timestep: 'time step' = 1/24, maxradius: 'approximate domain size' = 25., seed: 'random seed' = 0, endtime: 'end time' = numpy.inf): elemangle = 2 * numpy.pi / nelems melems = int(numpy.log(2*maxradius) / elemangle + .5) nutils.log.info('creating {}x{} mesh, outer radius {:.2f}'.format(melems, nelems, .5*numpy.exp(elemangle*melems))) domain, geom = nutils.mesh.rectilinear([melems, nelems], periodic=(1,)) domain = domain.withboundary(inner='left', outer='right') ns = nutils.function.Namespace() ns.uinf = 1, 0 ns.r = .5 * nutils.function.exp(elemangle * geom[0]) ns.Re = reynolds ns.phi = geom[1] * elemangle # add small angle to break element symmetry ns.x_i = 'r _i' ns.J = ns.x.grad(geom) ns.unbasis, ns.utbasis, ns.pbasis = nutils.function.chain([ # compatible spaces domain.basis('spline', degree=(degree,degree-1), removedofs=((0,),None)), domain.basis('spline', degree=(degree-1,degree)), domain.basis('spline', degree=degree-1), ]) / nutils.function.determinant(ns.J) ns.ubasis_ni = 'unbasis_n J_i0 + utbasis_n J_i1' # piola transformation ns.u_i = 'ubasis_ni ?lhs_n' ns.p = 'pbasis_n ?lhs_n' ns.sigma_ij = '(u_i,j + u_j,i) / Re - p δ_ij' ns.h = .5 * elemangle ns.N = 5 * degree / ns.h ns.rotation = rotation ns.uwall_i = '0.5 rotation <-sin(phi), cos(phi)>_i' inflow = domain.boundary['outer'].select(-ns.uinf.dotnorm(ns.x), ischeme='gauss1') # upstream half of the exterior boundary sqr = inflow.integral('(u_i - uinf_i) (u_i - uinf_i)' @ ns, degree=degree*2) cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15) # constrain inflow semicircle to uinf sqr = domain.integral('(u_i - uinf_i) (u_i - uinf_i) + p^2' @ ns, degree=degree*2) lhs0 = nutils.solver.optimize('lhs', sqr) # set initial condition to u=uinf, p=0 numpy.random.seed(seed) lhs0 *= numpy.random.normal(1, .1, lhs0.shape) # add small velocity noise res = domain.integral('(ubasis_ni u_i,j u_j + ubasis_ni,j sigma_ij + pbasis_n u_k,k) d:x' @ ns, degree=9) res += domain.boundary['inner'].integral('(N ubasis_ni - (ubasis_ni,j + ubasis_nj,i) n_j) (u_i - uwall_i) d:x / Re' @ ns, degree=9) inertia = domain.integral('ubasis_ni u_i d:x' @ ns, degree=9) bbox = numpy.array([[-2,46/9],[-2,2]]) # bounding box for figure based on 16x9 aspect ratio bezier0 = domain.sample('bezier', 5) bezier = bezier0.subset((bezier0.eval((ns.x-bbox[:,0]) * (bbox[:,1]-ns.x)) > 0).all(axis=1)) interpolate = nutils.util.tri_interpolator(bezier.tri, bezier.eval(ns.x), mergetol=1e-5) # interpolator for quivers spacing = .05 # initial quiver spacing xgrd = nutils.util.regularize(bbox, spacing) for istep, lhs in nutils.log.enumerate('timestep', nutils.solver.impliciteuler('lhs', residual=res, inertia=inertia, lhs0=lhs0, timestep=timestep, constrain=cons, newtontol=1e-10)): t = istep * timestep x, u, normu, p = bezier.eval(['x_i', 'u_i', 'sqrt(u_k u_k)', 'p'] @ ns, lhs=lhs) ugrd = interpolate[xgrd](u) with nutils.export.mplfigure('flow.png', figsize=(12.8,7.2)) as fig: ax = fig.add_axes([0,0,1,1], yticks=[], xticks=[], frame_on=False, xlim=bbox[0], ylim=bbox[1]) im = ax.tripcolor(x[:,0], x[:,1], bezier.tri, p, shading='gouraud', cmap='jet') import matplotlib.collections ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='k', linewidths=.1, alpha=.5)) ax.quiver(xgrd[:,0], xgrd[:,1], ugrd[:,0], ugrd[:,1], angles='xy', width=1e-3, headwidth=3e3, headlength=5e3, headaxislength=2e3, zorder=9, alpha=.5) ax.plot(0, 0, 'k', marker=(3,2,t*rotation*180/numpy.pi-90), markersize=20) cax = fig.add_axes([0.9, 0.1, 0.01, 0.8]) cax.tick_params(labelsize='large') fig.colorbar(im, cax=cax) if t >= endtime: break xgrd = nutils.util.regularize(bbox, spacing, xgrd + ugrd * timestep) return lhs0, lhs 

If the script is executed (as opposed to imported), nutils.cli.run() calls the main function with arguments provided from the command line.

 102 103 if __name__ == '__main__': nutils.cli.run(main) 

Once a simulation is developed and tested, it is good practice to save a few strategic return values for regression testing. The nutils.testing module, which builds on the standard unittest framework, facilitates this by providing nutils.testing.TestCase.assertAlmostEqual64() for the embedding of desired results as compressed base64 data.

 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 class test(nutils.testing.TestCase): @nutils.testing.requires('matplotlib', 'scipy') def test_rot0(self): lhs0, lhs = main(nelems=6, reynolds=100, timestep=.1, endtime=.05, rotation=0) with self.subTest('initial condition'): self.assertAlmostEqual64(lhs0, ''' eNqtjD8OwWAcQJ/JNSQ20Tbf135RkUjEZO8RJA7gChYXsDgEkZjN+k/zbQYDCU06Y2Co3yG86S3vtb27 C8fiXMDM0Q7s7MHCRHUUpPkqh42eaxhlvQzKQAewTMIEQjM2MEyuMUylrOxDykt3Id63Rrzprj0YFJ8T 7L2vHMlfcqlU6UMrjVJ4+8+gwS3exnUdye8//AB+zDQQ''', atol=2e-13) with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, ''' eNoB2AAn/5A0jjV/MIDKj8rFMoE4Rjcwz4LI7sery545+Dm5MwTGEsa8NVY8pjtWNSzE18OpyXI9VD02 M5zCnsJazE0+Hj76NsPByMH/yhA3izOMyGPIyC+gN5Y4JcofyEbI+csJOGk4OzXZxrTGLzIKOXo7Acj2 xOjENMk3O8Y85DcZwyTDAzjaPFY+sMfJwavBhDNPPvbFV8cxOKk3ADtFOFI86zqjN9o8D8hcNFjCfsXV Pd47Vj/qPdZBa0F5QUZD7UEJQYi527zjROVETUeVRfZIfrfvRKZKs7s6SVXLZ9k=''') @nutils.testing.requires('matplotlib', 'scipy') def test_rot1(self): lhs0, lhs = main(nelems=6, reynolds=100, timestep=.1, endtime=.05, rotation=1) with self.subTest('initial condition'): self.assertAlmostEqual64(lhs0, ''' eNqtjD8OwWAcQJ/JNSQ20Tbf135RkUjEZO8RJA7gChYXsDgEkZjN+k/zbQYDCU06Y2Co3yG86S3vtb27 C8fiXMDM0Q7s7MHCRHUUpPkqh42eaxhlvQzKQAewTMIEQjM2MEyuMUylrOxDykt3Id63Rrzprj0YFJ8T 7L2vHMlfcqlU6UMrjVJ4+8+gwS3exnUdye8//AB+zDQQ''', atol=2e-13) with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, ''' eNoB2AAn/4s0kDW8MIHKjcq1MoE4RzdQz4PI7sely545+Dm6MwTGEsa8NVY8pjtWNSzE18OpyXI9VD02 M5zCnsJazE0+Hj76NsPByMH/yi03ODSmyHbI0jGyN5M4FcoayEHI2MsEOGs4PjXZxrXGXTILOXo7AMj2 xOfEMsk3O8Y85DcZwyTDAzjaPFY+sMfJwavBhDNPPvTFXMc6OK43/zo7OFI87DqpN9o8Dcg2NFfCgcXX Pd87Vj/pPdZBbEF5QUZD7UEIQYe527zjROVETUeVRfZIfrfvRKZKsrs6ScqLaQk=''')