''' Channel flow past a cylindrical obstacle, using a LB method Copyright (C) 2006 Jonas Latt Address: Rue General Dufour 24, 1211 Geneva 4, Switzerland E-mail: Jonas.Latt@cui.unige.ch ''' import numpy # GENERAL FLOW CONSTANTS lx = 250 ly = 51 obst_x = lx/5.+1 # position of the cylinder; (exact obst_y = ly/2.+1 # y-symmetry is avoided) obst_r = ly/10.+1 # radius of the cylinder uMax = 0.02 # maximum velocity of Poiseuille inflow Re = 100 # Reynolds number nu = uMax * 2.*obst_r / Re # kinematic viscosity omega = 1. / (3*nu+1./2.) # relaxation parameter maxT = 4000 # total number of iterations tPlot = 5 # cycles # D2Q9 LATTICE CONSTANTS # matlab: t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36]; # matlab: cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1]; # matlab: cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1]; # matlab: opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7]; # matlab: col = [2:(ly-1)]; t = numpy.array([4/9., 1/9.,1/9.,1/9.,1/9., 1/36.,1/36.,1/36.,1/36.]) cx = numpy.array([ 0, 1, 0, -1, 0, 1, -1, -1, 1]) cy = numpy.array([ 0, 0, 1, 0, -1, 1, 1, -1, -1]) opp = numpy.array([ 1, 4, 5, 2, 3, 8, 9, 6, 7]) # CHANGED: col = numpy.arange(2,(ly-1)) # Because matlab includes the last index, first the direct translation # is numpy.arange(2, ly). Because second numpy uses zero-based indices, # and {col} is used for indicing, it's maybe correct to write: col = numpy.arange(1, ly - 1) # matlab: [y,x] = meshgrid(1:ly,1:lx); # matlab: obst = (x-obst_x).^2 + (y-obst_y).^2 <= obst_r.^2; # matlab: obst(:,[1,ly]) = 1; # matlab: bbRegion = find(obst); y,x = numpy.meshgrid(numpy.arange(ly),numpy.arange(lx)) obst = (x-obst_x)**2 + (y-obst_y)**2 <= obst_r**2 obst[:,0] = obst[:,ly-1] = 1 bbRegion = numpy.nonzero(obst) # INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i) # matlab: fIn = reshape( t' * ones(1,lx*ly), 9, lx, ly); # CHANGED: fIn = numpy.ones((lx,ly,9)) # fIn [:,:] = t # Duplicate the t array in x and y dimensions: fIn = t[:, numpy.newaxis, numpy.newaxis].\ repeat(lx, axis = 1).\ repeat(ly, axis = 2) # MAIN LOOP (TIME CYCLES) # matlab: for cycle = 1:maxT for cycle in range(maxT): # MACROSCOPIC VARIABLES # matlab: rho = sum(fIn); # matlab: ux = reshape ( ... # matlab: (cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho; # matlab: uy = reshape ( ... # matlab: (cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho; # CHANGED: rho = fIn.sum(-1) # Sum in "c" dimension: rho = fIn.sum(axis = 0) # CHANGED: ux = (cx*fIn).sum(-1)/rho # Weight with cx array and sum along c axis: ux = (cx[:, numpy.newaxis, numpy.newaxis] * fIn).sum(axis = 0) / rho # CHANGED: uy = (cx*fIn).sum(-1)/rho # Weight with cy array and sum along c axis: uy = (cy[:, numpy.newaxis, numpy.newaxis] * fIn).sum(axis = 0) / rho # MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS # Inlet: Poiseuille profile # matlab: L = ly-2; y = col-1.5; # matlab: ux(:,1,col) = 4 * uMax / (L*L) * (y.*L-y.*y); # matlab: uy(:,1,col) = 0; # matlab: rho(:,1,col) = 1 ./ (1-ux(:,1,col)) .* ( ... # matlab: sum(fIn([1,3,5],1,col)) + ... # matlab: 2*sum(fIn([4,7,8],1,col)) ); L = ly-2.0 y = col-1.5 # Is that right ? # Probably not. # CHANGED: ux[0,1:-1] = 4 * uMax / (L*L) * (y.*L-y.*y) # I'm not shure about the following. It's just a translation of the # matlab code, given the type of the variables and assumed that # indices are set correctly. ux[:, 1, col] = 4 * uMax / (L ** 2) * (y * L - y ** 2) # Is that right ? # CHANGED: uy[0,1:-1] = 0 uy[:, 1, col] = 0 # CHANGED: rho[0,1:-1] = 1./(1-ux[1,1:-1])*(sum(fIn[ ([1,3,5],1,col)) + 2*sum(fIn([4,7,8],1,col))) # TODO I do not understand how one can give to a two-dimensional array # {rho} three indices (matlab code). I leave out the first : instead. rho[0, col] = 1 / (1 - ux[0, col]) * \ (fIn[[1, 3, 5]][:, 0][:, col].sum(axis = 0) + \ 2 * fIn[[4, 7, 8]][:, 1][:, col].sum(axis = 0)) # Outlet: Zero gradient on rho/ux # matlab: rho(:,lx,col) = rho(:,lx-1,col); # matlab: uy(:,lx,col) = 0; # matlab: ux(:,lx,col) = ux(:,lx-1,col); # % Outlet: Zero gradient on rho/ux # rho(:,lx,col) = rho(:,lx-1,col); # uy(:,lx,col) = 0; # ux(:,lx,col) = ux(:,lx-1,col); rho[lx - 1, col] = rho[lx - 2, col] uy[lx - 1, col] = 0 ux[lx - 1, col] = ux[:, lx - 2, col] # % COLLISION STEP # for i=1:9 # cu = 3*(cx(i)*ux+cy(i)*uy); # fEq(i,:,:) = rho .* t(i) .* ( 1 + cu + 1/2*(cu.*cu) - 3/2*(ux.^2+uy.^2) ); # fOut(i,:,:) = fIn(i,:,:) - omega .* (fIn(i,:,:)-fEq(i,:,:)); # end fEq = numpy.zeros((9, lx, ly)) fOut = numpy.zeros((9, lx, ly)) for i in xrange(0, 9): cu = 3 * (cx[i] * ux + cy[i] * uy) fEq[i] = rho * t[i] * (1 + cu + 0.5 * cu ** 2 - \ 1.5 * (ux ** 2 + uy ** 2)) fOut[i] = fIn[i] - omega * (fIn[i] - fIn[i]) # % MICROSCOPIC BOUNDARY CONDITIONS # for i=1:9 # % Left boundary # fOut(i,1,col) = fEq(i,1,col) + 18*t(i)*cx(i)*cy(i)* ( fIn(8,1,col) - fIn(7,1,col)-fEq(8,1,col)+fEq(7,1,col) ); # % Right boundary # fOut(i,lx,col) = fEq(i,lx,col) + 18*t(i)*cx(i)*cy(i)* ( fIn(6,lx,col) - fIn(9,lx,col)-fEq(6,lx,col)+fEq(9,lx,col) ); # % Bounce back region # fOut(i,bbRegion) = fIn(opp(i),bbRegion); # NOTE: Changed indices from matlab to python: for i in xrange(0, 9): # Left boundary: fOut[i, 1, col] = fEq[i, 0, col] + 18 * t[i] * cx[i] * cy[i] * \ (fIn[7, 0, col] - fIn[6, 0, col] - fEq[7, 0, col] + \ fEq[6, 0, col]) fOut[i, lx - 1, col] = \ fEq[i, lx - 1, col] + 18 * t[i] * cx[i] * cy[i] * \ (fIn[5, lx - 1, col] - fIn[8, lx - 1, col] - \ fEq[5, lx - 1, col] + fEq[8, lx - 1, col]) # TODO Not shure about opp(i). fOut[i, bbRegion] = fIn[(9 - 1) - i, bbRegion] # % STREAMING STEP # for i=1:9 # fIn(i,:,:) = circshift(fOut(i,:,:), [0,cx(i),cy(i)]); # TODO Don't know how to implement circshift().