Immersed Boundary Method (IBM)

Immersed Boundary Method (IBM)#

Starting from version 5.0, DALES has a grid-conforming Immersed Boundary Method (IBM) implemented. This implementation is a modified version of an older implementation (see Tomas, 2016). This section explains the IBM implementation, how to prepare a run with immersed boundaries, and how to start it.

General description of the IBM#

The Immersed Boundary Method implemented in DALES is grid conforming, meaning that a computational grid cell is either completely set as free fluid or solid object. This effectively changes complex shapes (e.g., inclined roofs) to a cubical representation, which may be innaccurate when the resolution is much coarser than the size of such shapes. The benefit is however its ease of implementation, incl. control over the fluxes and velocity values. Boundaries between fluid and solid always coincide with a cubic face at which the corresponding velocity component is set to zero, thereby ensuring no advection over such boundaries occurs. Momentum and scalar values adjacent to solid cells are corrected by locally adapting the tendencies in two steps:

regular contributions to the diffusive flux over faces (as if no solid were present) are substracted;
new contributions as a result of wall drag or wall flux are calculated and added.

No correction for subgrid TKE is applied as diffusive flux contribution over walls is assumed to be negligble (which may require further testing). Tendencies for values inside solids are every substep set to correct for any drift, for example, for temperature thlp(i,j,k) = (thlibm - thl0(i,j,k))*rk3coefi.

<< still add sketch >>

Currently, boundary conditions for solid walls are still limited. They are listed below:

temperature: Dirichlet condition (fixed temperature) for or zero heat flux across vertical walls and roofs. These values are set to separate constant values for walls thlwall and roofs thlroof, allowing some thermal response between fluid and solid.
moisture: No flux across walls and roofs.
scalars: No flux across walls and roofs.

Preparing the IBM file#

Due to the grid conformity of the IBM, only obstacle heights are required to use the IBM. That implies overhanging structures are not allowed, meaning no fluid cell can be positioned vertically below a solid cell. Obstacle heights should be supplied as a list of values in a file named ibm.inp.<iexpnr>, preceded by seven header lines (see example below). These values are ordered in a row-major ordering, with the first value corresponding to the top-left corner of the domain (xmin,ymax), the second value being one position to the right (xmin+1*dx,ymax), and the last value corresponding to the bottom-right corner (xmax,ymin).

# EXPERIMENT: ibm example
# Horizontal grid points Nx =  256, Ny =  256
# Horizontal grid size Delta x =   5.00, Delta y =   5.00
#
#
#
#
   0.0
   0.0
   5.7
   5.7
   6.8
   5.0
   0.0
   .
   .

Although the structure of this input file is straightforward, creating it requires some steps. Data of the desired obstacles first need to be acquired and then processed to get the correct input file.

Using Wavefront OBJ format#

The following example assumes that your buildings are represented in the Wavefront OBJ format. Building representation in OBJ format can be obtained directly via https://3dbag.nl for the Netherlands or be created via the open-source software City4CFD. Note that the second option is preferred when there is also topography involved and a high degree of accuracy is required. Further, the example assumes that the coordinates are specified in the Dutch RD coordinate system (Rijksdriehoeksmeting). Extension to different coordinates is trivial as long as an orthogonal coordinate system is used.

##### IMPORT PACKAGES #####

import numpy as np
import matplotlib.pyplot as plt
import math 
import trimesh # <-- this package is required to handle OBJ files!

##### SPECIFY GENERAL SETTINGS #####

iexpnr = '001' 
building_file = 'Buildings.obj'
terrain_file  = 'Terrain.obj'

dx = 5.     # horizontal grid spacing in meters
Nx = 256    # number of grid cells in x-direction
Ny = 256    # number of grid cells in y-direction
bz = 8      # number of blend-out cells around edge, where no buildings will be placed

origin_obj  = [139751, 455487]  # RD-COORDINATES of the center of domain (as in OBJ file)
                                # this example is situated on the Utrecht University campus
                                # alternatively, if these are lower-left of the OBJ file, 
                                # origin_shift below should be set to zero
float_type = np.float64

##### PROCESSING DATA #####

origin_shift = [-dx*Nx/2, -dx*Ny/2] # shift from origin_obj to lower left corner of domain
origin_les   = np.asarray(origin_obj) + np.asarray(origin_shift)

# Import OBJ files containing meshes and shift them
mesh1 = trimesh.load(building_file)

mesh1.vertices[:,0] = mesh1.vertices[:,0] - origin_shift[0]
mesh1.vertices[:,1] = mesh1.vertices[:,1] - origin_shift[1]

mesh2 = trimesh.load(terrain_file)

mesh2.vertices[:,0] = mesh2.vertices[:,0] - origin_shift[0]
mesh2.vertices[:,1] = mesh2.vertices[:,1] - origin_shift[1]

# Combine building and terrain meshes in one mesh
mesh = trimesh.boolean.union([mesh1, mesh2], engine='blender', check_volume=False)

# Voxelize the OBJ-mesh to computational grid spacing dx
voxgrid = mesh.voxelized(dx).hollow()

# Create region of interest for LES
end = [Nx*dx, Ny*dy]
xh_les = np.linspace(dx/2, end[0]-dx/2, Nx) 
yh_les = np.linspace(dy/2, end[1]-dy/2, Ny) 

# Prepare arrays for obstacle heights. Mind the orientation:
#   Voxgrid structure orders data as increasing x, then y (xy-ordering); start bottom-left
#   IBM input file should however start at top-left (largest y, smallest x; row-major) 
#   Column indices 'j' correspond to the x-dimension and row indices 'i' to y-dimension 
#   DALES itself corrects for this again by reading the file upside-down
#   There is potential to correct for this code-wise (which has not been done yet)

data_les = np.zeros((Ny,Nx), dtype=float)
mask_2d  = np.zeros((Ny,Nx), dtype=bool)

for bb in range(voxgrid.points.shape[0]):
    xpos = voxgrid.points[bb,0]
    ypos = voxgrid.points[bb,1]
    zhei = voxgrid.points[bb,2]

    ii_pos = Ny - int(np.floor(ypos/dy))
    jj_pos = int(np.floor(xpos/dx))

    if(bz<=jj_pos<Ny-bz and bz<=ii_pos<Ny-bz):
        data_les[ii_pos,jj_pos] = zhei 
        mask_2d [ii_pos,jj_pos] = (zhei >= dx)

print('Maximum height in area: {:>3.1f}'.format(np.max(data_les)))

##### WRITE IBM INPUT FILE #####
do = data_les.ravel() # first meshgridded and then converted to list in right direction

f = open('ibm.inp.'+str(iexpnr),'w')

# Ensure 7 lines of text heade
f.write('# EXPERIMENT: ibm example)
f.write('# Horizontal grid points Nx = {:>4d}, Ny = {:>4d}\n'.format(Nx,Ny))
f.write('# Horizontal grid size Delta x = {:>6.2f}, Delta y = {:>6.2f}\n'.format(dx,dx))
f.write('#\n')
f.write('#\n')
f.write('#\n')
f.write('#\n')

for nn in range(do.shape[0]):
    line = '{:>6.1f}\n'.format(do[nn])
    f.write(line)

f.close() 

##### OPTIONAL: PLOT VOXELIZED BUILDINGS #####
x_mesh, y_mesh = np.meshgrid(xh_les,yh_les,indexing='xy') 
y_mesh         = np.flipud(y_mesh)

xp = x_mesh[mask_2d].ravel()   # like flatten, make continuous vector of xdata
yp = y_mesh[mask_2d].ravel()   
zp = data_les[mask_2d].ravel() 
b0 = np.zeros(zp.shape[0])

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.bar3d(xp, yp, b0, dx, dy, zp, shade="True")
ax.set_xlim3d(origin_les[0],origin_les[0]+dx*Nx)
ax.set_ylim3d(origin_les[1],origin_les[1]+dx*Ny)
ax.set_aspect('equal')
plt.axis('off')
ax.view_init(elev=30, azim=-0)

plt.show()

Alternative: directly from LiDAR#

Sometimes OBJ files are not readily available and creating them yourself with City4CFD is not feasible. Luckily, many authorities nowadays have made accurate LiDAR scans (point clouds) publicly available. For example, tiled LiDAR scans for the Netherlands (AHN) can accessed via GeoTiles and tiled scans for France can be accessed via IGN France. Such LiDAR scans can be quickly voxelized for use with our IBM. Further, most LiDAR return points have a classification, which makes it possible to separate buildings from ground points. The table below shows the typical classification used by the Dutch AHN scans (which may vary among different institutes).

Classification value	Meaning
`1`	unclassified (mostly vegetation)
`2`	ground points
`6`	building points
`9`	water points
`26`	reserved for ASPRS definition

The example below shows the basics of processing such LiDAR scans. Note that the size of the domain is inferred from the dimensions of the LiDAR tile (assumed to be in meters). If you want only part of the domain, you will have to subsample data_les yourself. Likewise, if you want to simulate a larger region, you have to merge LiDAR tiles yourself.

Thanks to Mahmoud Ahmed (TUDelft) for sharing his Python script.

##### IMPORT PACKAGES #####

import math
import matplotlib.pyplot as plt
import open3d as o3d
import numpy as np
import laspy as lp      # <-- this package should be installed with laz backend
                        # pip3 install "laspy[lazrs,laszip]"

##### SPECIFY INPUT FILE AND VOXEL SIZE (=GRID SIZE) #####

file = '37EN2_11.LAZ'   # AHN3 Wippolder area in Delft, downloaded from GeoTiles 
dx   = 10.

##### PROCESS LIDAR POINT CLOUDS #####

point_cloud = lp.read(file)

# Optional: only extract one or two classes. In Dutch AHN, buildings correspond to class 6
point_cloud = point_cloud.points[point_cloud.classification == 6]

pcd = o3d.geometry.PointCloud()

pcd.points = o3d.utility.Vector3dVector(
    np.vstack((point_cloud.x, point_cloud.y, point_cloud.z)).transpose())

pcd.colors = o3d.utility.Vector3dVector(
    np.vstack((point_cloud.red, point_cloud.green, point_cloud.blue))
    .transpose()/255)

# Create voxels from point cloud
voxel_grid  = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd, voxel_size=dx)
voxels      = voxel_grid.get_voxels()

# Determine number of voxels in x- and y-directions
Lx = voxel_grid.get_oriented_bounding_box().extent[0]
Ly = voxel_grid.get_oriented_bounding_box().extent[1]

Nx = math.floor(Lx/dx) if ( math.floor(Lx/dx) % 2 == 0 ) else math.floor(Lx/dx)-1
Ny = math.floor(Ly/dx) if ( math.floor(Ly/dx) % 2 == 0 ) else math.floor(Ly/dx)-1

# Prepare arrays for obstacle heights and coordinate meshes
data_les = np.zeros((Ny,Nx), dtype=float)
mask_2d  = np.zeros((Ny,Nx), dtype=bool)

xh_les = np.linspace(dx/2, Nx*dx-dx/2, Nx) 
yh_les = np.linspace(dx/2, Ny*dx-dx/2, Ny) 

for v in voxels:
    if (v.grid_index[0] >= Nx or v.grid_index[1] >= Ny): continue
    ii_pos = (Ny - 1) - v.grid_index[1]
    jj_pos = v.grid_index[0] 
    zhei   = (v.grid_index[2] + 1) * dx # dx is also voxel size

    data_les[ii_pos,jj_pos] = max(data_les[ii_pos,jj_pos], zhei)
    mask_2d [ii_pos,jj_pos] = True

##### OPTIONAL: PLOT VOXELIZED POINT CLOUD #####

# Visualize voxels
o3d.visualization.draw_geometries([voxel_grid])

##### OPTIONAL: EXPORT THE MESH AS .PLY #####
# vox_mesh = o3d.geometry.TriangleMesh()

# for v in voxels:
#     cube = o3d.geometry.TriangleMesh.create_box(width=1, height=1, depth=1)
#     cube.paint_uniform_color(v.color)
#     cube.translate(v.grid_index, relative=False)
#     vox_mesh += cube

# # Export the mesh as an .obj or .ply file
# vox_mesh.translate([0.5, 0.5, 0.5], relative=True)
# vox_mesh.scale(dx, [0, 0, 0])
# vox_mesh.translate(voxel_grid.origin, relative=True)

# o3d.io.write_triangle_mesh(file+"voxel_mesh_h.ply", vox_mesh)

From this point, writing of the IBM input file is the same as in the earlier example.

Setting &NAMIBM options#

Now the input file containing the obstacle heights has been prepared, you are almost ready to use the IBM. The final step is to switch on the IBM and supply general settings in the &NAMIBM section of the namoptions.<iexpnr> file. The available settings are show in table below.

Options	Default	Allowed Values	Description / Remark
`lapply_ibm`	`.false.`	`.true.` or `.false.`	switch to enable immersed boundary method
`lwallheat`	`.false.`	`.true.` or `.false.`	switch to apply lateral heat flux from buildings
`lpoislast`	`.true.`	`.true.` or `.false.`	switch to apply Poisson solver after the IBM
`thlwall`	\(293\,\mathrm{K}\)	\(>0\)	potential temperature of vertical walls (only if `lwallheat=.true.`)
`thlroof`	\(293\,\mathrm{K}\)	\(>0\)	potential temperature of roof faces
`thlibm`	\(293\,\mathrm{K}\)	\(>0\)	potential temperature inside of obstacles
`qtibm`	\(0\,\mathrm{kg/kg}\)	\(>0\)	moisture inside of obstacles
`z0m_wall`	\(0.03\,\mathrm{m}\)	\(>0\) and \(<\) `zf(1)`	roughness length for momentum at walls/roofs
`z0h_wall`	\(0.03\,\mathrm{m}\)	\(>0\) and \(<\) `zf(1)`	roughness length for heat at walls/roofs

At runtime, DALES will check for potentially conflicting options from different modules and throw ERROR messages when required.

TODO: How to visualize IBM runs#

…