| 314 | | = References = |
| | 319 | == Non-cyclic lateral boundary conditions == |
| | 320 | |
| | 321 | Figure 6 shows the grid structure for non-cyclic boundary conditions at the left/right boundary '''LB/RB''' ([../../app/inipar/#bc_lr bc_lr]) and figure 7 for non-cyclic boundary conditions at the north/south boundary '''NB/SB''' ([../../app/inipar/#bc_ns bc_ns]). |
| | 322 | The indices (''i,j,k'') represent the directions (''x,y,z''). The model domain extends from ''-1:nx+1'' in the ''x''-direction, from ''-1:ny+1'' in the ''y''-direction and from ''0:nzt+1'' in the ''z''-direction. For the advection scheme of Wicker and Skamarock, two more grid points are added at the lateral boundaries which are not needed for non-cyclic boundary conditions. The figures display the grid layer of the horizontal velocity components ''u'' and ''v'', and scalar ''s''. The grid points of the vertical velocity w are defined at the scalar position but shifted by one half grid spacing in vertical direction (not shown, detailed information about the grid structure in PALM can be found [wiki:doc/tec/discret here]). The prognostic equations are solved at all inner grid points which are marked black. The grid points at the respective non-cyclic boundaries (blue) are treated as follows. '''LB''' is defined at ''i = -1'' for ''v'', ''w'', ''s'' and at ''i = 0'' for ''u''. '''SB''' is defined at ''j = -1'' for ''u'', ''w'', ''s'' and at ''j = 0'' for ''v''. '''RB''' is defined at ''i = nx + 1'' and '''NB''' at ''j = ny + 1'' for all quantities. '''LB''' and '''SB''' are treated this way so that the order and number of grid points for the streamwise velocity component and scalars is the same, independent of the flow direction. |
| | 323 | |
| | 324 | |
| | 325 | For technical reasons, the prognostic equations are first solved for ''u'' at ''i = 0'' (''v'' at ''j = 0''), since these grid points technically belong to the inner grid, but afterwards, these results at ''i = 0'' (''j = 0'') are replaced by the respective boundary condition in routine ''boundary_conds.f90''. In case of a Dirichlet condition, the values at ''i = 0'' (''j = 0'') are taken from ''i = -1'' (''j = -1''). In case of a radiation boundary condition, the solution of the Sommerfeld equation overwrites the prognostic values at ''i = 0'' (''j = 0''). |
| | 326 | |
| | 327 | |
| | 328 | For non-cyclic lateral boundary conditions, the parameter [../../app/inipar/#psolver psolver] has to be set to '' 'multigrid' '' because the default FFT-solver can only be applied for cyclic boundary conditions. |
| | 329 | |
| | 330 | [[Image(grid_lr.png, 700px, border=1)]] |
| | 331 | |
| | 332 | Figure 6: Grid structure at the lateral boundaries with non-cyclic lateral boundary conditions along the left-right direction. |
| | 333 | |
| | 334 | [[Image(grid_ns.png, 700px, border=1)]] |
| | 335 | |
| | 336 | Figure 7: Grid structure of the lateral boundaries with non-cyclic lateral boundary conditions along the north-south direction. |
| | 337 | |
| | 338 | |
| | 339 | === Inflow boundary === |
| | 340 | |
| | 341 | At the inflow boundary, Dirichlet conditions are used for the three velocity components ''ψ = {u,v,w}'' as well as for all scalar quantities s and are implemented as follows (here e.g. for ''s'' and a flow in positive ''x'' direction): |
| | 342 | {{{ |
| | 343 | #!Latex |
| | 344 | \begin{align*} |
| | 345 | s^{t + \Delta t}(k,j,-1) = s_{init}(k) \; . |
| | 346 | \end{align*} |
| | 347 | }}} |
| | 348 | ''t'' denotes the time, Δt the time step and ''s'',,init,, the initialization profile of the scalar quantities which is constant in time. |
| | 349 | The quantities at the inflow are set by the initial vertical profiles (see [../../app/inipar/#initializing_actions initializing_actions]). |
| | 350 | A Neumann condition is used for the subgrid-scale turbulent kinetic energy ''e'' (here e.g. for a left-right flow): |
| | 351 | {{{ |
| | 352 | #!Latex |
| | 353 | \begin{align*} |
| | 354 | e^{t + \Delta t}(k,j,-1) = e^{t + \Delta t}(k,j,0) \; . |
| | 355 | \end{align*} |
| | 356 | }}} |
| | 357 | To prevent gravity waves from being reflected at the inflow, a relaxation term can be added to the prognostic equations for the potential temperature ''θ'' ([#davies1976 Davies, 1976]): |
| | 358 | {{{ |
| | 359 | #!Latex |
| | 360 | \begin{align*} |
| | 361 | \theta^{t+1}(d) = ... - \Delta t \cdot K(d) \cdot \left( \theta^{t}(d) - \theta_{init} \right) \; . |
| | 362 | \end{align*} |
| | 363 | }}} |
| | 364 | Here, ''d'' is the distance normal to the wall and ''θ'',,init,, the initial value of the potential temperature, which corresponds to the value at the inflow boundary. The damping or relaxation function ''K'' depends only on the distance ''d'' to the inflow. ''K'' is calculated by |
| | 365 | {{{ |
| | 366 | #!Latex |
| | 367 | \begin{align*} |
| | 368 | K(d) = |
| | 369 | \begin{cases} |
| | 370 | d_f \sin^2\left( \frac{\pi}{2} \frac{d_w - d}{d_w} \right) , \text{for } d < d_w \\ \qquad\quad 0 \qquad \quad \;\;\; , \text{for } d \ge d_w |
| | 371 | \end{cases} . |
| | 372 | \end{align*} |
| | 373 | }}} |
| | 374 | ''d'',,f,, is a damping factor to control the damping intensity, and ''d'',,w,, is the width of the relaxation region extending from the inflow. Quantities ''d'',,f,, and ''d'',,w,, can be set with parameters [../../app/inipar/#pt_damping_factor pt_damping_factor] and [../../app/inipar/#pt_damping_width pt_damping_width], respectively. Both parameters have to be set by the user and must be adjusted case-by-case, because both parameters depend on the numerical and physical conditions, so that application of universal default values is not possible. |
| | 375 | So far, we have experience with gravity waves in case of cold air outbreaks, which grow in amplitude up to quite extreme values, if no damping is applied. In the respective simulations, we used typical values for [../../app/inipar/#pt_damping_factor pt_damping_factor] of 0.05 and for [../../app/inipar/#pt_damping_width pt_damping_width] of 25 km in order to prevent the gravity waves from growing. |
| | 376 | |
| | 377 | |
| | 378 | === Outflow boundary === |
| | 379 | |
| | 380 | At the outflow, an open boundary condition is needed to ensure that disturbances of the mean flow can exit the model domain without effecting the flow upstream. For the scalar quantities, Neumann boundary conditions are used at the outflow boundary which is the simplest way. For the velocity components, a Neumann condition would require to be considered in the solution of the Poisson equation for perturbation pressure, which has not been realized so far, because it requires some technical effort. Instead, PALM offers two types of radiation boundary conditions for the velocity components, which are not in conflict with the pressure solver (see [../../app/inipar/#use_cmax use_cmax], [../../app/inipar/#bc_lr bc_lr] and [../../app/inipar/#bc_ns bc_ns]). For the radiation condition, the Sommerfeld radiation equation is solved at the outflow |
| | 381 | {{{ |
| | 382 | #!Latex |
| | 383 | \begin{align*} |
| | 384 | \partial_t \psi + c_{\psi} \partial_n \psi = 0 \; , |
| | 385 | \end{align*} |
| | 386 | }}} |
| | 387 | which considers flow disturbances propagating with the mean flow and by waves. Here ''ψ'' is the transported quantity and ''∂'',,n,, is the derivative normal to the outflow boundary. In PALM, based on the equation above, the radiation boundary condition is realized in two ways as follows. |
| | 388 | |
| | 389 | |
| | 390 | ==== Variable Phase velocity ==== |
| | 391 | |
| | 392 | The phase velocity ''c'',,ψ,, (also called convective velocity) is calculated after [#orlanski1976 Orlanski (1976)]: |
| | 393 | {{{ |
| | 394 | #!Latex |
| | 395 | \begin{align*} |
| | 396 | c_{\psi} = - \dfrac{\partial_t \psi}{\partial_n \psi} \; . |
| | 397 | \end{align*} |
| | 398 | }}} |
| | 399 | If the outflow is defined at the right boundary (''i = nx + 1'', Left-right flow, see Fig. 6), the phase velocity for each velocity component is calculated by |
| | 400 | {{{ |
| | 401 | #!Latex |
| | 402 | \begin{align*} |
| | 403 | c_{\psi} = - c_{max} \dfrac{\psi^t_{nx} - \psi^{t - \Delta t}_{nx} }{ \psi^{t-\Delta t}_{nx} - \psi^{t - \Delta t}_{nx-1} } \; , |
| | 404 | \end{align*} |
| | 405 | }}} |
| | 406 | with the maximum phase velocity |
| | 407 | {{{ |
| | 408 | #!Latex |
| | 409 | \begin{align*} |
| | 410 | c_{max} = \dfrac{\Delta x}{\Delta t} \; . |
| | 411 | \end{align*} |
| | 412 | }}} |
| | 413 | The phase velocity has to be in the range of ''0 ≤ c,,ψ,, < c'',,max,, because negative values propagate waves back to the inner domain. ''c'',,max,, represents the maximum phase velocity that ensures numerical stability (Courant-Friedrichs-Lewy condition). Both conditions are not ensured by the second equation of c,,ψ,,, hence c,,ψ,, is set to zero for negative values and to ''c'',,max,, for larger values than ''c'',,max,,. The velocity components ψ at the outflow boundary are then calculated by |
| | 414 | {{{ |
| | 415 | #!Latex |
| | 416 | \begin{align*} |
| | 417 | \psi^{t+\Delta t}_{nx+1} = \psi^{t}_{nx+1} - \dfrac{\overline{c}_{\psi}}{c_{max}} (\psi^{t}_{nx+1} - \psi^{t}_{nx}) \; , |
| | 418 | \end{align*} |
| | 419 | }}} |
| | 420 | with the phase velocity averaged parallel to the outflow: |
| | 421 | {{{ |
| | 422 | #!Latex |
| | 423 | \begin{align*} |
| | 424 | \overline{c}_{\psi} = \dfrac{1}{ny+1} \sum_{j=0}^{ny} c_{\psi, j} \; . |
| | 425 | \end{align*} |
| | 426 | }}} |
| | 427 | In Orlanskis work, the phase velocity ''c'',,ψ,, was not averaged along the outflow, which is sufficient for simplified flows as shown by [#yoshida2010 Yoshida and Watanabe (2010)]. However, we found that in simulations of shear driven convective boundary layer with a strong wind component along the outflow, the approach of [#orlanski1976 Orlanski (1976)] was unstable. [#raymond1984 Raymond and Kuo (1984)] argued that locally calculated phase velocities leads to a large variation of phase speeds which can lead to numerical instabilities. Other works used a constant phase velocity or an average over the whole outflow area (see e.g. [#frölich2006 Fröhlich 2006, p.216) instead of the approach of [#orlanski1976 Orlanski (1976)]. With this in mind, we tested the average phase velocity calculated from the last mentioned equation, and the solution at the outflow becomes stable. We did not additionally average the phase velocity along the vertical direction, because the background wind usually varies with height. |
| | 428 | |
| | 429 | For the other flow directions, the indices of the above equations have to be replaced as follows: |
| | 430 | {{{ |
| | 431 | #!Latex |
| | 432 | \begin{tabular}{|c |c |c |c| c|} |
| | 433 | \hline |
| | 434 | & Right-left flow & \multicolumn{2}{c|}{ North-south flow} & South-north flow \\ |
| | 435 | \hline |
| | 436 | & $(nx + 1) \rightarrow 0$ & $(nx + 1) \rightarrow -1$ & & \\ |
| | 437 | $\psi = u$ & $ nx \rightarrow 1$ & $ nx \rightarrow 0$ & & \\ |
| | 438 | & $(nx - 1) \rightarrow 2$ & $(nx - 1) \rightarrow 1$ & & \\ |
| | 439 | \cline{1-3} |
| | 440 | & $(nx + 1) \rightarrow -1$ & $(nx + 1) \rightarrow 0$ & & $i \rightarrow j$ \\ |
| | 441 | $\psi = v$ & $ nx \rightarrow 0$ & $ nx \rightarrow 1$ & $i \rightarrow j$ & \\ |
| | 442 | & $(nx - 1) \rightarrow 1$ & $(nx - 1) \rightarrow 2$ & & $nx \rightarrow ny$ \\ |
| | 443 | \cline{1-3} |
| | 444 | & $(nx + 1) \rightarrow -1$ & $(nx + 1) \rightarrow -1$ & & \\ |
| | 445 | $\psi = w$ & $ nx \rightarrow 0$ & $ nx \rightarrow 0$ & & \\ |
| | 446 | & $(nx - 1) \rightarrow 1$ & $(nx - 1) \rightarrow 1$ & & \\ |
| | 447 | \hline |
| | 448 | \end{tabular} |
| | 449 | }}} |
| | 450 | |
| | 451 | |
| | 452 | ==== Constant Phase velocity ==== |
| | 453 | |
| | 454 | Setting ''c'',,ψ,, = ''c'',,max,, leads to a more simple radiation boundary condition (here e.g. for a left-right flow along positive ''x''-direction): |
| | 455 | {{{ |
| | 456 | #!Latex |
| | 457 | \begin{align*} |
| | 458 | \psi^{t + \Delta t}(k,j,nx+1) = \psi^{t}(k,j,nx) \; , |
| | 459 | \end{align*} |
| | 460 | }}} |
| | 461 | with ''ψ = {u,v,w}''. This formulation of the radiation boundary condtions saves computational time compared to the formulation of a variable Phase velocity. Although, [#orlanski1976 Orlanski (1976)] suggested that this approach of radiation boundary condition leads to reflection for waves smaller than ''c'',,max,, which may occur in complex geophysical flows, our simulations of a convective boundary layer with background wind have been stable so far. |
| | 462 | |
| | 463 | |
| | 464 | === Mass flux correction === |
| | 465 | |
| | 466 | PALM offers the possibility of a mass flux correction at the outflow (e.g. [#tian2004 Tian, 2004]). If parameter [../../app/inipar/#conserve_volume_flow conserve_volume_flow] is set true, the mass flux at the inflow and outflow is calculated by: |
| | 467 | {{{ |
| | 468 | #!Latex |
| | 469 | \begin{align*} |
| | 470 | \dot{m} = \sum_{k=1}^{nz-1} \Delta z(k) \sum_{l=0}^{nx_i} \psi(l,k) \Delta x_i \; . |
| | 471 | \end{align*} |
| | 472 | }}} |
| | 473 | where ''Δx,,i,,'' and ''ψ'' is equal to ''Δy'' (''Δx'') and ''u'' (''v'') in case of [../../app/inipar/#bc_lr bc_lr] ([../../app/inipar/#bc_ns bc_ns]). The correction factor for the outflow velocity, which is necessary due to different mass fluxes at inflow and outflow, can be calculated by |
| | 474 | {{{ |
| | 475 | #!Latex |
| | 476 | \begin{align*} |
| | 477 | \psi_{corr} = \dfrac{\dot{m}_{inflow} - \dot{m}_{outflow}}{A} \; , |
| | 478 | \end{align*} |
| | 479 | }}} |
| | 480 | where A is the area of the boundary |
| | 481 | {{{ |
| | 482 | #!Latex |
| | 483 | \begin{align*} |
| | 484 | A = \sum_1^{nz-1} \Delta z \sum_0^{nx_i} \Delta x_i \; . |
| | 485 | \end{align*} |
| | 486 | }}} |
| | 487 | The streamwise velocity at the outflow is corrected by adding ''ψ'',,corr,, at each grid point of the outflow between bottom and top (''k=1:nz''), in order to guarantee that the mass leaving the domain exactly balances the one entering it. |
| | 488 | |
| | 489 | |
| | 490 | == References == |
