LES with Terrain

In this walkthrough, we discuss the steps to setup a terrain simulation using the newly implemented immersed boundary forcing method (IBFM). The theory for the technique can be found at this link.

The setup for the terrain follows the typical simulation of the atmospheric boundary layer (ABL) using large eddy simulation or Reynolds-averaged Navier Stokes turbulence models. The IBFM can be used with periodic or inflow-outflow boundary conditions with few modifications.

The first step in including the terrain is to set the terrain variables. This is accomplished by modifying the ABL physics to include the TerrainDrag flow physics: incflo.physics = ABL TerrainDrag. This looks for the terrain.amrwind text file in the case folder (this is the default name, the user can modify the file it searches for by specifying the TerrainDrag.terrain_file input parameter). The file contains the terrain height as a single column organized as: nx, ny, x values (of length nx), y values (of length ny), terrain height values (of length nx x ny).

The second step is the inclusion of the terrain forcing in the momentum and energy equations. This is accomplished by adding DragForcing and DragTempForcing terms to ICNS.source_terms and Temperature.source_terms, respectively. The terrain simulations requires adding a sponge layer at the outflow and Rayleigh damping at top of the domain. Rayleigh damping is already available from the existing forcing terms and can be used directly. The sponge-layer is implemented by specifying the boundary and the span. For example, a sponge layer of size 1000 m at the east (+x) boundary, we need to include DragForcing.sponge_east=1 and DragForcing.sponge_distance_east=1000 in the input file. The sponge layer is not required for periodic boundary boundary conditions. The only input recommended for the energy equation source term is the specification of the internal temperature of the terrain. This is set as DragTempForcing.reference_temperature=300. It is recommended that the reference_temperature be set to the value of the reference temperature used in the Boussinesq term. The current terrain setup can only be used for the simulation of neutral ABL. A future release will update this calculation to automatically use the values from a precursor simulation for both neutral and non-neutral stratification. The terrain can be visualized by including io.int_outputs = terrain_blank in the input file.

It is recommended to use the ProbeSampler to create the terrain-aligned output planes. The easiest method to generate the text file for ProbeSampler is to write the STL as a text file and then use offsets in postprocessing to write the planes at different heights above the terrain. The terrain-aware output can also be used with FAST.Farm and FLORIS.

An example paraview visualization of the terrain is shown below (with three levels of refinement):

../_images/terrain_blank.png

Here is a sample content of precursor and inflow-outflow input files to drive terrain simulations:

 1# Generating the precursor file
 2# Geometry
 3geometry.prob_lo 			 = 708751 5.00187e+06 446.15 
 4geometry.prob_hi 			 = 723151 5.016e+06 2041.56 
 5geometry.is_periodic 			 = 1 1 0
 6# Grid 
 7amr.n_cell 			 = 232 224 40
 8amr.max_level 			 = 0
 9time.stop_time 			 = -1
10time.max_step 			 = 10000
11time.initial_dt 			 = 0.1
12time.fixed_dt 			 = -1
13time.cfl 			 = 0.9
14time.plot_interval 			 = 5000
15time.checkpoint_interval 			 = 2000
16# incflo 
17incflo.physics 			 = ABL
18incflo.density 			 = 1.225
19incflo.gravity 			 = 0.  0. -9.81  # Gravitational force (3D)
20incflo.velocity 			 = 10 0 0 
21incflo.verbose  			 = 0
22incflo.initial_iterations 			 = 8
23incflo.do_initial_proj 			 = true
24incflo.constant_density 			 = true
25incflo.use_godunov 			 = true
26incflo.godunov_type 			 = "weno_z"
27incflo.diffusion_type 			 = 2
28# transport equation parameters 
29transport.model 			 = ConstTransport
30transport.viscosity 			 = 1e-5
31transport.laminar_prandtl 			 = 0.7
32transport.turbulent_prandtl 			 = 0.333
33# turbulence equation parameters 
34turbulence.model 			 = Kosovic
35Kosovic.refMOL 			 = -1e30
36# Atmospheric boundary layer
37ABL.Uperiods 			 = 72
38ABL.Vperiods 			 = 72
39ABL.cutoff_height 			 = 50.0
40ABL.deltaU 			 = 1.0
41ABL.deltaV 			 = 1.0
42ABL.perturb_ref_height 			 = 50.0
43ABL.perturb_velocity 			 = true
44ABL.perturb_temperature 			 = false
45ABL.kappa 			 = .41
46ABL.normal_direction 			 = 2
47ABL.reference_temperature 			 = 300
48ABL.stats_output_format 			 = netcdf
49ABL.surface_roughness_z0 			 = 0.1
50ABL.temperature_heights 			 = 0  800 900 1900
51ABL.temperature_values  			 = 300 300 308 311
52ABL.wall_shear_stress_type 			 = local
53ABL.surface_temp_flux 			 = 0
54ABL.bndry_file 			 = "bndry_files"
55ABL.bndry_write_frequency 			 = 100
56ABL.bndry_io_mode 			 = 0
57ABL.bndry_planes 			 = xlo ylo 
58ABL.bndry_output_start_time 			 = 434.028
59ABL.bndry_var_names 			 = velocity temperature
60ABL.bndry_output_format 			 = native
61# Source
62ICNS.source_terms 			 = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing RayleighDamping NonLinearSGSTerm
63BoussinesqBuoyancy.reference_temperature 			 = 300
64BoussinesqBuoyancy.thermal_expansion_coeff 			 = 0.00333333
65CoriolisForcing.east_vector 			 = 1.0 0.0 0.0 
66CoriolisForcing.north_vector 			 = 0.0 1.0 0.0 
67CoriolisForcing.latitude 			 = 90 
68CoriolisForcing.rotational_time_period 			 = 125664 
69GeostrophicForcing.geostrophic_wind 			 = 10 0 0
70RayleighDamping.reference_velocity 			 = 10 0 0
71RayleighDamping.length_sloped_damping 			 = 400
72RayleighDamping.length_complete_damping 			 = 200
73RayleighDamping.time_scale 			 = 5.0
74# BC 
75zhi.type 			 = "slip_wall"
76zhi.temperature_type 			 = "fixed_gradient"
77zhi.temperature 			 =  0.003
78zlo.type 			 = "wall_model"
79mac_proj.mg_rtol 			 = 1.0e-4 
80mac_proj.mg_atol 			 = 1.0e-8 
81mac_proj.maxiter 			 = 360 
82nodal_proj.mg_rtol 			 = 1.0e-4 
83nodal_proj.mg_atol 			 = 1.0e-8 
84diffusion.mg_rtol 			 = 1.0e-4 
85diffusion.mg_atol 			 = 1.0e-8 
86temperature_diffusion.mg_rtol 			 = 1.0e-4 
87temperature_diffusion.mg_atol 			 = 1.0e-8 
88nodal_proj.maxiter 			 = 360 

 1# Generating the terrain file
 2# Geometry
 3geometry.prob_lo 			 = 708751 5.00187e+06 446.15 
 4geometry.prob_hi 			 = 723151 5.016e+06 2041.56 
 5geometry.is_periodic 			 = 0 0 0
 6# Grid 
 7amr.n_cell 			 = 232 224 40
 8amr.max_level 			 = 0
 9time.stop_time 			 = -1
10time.max_step 			 = 10000
11time.initial_dt 			 = 0.1
12time.fixed_dt 			 = -1
13time.cfl 			 = 0.9
14time.plot_interval 			 = 5000
15time.checkpoint_interval 			 = 2000
16# incflo 
17incflo.physics 			 = ABL TerrainDrag
18incflo.density 			 = 1.225
19incflo.gravity 			 = 0.  0. -9.81  # Gravitational force (3D)
20incflo.velocity 			 = 10 0 0 
21incflo.verbose  			 = 0
22incflo.initial_iterations 			 = 8
23incflo.do_initial_proj 			 = true
24incflo.constant_density 			 = true
25incflo.use_godunov 			 = true
26incflo.godunov_type 			 = "weno_z"
27incflo.diffusion_type 			 = 2
28# transport equation parameters 
29transport.model 			 = ConstTransport
30transport.viscosity 			 = 1e-5
31transport.laminar_prandtl 			 = 0.7
32transport.turbulent_prandtl 			 = 0.333
33# turbulence equation parameters 
34turbulence.model 			 = Kosovic
35Kosovic.refMOL 			 = -1e30
36# Atmospheric boundary layer
37ABL.kappa 			 = .41
38ABL.normal_direction 			 = 2
39ABL.reference_temperature 			 = 300
40ABL.stats_output_format 			 = netcdf
41ABL.surface_roughness_z0 			 = 0.1
42ABL.temperature_heights 			 = 0  800 900 1900
43ABL.temperature_values  			 = 300 300 308 311
44ABL.wall_shear_stress_type 			 = local
45ABL.surface_temp_flux 			 = 0
46ABL.bndry_file 			 = "../precursor/bndry_files"
47ABL.bndry_io_mode 			 = 1
48ABL.bndry_var_names 			 = velocity temperature
49ABL.bndry_output_format 			 = native
50# Source
51ICNS.source_terms 			 = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing RayleighDamping NonLinearSGSTerm DragForcing
52BoussinesqBuoyancy.reference_temperature 			 = 300
53BoussinesqBuoyancy.thermal_expansion_coeff 			 = 0.00333333
54CoriolisForcing.east_vector 			 = 1.0 0.0 0.0 
55CoriolisForcing.north_vector 			 = 0.0 1.0 0.0 
56CoriolisForcing.latitude 			 = 90 
57CoriolisForcing.rotational_time_period 			 = 125664 
58GeostrophicForcing.geostrophic_wind 			 = 10 0 0
59RayleighDamping.reference_velocity 			 = 10 0 0
60RayleighDamping.length_sloped_damping 			 = 400
61RayleighDamping.length_complete_damping 			 = 200
62RayleighDamping.time_scale 			 = 5.0
63# BC 
64xlo.type 			 = "mass_inflow"
65xlo.density 			 = 1.225
66xlo.temperature 			 = 300
67xhi.type 			 = "pressure_outflow"
68ylo.type 			 = "mass_inflow"
69ylo.density 			 = 1.225
70ylo.temperature 			 = 300
71yhi.type 			 = "pressure_outflow"
72zhi.type 			 = "slip_wall"
73zhi.temperature_type 			 = "fixed_gradient"
74zhi.temperature 			 =  0.003
75zlo.type 			 = "wall_model"
76mac_proj.mg_rtol 			 = 1.0e-4 
77mac_proj.mg_atol 			 = 1.0e-8 
78mac_proj.maxiter 			 = 360 
79nodal_proj.mg_rtol 			 = 1.0e-4 
80nodal_proj.mg_atol 			 = 1.0e-8 
81diffusion.mg_rtol 			 = 1.0e-4 
82diffusion.mg_atol 			 = 1.0e-8 
83temperature_diffusion.mg_rtol 			 = 1.0e-4 
84temperature_diffusion.mg_atol 			 = 1.0e-8 
85nodal_proj.maxiter 			 = 360 
86#io 
87io.restart_file 			 = "../precursor/chk02000"

Setup using Python Tools

The setup of the terrain files can be cumbersome to do by hand. A set of python tools are made available at amrTerrain. A more comprehensive set of tools will be available in future at: windtools.

The python code is executed as follows:

python backendinterface.py nameofyamlfile.yaml

Sample input files are available in the GitHub repository. A typical sample file looks as follows:

 1solver: "amrWind"
 2caseParent: "/Users/hgopalan/Documents/P101_AMR-Wind/Data/tempGUI"
 3caseFolder: "WFIP2_test"
 4caseType: "terrainTurbine"
 5caseInitial: "amr"
 6centerLat: 45.63374
 7centerLon: -120.66047
 8refHeight: 2184
 9west: 5000
10east: 5000
11south: 5000
12north: 5000
13cellSize: 128
14verticalAR: 4
15timeMethod: "step"
16numOfSteps: 5000
17plotOutput: 1000
18restartOutput: 1000
19forcingHeight: 10.0
20windX: 13.5
21windY: 0.0 
22windZ: 0.0 
23refTemperature: 300.0
24refRoughness: 0.1
25refHeatflux: 0.0 
26refLat: 45.63374
27refPeriod: 125663.706143592
28includeCoriolis: True
29turbineMarkType: "database"
30turbineType: "UniformCtDisk"

The variable caseType takes three kinds of inputs: precursor or terrain or terrainTurbine. For running terrain simulations, it is recommended to use caseType:terrain. The use of caseType:terrainTurbine also creates turbines aligned with the terrain height using the turbine latitude and longitude in the file turbine.csv.

The python code first reads the centerLat and centerLon and creates a domain of size specified by west, east, south, and north. For the example shown above, a domain size of 10 km is created around centerLat and centerLon. The terrain module uses the SRTM 30 m database to create the terrain. It is possible to add a user-defined file to define the terrain by modifying the python code.

The cellSize: 128 sets a grid resolution of 128 m at level 0.The variable verticalAR: 4 sets dz=4dx=4dy. You do not need Hypre to run the high aspect ratio simulations. User has to manually edit the input file to create refinement regions in area of interest around the terrain.

All other inputs in the yaml file are for creating dummy inputs to the amr-wind simulations and user can modify them manually to fit their needs. The inputs caseType: "terrainTurbine" and turbineType: "UniformCtDisk" are useful for aligning the turbine vertically with the terrain. The file includes all the turbines within the continental US and have to be modified for other locations. The turbine type information is ad-hoc and has to be manually modified by the user for the specific turbine type. A future update to the code will include options to specify the turbine information from a text file.