1%------------------------------------------------------------------------------- 2 3% This file is part of Code_Saturne, a general-purpose CFD tool. 4% 5% Copyright (C) 1998-2021 EDF S.A. 6% 7% This program is free software; you can redistribute it and/or modify it under 8% the terms of the GNU General Public License as published by the Free Software 9% Foundation; either version 2 of the License, or (at your option) any later 10% version. 11% 12% This program is distributed in the hope that it will be useful, but WITHOUT 13% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15% details. 16% 17% You should have received a copy of the GNU General Public License along with 18% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19% Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21%------------------------------------------------------------------------------- 22 23\nopagebreak 24%================================== 25%================================== 26\section{Introduction} 27%================================== 28%================================== 29 30This document is a practical user guide for \CS version \verscs. 31It is the result of the joint effort of 32all the members in the development team. 33 34This document provides practical information for the usage of \CS. 35For more details about the algorithms and their numerical implementation, 36please refer to the reports 37 \cite{ijvf}, and \cite{mechitoua98}, 38and to the theoretical documentation \cite{theory}. 39 40The latest updated version of this document is available on-line with the version of \CS 41and accessible through the command 42\texttt{code\_saturne info --guide theory}. 43 44This document presents some the necessary elements to run a calculation 45with \CS version \verscs. It then lists all the variables of the code 46which may be useful for more advanced users. 47The user subroutines of all the modules within the code are then documented. 48Eventually, for each keyword and user-modifiable parameter in the code, 49their definition, allowed values, default values and conditions for use are given. 50These keywords and parameters are grouped under headings 51based on their function. An alphabetical index is also given at the end of 52the document for easier reference. 53 54In addition to the present user guide, a complete 55\texttt{Doxygen} documentation is available with \CS. It can provide 56information about the implementation such as details on variables used 57throughout the solver and the user subroutines. It also provides an easily 58explorable set of user subroutine examples and Fortran-C naming references for 59quantities linked to the mesh or the physical fields. 60 61The user documentation is in the process of migration from this pdf documentation 62to the Doxygen documentation, so the user should first lok there. 63One can access the \texttt{Doxygen} main page through \doxygenfile{index.html}{this link} or from a terminal by typing the following command: 64\texttt{code\_saturne info --guide theory}. 65 66%================================== 67\section{Basic modelling setup} 68%================================== 69 70%================================== 71\subsection{Manage boundary conditions} 72%================================== 73\label{sec:prg_boundaryconditions} 74The boundary conditions can be specified in the Graphical User Interface (GUI) under the heading ``Boundary conditions'' 75or in the user subroutine \texttt{cs\_user\_boundary\_conditions} called every time step. 76With the GUI, each region and the type of boundary condition associated to it are defined in 77\figurename~\ref{fig:gui_bc_regions}. Then, the parameters of the boundary condition are specified 78 in \figurename~\ref{fig:gui_bc_parameters}. The colors of the boundary faces may be read directly from 79a ``preprocessor.log'' file created by the Preprocessor. This file can be generated directly by the interface 80under the heading ``Definition of boundary regions $\rightarrow$ Add from Preprocessor log $\rightarrow$ import groups and references from Preprocessor log'', see \figurename~\ref{fig:gui_bc_regions}. 81% 82\begin{figure}[!ht] 83\begin{center} 84\includegraphics[width=0.9\textwidth]{gui_bc_regions} 85\caption{Definition of the boundary conditions} 86\label{fig:gui_bc_regions} 87\end{center} 88\end{figure} 89% 90\begin{figure}[!ht] 91\begin{center} 92\includegraphics[width=13cm]{gui_bc_parameters} 93\caption{Parameters of the boundary conditions} 94\label{fig:gui_bc_parameters} 95\end{center} 96\end{figure} 97% 98\texttt{cs\_user\_boundary\_conditions} is the second compulsory subroutine for every calculation launched 99without interface (except in the case of specific physics where the 100corresponding boundary condition user subroutine must be used). 101 102When using the interface, only complex boundary conditions (input profiles, conditions varying in time, ...) 103need to be defined with \texttt{cs\_user\_boundary\_conditions}. 104In the case of a calculation launched without the 105interface, all the boundary conditions must appear in \texttt{cs\_user\_boundary\_conditions}. 106 107\texttt{cs\_user\_boundary\_conditions} is essentially constituted of loops on boundary 108face subsets. Several sequences 109of \verb+call getfbr+ \verb+('criterion', nlelt, lstelt)+ (cf. 110\S\ref{sec:fvm_selector}) allow selecting 111the boundary faces with respect to their group(s), their 112color(s) or geometric criteria. If needed, geometric and 113physical variables are also available to the user. These allow him 114to select the boundary faces using other criteria. 115 116For more details about the treatment of boundary conditions, the user 117may refer to the theoretical and computer documentation \cite{theory} of 118the subroutine \texttt{condli} (for wall conditions, see 119\texttt{clptur}) (to access this document on a workstation, use 120\mbox{\texttt{code\_saturne~info --guide theory}}). 121 122From the user point of view, the boundary conditions are fully 123defined by three arrays\footnote{Except with Lagrangian boundary condition}: 124\texttt{itypfb(nfabor)}\index{\texttt{itypfb}}, 125\texttt{icodcl(nfabor,nvar)}\index{\texttt{icodcl}} and 126\texttt{rcodcl(nfabor,nvar,3)}\index{\texttt{rcodcl}}. 127\begin{list}{-}{} 128\item \texttt{itypfb(ifac)} defines the type of the face \texttt{ifac} 129 (input, wall, ...). 130\item \texttt{icodcl(ifac,ivar)} defines the type of boundary 131 condition for the variable \texttt{ivar} on the face \texttt{ifac} 132 (Dirichlet, flux ...). 133\item \texttt{rcodcl(ifac,ivar,.)} gives the numerical values associated with the 134 type of boundary condition (value of the Dirichlet condition, of the flux ...). 135\end{list} 136 137In the case of standard boundary conditions (see 138\S\ref{sec:prg_clstandard}), it is sufficient to complete \texttt{itypfb(ifac)} and 139parts of the array \texttt{rcodcl}; the array \texttt{icodcl} and most of \texttt{rcodcl} are filled automatically. For non-standard boundary 140conditions (see \S\ref{sec:prg_clnonstandard}), the arrays \texttt{icodcl} and 141\texttt{rcodcl} must be fully completed. 142 143%================================== 144\subsubsection{Coding of standard boundary conditions} 145%================================== 146\label{sec:prg_clstandard}% 147The standard keywords used by the indicator \texttt{itypfb} are: 148\texttt{ientre\index{ientre}}, \texttt{iparoi\index{iparoi}}, 149\texttt{iparug\index{iparug}}, \texttt{isymet\index{isymet}}, 150\texttt{isolib\index{isolib}}, \texttt{ifrent\index{ifrent}}, \texttt{ifresf\index{ifresf}}, 151\texttt{i\_convective\_inlet\index{i\_convective\_inlet}} and \texttt{iindef\index{iindef}}. 152 153\begin{list}{$\bullet$}{} 154\item If \texttt{itypfb=ientre}: inlet face. 155 156\begin{list}{$\rightarrow$}{} 157\item Zero-flux condition for pressure and Dirichlet condition for all 158 other variables. The value of the Dirichlet condition must be given in 159 \texttt{rcodcl(ifac,ivar,1)} for every value of \texttt{ivar}, except for 160 \texttt{ivar=ipr}. The other values of \texttt{rcodcl} and 161 \texttt{icodcl} are filled automatically. 162\end{list} 163 164\item If \texttt{itypfb=iparoi}: smooth solid wall face, impermeable and with friction. 165 166\begin{list}{$\rightarrow$}{} 167\item the eventual sliding wall velocity of the face is 168 found in \texttt{rcodcl(ifac,ivar,1)} (\texttt{ivar} being 169 \texttt{iu}, \texttt{iv} or \texttt{iw}). The initial 170 values of \texttt{rcodcl(ifac,ivar,1)} are zero for 171 the three velocity components (and therefore are to be specified 172 only if the velocity is not equal to zero). \\ 173{\em WARNING: the wall sliding velocity must belong to the boundary face 174 plane. For safety, the code only uses the projection of this 175 velocity on the face. As a consequence, if the velocity specified 176 by the user does not belong to the face plane, the wall sliding velocity really 177 taken into account will be different.} 178 179\item For scalars, two kinds of boundary conditions can be 180 defined: 181\begin{list}{$\rightsquigarrow$}{} 182\item Imposed value at the wall. The user must write\\ 183\hspace*{1cm}\texttt{icodcl(ifac,ivar)}=5\\ 184\hspace*{1cm}\texttt{rcodcl(ifac,ivar,1)}=imposed value\\ 185\item Imposed flux at the wall. The user must write\\ 186\hspace*{1cm}\texttt{icodcl(ifac,ivar)}=3\\ 187\hspace*{1cm}\texttt{rcodcl(ifac,ivar,3)}=imposed flux value (depending on the 188variable, the user may refer to the case \texttt{icodcl}=3 of \S~\ref{sec:prg_clnonstandard} for the flux definition). 189\item If the user does not fill these arrays, the default condition 190 is zero flux. 191\end{list} 192\end{list} 193 194\item If \texttt{itypfb=iparug}: rough solid wall face, impermeable and with friction. 195 196\begin{list}{$\rightarrow$}{} 197\item the eventual moving velocity of the wall tangent to the face is 198 given by \texttt{rcodcl(ifac,ivar,1)} (\texttt{ivar} being 199 \texttt{iu}, \texttt{iv} or \texttt{iw}). The initial 200 value of \texttt{rcodcl(ifac,ivar,1)} is zero for 201 the three velocity components (and therefore must be specified 202 only in the case of the existence of a slipping velocity). \\ 203{\em WARNING: the wall moving velocity must be in the boundary face 204 plane. By security, the code uses only the projection of this 205 velocity on the face. As a consequence, if the velocity specified 206 by the user is not in the face plane, the wall moving velocity really 207 taken into account will be different.} 208\item The dynamic roughness must be specified in \texttt{rcodcl(ifac,iu,3)}. 209 The values of \texttt{rcodcl(ifac,iv,3)} stores the thermal and scalar roughness. 210 The values of \texttt{rcodcl(ifac,iw,3)} is not used. 211\item For scalars, two kinds of boundary conditions can be defined: 212\begin{list}{$\rightsquigarrow$}{} 213\item Imposed value at the wall. The user must write\\ 214\hspace*{1cm}\texttt{icodcl(ifac,ivar)}=6\\ 215\hspace*{1cm}\texttt{rcodcl(ifac,ivar,1)}=imposed value\\ 216\item Imposed flux at the wall. The user must write\\ 217\hspace*{1cm}\texttt{icodcl(ifac,ivar)}=3\\ 218\hspace*{1cm}\texttt{rcodcl(ifac,ivar,3)}= imposed flux value (definition 219 of the flux condition according to the variable, the user can refer to the 220 case \texttt{icodcl}=3 of the paragraph \ref{sec:prg_clnonstandard}). 221\item If the user does not complete these arrays, the default condition 222 is zero flux. 223\end{list} 224\end{list} 225\item If \texttt{itypfb=isymet}: symmetry face (or wall without friction). 226\begin{list}{$\rightarrow$}{} 227\item Nothing to be writen in \texttt{icodcl} and \texttt{rcodcl}. 228\end{list} 229 230\item If \texttt{itypfb=isolib}: free outlet face (or more precisely free 231 inlet/outlet with forced pressure) 232\begin{list}{$\rightarrow$}{} 233\item The pressure is always treated with a Dirichlet condition, calculated 234 with the constraint $\displaystyle \frac{\partial }{\partial n}\left(\frac{ \partial P}{\partial \tau}\right)=0$. 235 The pressure is set to $P_0$ at the first \texttt{isolib} face met. 236 The pressure calibration is always done on a single face, even if there are 237 several outlets. 238\item If the mass flow is coming in, the velocity is set to zero 239 and a Dirichlet condition for the scalars and the turbulent quantities is used 240 (or zero-flux condition if no Dirichlet value has been specified). 241\item If the mass flow is going out, zero-flux condition are set for the velocity, 242 the scalars and the turbulent quantities. 243\item Nothing is written in \texttt{icodcl} or \texttt{rcodcl} for the pressure or 244 the velocity. An optional Dirichlet condition can be specified for the scalars 245 and turbulent quantities. 246\end{list} 247 248\item If \texttt{itypfb=ifrent}: free outlet, free inlet (based on Bernoulli relationship) face. 249 250\begin{list}{$\rightarrow$}{} 251\item if outlet, the equivalent to standard outlet. 252 In case of ingoing flux, the Benoulli relationship which links pressure and velocity is used (see the thory guide for more information). An additional head loss modelling what is going on outward of the domain can be added by the user. 253\end{list} 254 255\item If \texttt{itypfb=ifresf}: free-surface boundary condition. 256 257\item If \texttt{itypfb=i\_convective\_inlet}: inlet with zero diffusive flux for all transported variables (species and velocity). This allows to exactly impose the ingoing flux. 258 259\item If \texttt{itypfb=iindef}: undefined type face (non-standard case). 260\begin{list}{$\rightarrow$}{} 261\item Coding is done in a non-standard way by filling both arrays \texttt{rcodcl} and 262 \texttt{icodcl} (see \S~\ref{sec:prg_clnonstandard}). 263\end{list} 264\end{list} 265 266\minititre{Notes} 267 268$\bullet\ $ Whatever is the value of the indicator \texttt{itypfb(ifac)}, if 269the array \texttt{icodcl(ifac,ivar)} is modified by the user ({\em i.e.} filled 270with a non-zero value), the code will not use the default 271conditions for the variable \texttt{ivar} at the face \texttt{ifac}. It will 272take into account only the values of \texttt{icodcl} and \texttt{rcodcl} provided by the 273user (these arrays must then be fully completed, like in the non-standard case). \\ 274For instance, for a normal symmetry face where scalar 1 is associated with a 275Dirichlet condition equal to 23.8 (with an infinite exchange 276coefficient):\\ 277\hspace*{2cm}\texttt{itypfb(ifac)=isymet}\\ 278\hspace*{2cm}\texttt{icodcl(ifac,isca(1))=1}\\ 279\hspace*{2cm}\texttt{rcodcl(ifac,isca(1),1)=23.8}\\ 280(\texttt{rcodcl(ifac,isca(1),2)=rinfin} is the default value, therefore it is 281not necessary to specify a value)\\ 282The boundary conditions for the other variables are automatically 283defined. 284 285\noindent 286$\bullet\ $The user can define new types of boundary faces. He only must 287choose a value $N$ and to fully specify the boundary conditions (see 288\S\ref{sec:prg_clnonstandard}). He must specify 289\texttt{itypfb(ifac)}=$N$ where $N$ range is 1 to 290\texttt{ntypmx\index{ntypmx}} (maximum number of boundary face types), and of 291course different from the values \texttt{ientre}, \texttt{iparoi}, 292\texttt{iparug}, \texttt{isymet}, \texttt{isolib} and \texttt{iindef} (the values 293of these variables are given in the \texttt{paramx} module). This allows to 294easily isolate some boundary faces, in order for instance to calculate balances. 295 296\minititre{Boundary condition types} 297 298The \textbf{gradient} boundary conditions in \CS boil down to determine a 299value for the current variable $\varia$ at the boundary faces $\fib$, that is to 300say $\varia_\fib$, value expressed as a function of $\varia_{\centip}$, value 301of $\varia$ in $\centip$, projection of the center of the adjacent cell on the 302straight line perpendicular to the boundary face and crossing its center: 303\begin{equation} 304\varia_\fib=A_{\fib}^g +B_{\fib}^g \varia_{\centip}. 305\end{equation} 306 307For a face \texttt{ifac}, the pair of coefficients $A_{\fib}^g , \, B_{\fib}^g$ is 308may be accessed using the \texttt{field\_get\_coefa\_s} and 309\texttt{field\_get\_coefb\_s} functions, replacing \texttt{s} with \texttt{v} 310for a vector. 311 312The \textbf{flux} boundary conditions in \CS boil down to determine the 313value of the diffusive flux of the 314current variable $\varia$ at the boundary faces $\fib$, that is to say 315 $D_{\ib} \left(K_\fib, \, \varia \right)$, 316value expressed as a function of $\varia_{\centip}$, value of $\varia$ in $\centip$, 317projection of the center of the adjacent cell on the straight line 318perpendicular to the boundary face and crossing its center: 319\begin{equation} 320D_{\ib} \left(K_\fib, \, \varia \right) = A_{\fib}^f +B_{\fib}^f \varia_{\centip}. 321\end{equation} 322 323For a face \texttt{ifac}, the pair of coefficients $A_{\fib}^f , \, B_{\fib}^f$ 324may be accessed using the \texttt{field\_get\_coefaf\_s} and 325\texttt{field\_get\_coefbf\_s} functions, replacing \texttt{s} with \texttt{v} 326for a vector. 327 328The \textbf{divergence} boundary conditions in \CS boil down to determine a value for the 329current variable $\varia$ (mainly the Reynolds stress components, the divergence $\divv \left(\tens{R} \right)$ used in the calculation of the momentum equation) at the boundary 330faces $\fib$, that is to say $\varia_\fib$, 331value expressed as a function of $\varia_{\centip}$, value of $\varia$ in $\centip$, 332projection of the center of the adjacent cell on the straight line 333perpendicular to the boundary face and crossing its center: 334\begin{equation} 335\varia_\fib=A_{\fib}^d +B_{\fib}^d \varia_{\centip}. 336\end{equation} 337 338For a face \texttt{ifac}, the pair of coefficients $A_{\fib}^d , \, B_{\fib}^d$ 339may be accessed using the \texttt{field\_get\_coefad\_s} and 340\texttt{field\_get\_coefbd\_s} functions, replacing \texttt{s} with \texttt{v} 341for a vector. 342 343%================================== 344\subsubsection{Coding of non-standard boundary conditions} 345%================================== 346\label{sec:prg_clnonstandard}% 347Ifa face does not correspond to a standard type, the user 348must completely fill the arrays \texttt{itypfb}, \texttt{icodcl} and 349\texttt{rcodcl}. \texttt{itypfb(ifac)} is then equal to \texttt{iindef} 350or another value defined by the user (see note at the end of 351\S~\ref{sec:prg_clstandard}). The arrays \texttt{icodcl} and \texttt{rcodcl} 352must be filled as follows: 353 354\begin{list}{$\bullet$}{} 355\item If \texttt{icodcl(ifac,ivar)}=1: Dirichlet condition at the face 356 \texttt{ifac} for the variable \texttt{ivar}. 357 358\begin{list}{$\rightarrow$}{} 359\item \texttt{rcodcl(ifac,ivar,1)} is the value of the variable \texttt{ivar} 360 at the face \texttt{ifac}. 361 362\item \texttt{rcodcl(ifac,ivar,2)} is the value of the exchange coefficient 363 between the outside and the fluid for the variable \texttt{ivar}. An 364 infinite value (\texttt{rcodcl(ifac,ivar,2)=rinfin}) indicates an 365 ideal transfer between the outside and the fluid (default case). 366 367\item \texttt{rcodcl(ifac,ivar,3)} is not used. 368 369\item \texttt{rcodcl(ifac,ivar,1)} has the units of the variable 370 \texttt{ivar}, {\em i.e.}: 371\begin{list}{$\rightsquigarrow$}{} 372\item $m/s$ for the velocity 373 374\item $m^2/s^2$ for the Reynolds stress 375 376\item $m^2/s^3$ for the dissipation 377 378\item $Pa$ for the pressure 379 380\item \degresC\ for the temperature 381 382\item $J.kg^{-1}$ for the enthalpy 383 384\item \degresC$^2$ for temperature fluctuations 385 386\item $J^2.kg^{-2}$ for enthalpy fluctuations 387\end{list} 388 389\item \texttt{rcodcl(ifac,ivar,2)} has the following units (defined in such way 390 that when multiplying the exchange coefficient by the variable, the 391 given flux has the same units as the flux defined below when 392 \texttt{icodcl=3}): 393 394\begin{list}{$\rightsquigarrow$}{} 395\item $kg.m^{-2}.s^{-1}$ for the velocity 396 397\item $kg.m^{-2}.s^{-1}$ for the Reynolds stress 398 399\item $s.m^{-1}$ for the pressure 400 401\item $W.m^{-2}.\mbox{\degresC}^{-1}$ for the temperature 402 403\item $kg.m^{-2}.s^{-1}$ for the enthalpy 404\end{list} 405 406\end{list} 407 408\item If \texttt{icodcl(ifac,ivar)=2}: radiative outlet at the face \texttt{ifac} 409 for the variable \texttt{ivar}. It reads $ \dfrac{\partial \varia }{\partial t} + C \dfrac{\partial \varia}{\partial n} = 0 $, where $C$ is a to be defined celerity of radiation. 410 411\begin{list}{$\rightarrow$}{} 412\item \texttt{rcodcl(ifac,ivar,3)} is not used. 413 414\item \texttt{rcodcl(ifac,ivar,1)} is the flux value of \texttt{ivar} at the cell center $\centip$, 415 projection of the center of the adjacent cell on the straight line 416 perpendicular to the boundary face and crossing its center, 417 at the previous time step. 418 It corresponds to: 419\item \texttt{rcodcl(ifac,ivar,2)} is CFL number based on the parameter $C$, 420 the distance to the boundary $\centip \centf$ and the time step: 421 $CFL = \dfrac{C dt }{\centip \centf}$, 422 423\end{list} 424 425\item If \texttt{icodcl(ifac,ivar)=3}: flux condition at the face \texttt{ifac} 426 for the variable \texttt{ivar}. 427 428\begin{list}{$\rightarrow$}{} 429\item \texttt{rcodcl(ifac,ivar,1)} and \texttt{rcodcl(ifac,ivar,2)} are not used. 430 431\item \texttt{rcodcl(ifac,ivar,3)} is the flux value of \texttt{ivar} at the 432 wall. This flux is negative if it is a source for the fluid. It corresponds to: 433\begin{list}{$\rightsquigarrow$}{} 434\item 435$\displaystyle -(\lambda_T+C_p\frac{\mu_t}{\sigma_T})\grad T\cdot\vect{n}$ for a temperature (in $W/m^2$) 436 437$\displaystyle -(\frac{\lambda_T}{C_p}+\frac{\mu_t}{\sigma_h})\grad h\cdot\vect{n}$ 438 for an enthalpy (in $W/m^2$). 439 440$\displaystyle -(\lambda_\varphi+\frac{\mu_t}{\sigma_\varphi})\grad\varphi\cdot\vect{n}$ in the case of another scalar $\varphi$ (in $kg.m^{-2}.s^{-1}.[\varphi]$, where $[\varphi]$ are the units of $\varphi$). 441 442\item $-\Delta t\ \grad P\cdot\vect{n}$ for the pressure (in $kg.m^{-2}.s^{-1}$). 443 444\item $-(\mu+\mu_t)\grad U_i\cdot\vect{n}$ for a velocity component (in $kg.m^{-1}.s^{-2}$). 445 446\item $-\mu\grad R_{ij}\cdot\vect{n}$ for a $R_{ij}$ tensor component (in $W/m^2$). 447\end{list} 448 449\end{list} 450 451\item If \texttt{icodcl(ifac,ivar)}=4: symmetry condition, for the symmetry 452 faces or wall faces without friction. This condition can only be 453 used for velocity components ($\vect{U}\cdot\vect{n}=0$) and 454 the $R_{ij}$ tensor components (for other variables, a zero-flux 455 condition type is usually used).\\ 456 457\item If \texttt{icodcl(ifac,ivar)}=5: friction condition, for wall faces 458 with friction. This condition can not be applied to the pressure. 459\begin{list}{$\rightsquigarrow$}{} 460\item For the velocity and (if necessary) the turbulent variables, the 461 values at the wall are calculated from theoretical profiles. In 462 the case of a sliding wall, the three components of the sliding 463 velocity are given by (\texttt{rcodcl(ifac,iu,1)}, 464 \texttt{rcodcl(ifac,iv,1)}, and \\\texttt{rcodcl(ifac,iw,1)}).\\ 465{\em WARNING: the wall sliding velocity must belong to the boundary face 466 plane. For safety, the code uses only the projection of this 467 velocity on the face. Therefore, if the velocity vector specified 468 by the user does not belong to the face plane, the wall sliding velocity really 469 taken into account will be different.} 470 471\item For other scalars, the condition \texttt{icodcl}=5 is similar to 472 \texttt{icodcl=1}, but with a wall exchange coefficient calculated from a 473 theoretical law. Therefore, the values of \\\texttt{rcodcl(ifac,ivar,1)} and 474 \texttt{rcodcl(ifac,ivar,2)} must be specified: see \cite{theory}. 475\end{list} 476 477\item If \texttt{icodcl(ifac,ivar)}=6: friction condition, for the rough-wall faces 478 with friction. This condition can not be used with the pressure. 479\begin{list}{$\rightsquigarrow$}{} 480\item For the velocity and (if necessary) the turbulent variables, the 481 values at the wall are calculated from theoretical profiles. In 482 the case of a sliding wall, the three components of the sliding 483 velocity are given by (\texttt{rcodcl(ifac,iu,1)}, 484 \texttt{rcodcl(ifac,iv,1)}, and \\\texttt{rcodcl(ifac,iw,1)}).\\ 485 {\em WARNING: the wall sliding velocity must belong to the boundary face 486 plane. For safety, the code uses only the projection of this 487 velocity on the face. Therefore, if the velocity vector specified 488 by the user does not belong to the face plane, the wall sliding velocity really 489 taken into account will be different.}\\ 490 The dynamic roughness height is given by \texttt{rcodcl(ifac,iu,3)} only. 491 492\item For the other scalars, the condition \texttt{icodcl}=6 is similar to 493 \texttt{icodcl}=1, but with a wall exchange coefficient calculated from a 494 theoretical law. The values of \texttt{rcodcl(ifac,ivar,1)} and 495 \texttt{rcodcl(ifac,ivar,2)} must therefore be specified: see \cite{theory}. 496 The thermal roughness height is then given by \texttt{rcodcl(ifac,ivar,3)}. 497\end{list} 498 499\item If \texttt{icodcl(ifac,ivar)}=9: free outlet condition for the 500 velocity. This condition is only applicable to velocity 501 components.\\ 502 If the mass flow at the face is negative, this condition is equivalent 503 to a zero-flux condition.\\ 504 If the mass flow at the face is positive, the velocity at the face is set to zero (but not the mass flow).\\ 505\texttt{rcodcl} is not used. 506 507\item If \texttt{icodcl(ifac,ivar)}=14: generalized symmetry boundary condition for vectors (Marangoni 508 effect for the velocity for instance). 509 This condition is only applicable to vectors and set a Dirichlet boundary condition on the normal 510 component and a Neumann condition on the tangential components.\\ 511 If the three components are \texttt{ivar1, ivar2, ivar3}, the required values are: 512 513\begin{list}{$\rightarrow$}{} 514 \item \texttt{rcodcl(ifac,ivar1,1)}: Dirichlet value in the $x$ direction. 515 \item \texttt{rcodcl(ifac,ivar2,1)}: Dirichlet value in the $y$ direction. 516 \item \texttt{rcodcl(ifac,ivar3,1)}: Dirichlet value in the $z$ direction. 517 \item \texttt{rcodcl(ifac,ivar1,3)}: flux value for the $x$ direction. 518 \item \texttt{rcodcl(ifac,ivar2,3)}: flux value for the $y$ direction. 519 \item \texttt{rcodcl(ifac,ivar3,3)}: flux value for the $z$ direction. 520\end{list} 521 Therefore, the code automatically computes the boundary condition to impose to the normal and to 522 the tangential components. 523 524\end{list} 525 526\minititre{Note} 527$\bullet\ $A standard \texttt{isolib} outlet face amounts to a Dirichlet 528condition (\texttt{icodcl}=1) for the pressure, a free outlet condition 529(\texttt{icodcl}=9) for the velocity and a Dirichlet condition 530(\texttt{icodcl}=1) if the user has specified a Dirichlet value or a zero-flux 531condition (\texttt{icodcl}=3) for the other variables.\\ 532 533%================================== 534\subsubsection{Checking of the boundary conditions} 535%================================== 536 537The code checks the main compatibilities between the boundary 538conditions. In particular, the following rules must be respected: \\ 539$\bullet\ $On each face, the boundary conditions of the three velocity components must belong to the same type. The same is true for the components of the $R_{ij}$ tensor.\\ 540$\bullet\ $If the boundary conditions for the velocity belong to the 541``sliding'' type (\texttt{icodcl}=4), the conditions for $R_{ij}$ must belong to 542the ``symmetry'' type (\texttt{icodcl}=4), and vice versa.\\ 543$\bullet\ $If the boundary conditions for the velocity belong to the 544``friction'' type (\texttt{icodcl}=5 or 6), the boundary conditions for the turbulent variables 545must belong to the ``friction'' type, too.\\ 546$\bullet\ $If the boundary condition of a scalar belongs to the 547``friction'' type, the boundary condition of the velocity must belong to 548the ``friction'' type, too. 549 550In case of mistakes, if the post-processing output is activated (which is the default setting), 551a special error output, similar to the mesh format, is produced in order to help 552correcting boundary condition definitions. 553 554%================================== 555\subsubsection{Sorting of the boundary faces} 556%================================== 557 558In the code, it may be necessary to have access to all the boundary 559faces of a given type. To ease this kind of search, an array made of 560sorted faces is automatically filled (and updated at each time step): 561\texttt{itrifb(nfabor)\index{itrifb}}.\\ 562\texttt{ifac=itrifb(i)} is the number of the i$^{\text{th}}$ face of type 5631.\\ 564\texttt{ifac=itrifb(i+n)} is the number of the i$^{\text{th}}$ face of type 5652, if there are $n$ faces of type 1.\\ 566... etc. 567 568Two auxiliary arrays of size \texttt{ntypmx} are also defined.\\ 569\texttt{idebty(ityp)\index{idebty}} is the index 570corresponding to the first 571face of type \texttt{ityp} in the array \texttt{itrifb}.\\ 572\texttt{ifinty(ityp)\index{ifinty}} is the index 573corresponding to the last 574face of type \texttt{ityp} in the array \texttt{itrifb}. 575 576Therefore, a value \texttt{ifac0} found between \texttt{idebty(ityp)} and 577\texttt{ifinty(ityp)} is associated to each face \texttt{ifac} of type 578\texttt{ityp=itypfb(ifac)}, so that \texttt{ifac=itrifb(ifac0)}. 579 580If there is no face of type \texttt{ityp}, the code set \\ 581\texttt{ifinty(ityp)=idebty(ityp)-1},\\ 582which enables to bypass, for all the missing \texttt{ityp}, the loops such as \\ 583\texttt{do ii=idebty(ityp),ifinty(ityp)}. 584 585The values of all these indicators are displayed at the beginning of the 586code execution log. 587 588%================================== 589\subsection{User source terms} 590%================================== 591\label{sec:prg_usersourceterms} 592 593Assume, for example, that the user source terms modify the equation of a 594variable $\varphi$ in the following way: 595\begin{displaymath} 596\rho\frac{\partial \varphi}{\partial t}+\ldots = \ldots + S_{impl}\times\varphi+S_{expl} 597\end{displaymath} 598The example is valid for a velocity component, for a turbulent variable ($k$, $\varepsilon$, $R_{ij}$, $\omega$, 599$\varphi$ or $\overline{f}$) and for a scalar (or for the average of the 600square of the fluctuations of a scalar), because the syntax of all the 601subroutines \texttt{ustsnv}, \texttt{cs\_user\_turbulence\_source\_terms} 602 and \texttt{ustssc} in the \texttt{cs\_user\_source\_terms} file is similar. 603 604In the finite volume formulation, the solved system is then modified as 605follows: 606\begin{displaymath} 607\left(\frac{\rho_i\Omega_i}{\Delta t_i}-\Omega_iS_{impl,i}\right) 608\left(\varphi_i^{(n+1)}-\varphi_i^{(n)}\right) 609+\ldots = \ldots + \Omega_iS_{impl,i}\varphi_i^{(n)} + \Omega_iS_{expl,i} 610\end{displaymath} 611The user needs therefore to provide the following values:\\ 612$\text{\texttt{crvimp}}_i=\Omega_iS_{impl,i}$\\ 613$\text{\texttt{crvexp}}_i=\Omega_iS_{expl,i}$ 614 615In practice, it is essential for the term 616$\displaystyle \left(\frac{\rho_i\Omega_i}{\Delta 617t_i}-\Omega_iS_{impl,i}\right)$ to be positive. To ensure this property, 618the equation really taken into account by the code is the following: 619\begin{displaymath} 620\left(\frac{\rho_i\Omega_i}{\Delta t_i}- 621\text{Min}(\Omega_iS_{impl,i};0)\right) 622\left(\varphi_i^{(n+1)}-\varphi_i^{(n)}\right) 623+\ldots = \ldots + \Omega_iS_{impl,i}\varphi_i^{(n)} + \Omega_iS_{expl,i} 624\end{displaymath} 625To make the ``implicitation'' effective, the source term decomposition 626between the implicit and explicit parts will be done by the user who must 627ensure that $\text{\texttt{crvimp}}_i=\Omega_iS_{impl,i}$ is always negative 628(otherwise the solved equation remains right, but there will not be 629``implicitation''). 630 631{\em WARNING: When the second-order in time is used along with the extrapolation of the 632source terms\footnote{indicator \texttt{isno2t} for the velocity, 633\texttt{isto2t} for the turbulence and \texttt{isso2t} for the scalars}, it is no longer possible to test the sign of $S_{impl,i}$, 634because of coherence reasons (for more details, the user may refer to 635the theoretical and computer documentation \cite{theory} of the 636subroutine \texttt{preduv}). The user must therefore make sure it is 637always positive (or take the risk to affect the calculation stability).} 638 639\minititre{Particular case of a linearised source term} 640 641In some cases, the added source term is not linear, but the user may 642want to linearise it using a first-order Taylor development, in order to 643make it partially implicit.\\ 644Consider an equation of the type: 645\begin{displaymath} 646\rho\frac{\partial\varphi}{\partial t}=F(\varphi) 647\end{displaymath} 648To make it implicit using the following method: 649\begin{eqnarray*} 650\frac{\rho_i\Omega_i}{\Delta t}\left(\varphi_i^{(n+1)}-\varphi_i^{(n)}\right) & = & 651\Omega_i\left[F(\varphi_i^{(n)})+\left(\varphi_i^{(n+1)}-\varphi_i^{(n)}\right) 652\frac{dF}{d\varphi}(\varphi_i^{(n)})\right]\\ 653& = & \Omega_i\frac{dF}{d\varphi}(\varphi_i^{(n)})\times\varphi_i^{(n+1)} 654+\Omega_i\left[F(\varphi_i^{(n)})-\frac{dF}{d\varphi}(\varphi_i^{(n)}) 655\times\varphi_i^{(n)}\right] 656\end{eqnarray*} 657 658The user must therefore specify:\\ 659$\displaystyle\text{\texttt{crvimp}}_i=\Omega_i\frac{dF}{d\varphi}(\varphi_i^{(n)})$\\ 660$\displaystyle\text{\texttt{crvexp}}_i= 661\Omega_i\left[F(\varphi_i^{(n)})-\frac{dF}{d\varphi}(\varphi_i^{(n)})\times\varphi_i^{(n)}\right]$ 662 663\underline{\em Example}:\\ 664If the equation is 665$\displaystyle \rho\frac{\partial\varphi}{\partial t}=-K\varphi^2$, 666the user must set:\\ 667$\text{\texttt{crvimp}}_i=-2K\Omega_i\varphi_i^{(n)}$\\ 668$\text{\texttt{crvexp}}_i=K\Omega_i[\varphi_i^{(n)}]^2$ 669 670%================================== 671\subsubsection{In Navier-Stokes} 672%================================== 673 674The source term in Navier-Stokes can be filled in thanks to the GUI or the 675\texttt{cs\_user\_source\_terms} user file. 676Without the GUI, the subroutine \texttt{ustsnv} is used to 677add user source terms to the Navier-Stokes equations (at each time step). 678 679\texttt{ustsnv} is called only once per time step; for each cell \texttt{iel}, 680the vector \texttt{crvexp(.,iel)} (explicit part) and the matrix \texttt{crvimp(.,.,iel)} 681(implicit part) must be filled in for the whole velocity vector. 682 683 684%================================== 685\subsubsection{For $k$ and $\varepsilon$} 686%================================== 687 688\noindent 689\textit{Subroutine called every time step, for the $k-\varepsilon$ and 690the v2f models.} 691 692The subroutine \texttt{cs\_user\_turbulence\_source\_terms} is used to add source terms to the transport equations 693related to the turbulent kinetics energy $k$ and to the turbulent 694dissipation $\varepsilon$. 695This subroutine is called every time step (the 696treatment of the two variables $k$ and $\varepsilon$ is made 697simultaneously). The user is expected to provide the arrays \texttt{crkimp} and 698\texttt{crkexp} for $k$, and \texttt{creimp} and \texttt{creexp} for 699$\varepsilon$. These arrays are similar to the arrays \texttt{crvimp} and 700\texttt{crvexp} given for the velocity in the user subroutine \texttt{ustsnv}. 701The way of making implicit the resulting source terms is the same as the one 702presented in \texttt{ustsnv}. For $\varphi$ and $\bar{f}$ 703in the v2f model, see \texttt{cs\_user\_turbulence\_source\_terms}, \S\ref{sec:prg_ustsv2}. 704 705%================================== 706\subsubsection{For $R_{ij}$ and $\varepsilon$} 707%================================== 708 709\noindent 710\textit{Subroutine called every time step, for the $R_{ij}-\varepsilon$ models.} 711 712The subroutine \texttt{cs\_user\_turbulence\_source\_terms} is used to add source terms to the transport equations 713related to the Reynolds stress variables $R_{ij}$ and to the turbulent 714dissipation $\varepsilon$. 715This subroutine is called 7 times every time step 716(once for each Reynolds stress component and once for the 717dissipation). The user must provide the arrays \texttt{crvimp} and 718\texttt{crvexp} for the field variable of index \texttt{f\_id} (referring successively to 719\texttt{ir11}, \texttt{ir22}, \texttt{ir33}, 720\texttt{ir12}, \texttt{ir13}, \texttt{ir23} and 721\texttt{iep}). These arrays are similar to the arrays \texttt{crvimp} 722and \texttt{crvexp} given for the velocity in the user subroutine 723\texttt{ustsnv}. The method for impliciting the resulting source terms is the 724same as that presented in \texttt{ustsnv}. 725 726%================================== 727\subsubsection{For $\varphi$ and $\overline{f}$} 728%================================== 729\label{sec:prg_ustsv2} 730 731\noindent 732\textit{Subroutine called every time step, for the v2f models.} 733 734The subroutine \texttt{cs\_user\_turbulence\_source\_terms} is used to add source terms to the transport equations 735related to the variables $\varphi$ and $\overline{f}$ of the v2f 736$\varphi$-model. This subroutine is called twice 737every time step (once for $\varphi$ and once for $\overline{f}$). 738The user is expected to provide the arrays \texttt{crvimp} and 739\texttt{crvexp} for \texttt{ivar} referring successively to \texttt{iphi} 740and \texttt{ifb}. Concerning $\varphi$, these arrays are similar to the arrays 741\texttt{crvimp} and \texttt{crvexp} given for the velocity in the user subroutine 742\texttt{ustsnv}. Concerning $\overline{f}$, the equation is slightly 743different: 744\begin{displaymath} 745L^2 div(\grad(\overline{f})) = \overline{f} + \ldots + S_{impl}\times\overline{f}+S_{expl} 746\end{displaymath} 747In the finite volume formulation, the solved system is written as: 748\begin{displaymath} 749\int_{\partial\Omega_i}\grad(\overline{f})^{(n+1)}dS=\frac{1}{L_i^2}\left( 750\Omega_i\overline{f}^{(n+1)}_i + \ldots + \Omega_iS_{impl,i}\overline{f}_i^{(n+1)} + 751\Omega_iS_{expl,i} \right) 752\end{displaymath} 753The user must then specify:\\ 754$\text{\texttt{crvimp}}_i=\Omega_iS_{impl,i}$\\ 755$\text{\texttt{crvexp}}_i=\Omega_iS_{expl,i}$ 756 757The way of making implicit the resulting source terms is the same as the 758one presented in \texttt{ustsnv}. 759 760%================================== 761\subsubsection{For $k$ and $\omega$} 762%================================== 763 764\noindent 765\textit{Subroutine called every time step, for the $k-\omega$ SST model.} 766 767The subroutine \texttt{cs\_user\_turbulence\_source\_terms} is used to add source terms to the transport equations 768related to the turbulent kinetics energy $k$ and to the specific 769dissipation rate $\omega$. This subroutine is 770called every time step (the treatment of the two 771variables $k$ and $\omega$ is made simultaneously). The user is expected 772to provide the arrays \texttt{crkimp} and \texttt{crkexp} for the variable 773$k$, and the arrays \texttt{crwimp} and \texttt{crwexp} for the variable $\omega$. 774These arrays are similar to the arrays \texttt{crvimp} and \texttt{crvexp} 775given for the velocity in the user subroutine \texttt{ustsnv}. The way of 776making implicit the resulting source terms is the same as the one presented in 777\texttt{ustsnv}. 778 779%================================== 780\subsubsection{For $\tilde{\nu}_t$} 781%================================== 782 783\noindent 784\textit{Subroutine called every time step, or the Spalart-Allmaras model.} 785 786The subroutine \texttt{cs\_user\_turbulence\_source\_terms} is used to add source terms to the transport equations 787related to the turbulent viscosity $\nu_t$ for the Spalart-Allmaras model. 788This subroutine is called every time step. The user is expected 789to provide the arrays \texttt{crkimp} and \texttt{crkexp} for the variable 790$\tilde{\nu}_t$. These arrays are similar to the arrays \texttt{crvimp} and \texttt{crvexp} 791given for the velocity in the user subroutine \texttt{ustsnv}. The way of 792making implicit the resulting source terms is the same as the one presented in 793\texttt{ustsnv}. 794 795%================================== 796\subsubsection{For user scalars} 797%================================== 798 799\noindent 800\textit{Subroutine called every time step.} 801 802The source terms in the transport equations related to the user scalars 803(passive or not, average of the square of the fluctuations of a scalar, ...) 804can be filled in thanks to the GUI or the \texttt{cs\_user\_source\_terms} user file. 805Without the GUI, the subroutine \texttt{ustssc} is used to add source terms to the 806transport equations related to the user scalars. In the same way as 807\texttt{ustsnv}, this subroutine is called every time step, once for 808each user scalar. The user must provide the arrays \texttt{crvimp} 809and \texttt{crvexp} related to each scalar. \texttt{cvimp} and \texttt{crvexp} 810must be set to 0 for the scalars on which it is not wished for the user source 811term to be applied (the arrays are initially set to 0 at each inlet in the subroutine). 812 813%================================== 814\section{Advanced modelling setup} 815%================================== 816 817%================================== 818\subsection{Use of a specific physics} 819%================================== 820\label{sec:prg_usppmo}% 821Specific physics such as dispersed phase, atmospheric flows, gas combustion, 822pulverised fuel combustion, electrical model and compressible model can be 823added by the user from the interface, or by using the subroutine \texttt{usppmo} of 824the \texttt{cs\_user\_parameters.f90} file (called only during the calculation initialisation). 825With the interface, when a specific physics is activated in \figurename~\ref{fig:5_GUI}, 826additional items or headings may appear (see for instance Sections \ref{sec:Ini-lag} 827and \ref{sec:Ini-coal}). 828 829When the interface is not used, \texttt{usppmo} is one of the three subroutines 830which must be completed by the user in order to use a specific physics module 831(only heavy fuel combustion is not available with the GUI). 832At the moment, \CS allows to use two ``pulverised coal'' modules 833(with Lagrangian coupling or not) and one ``pulverised heavy fuel'' module, 834two ``gas combustion'' modules, two ``electrical'' modules, 835a ``compressible'' module, and an ``atmospheric'' module. To activate one of 836these modules, the user must complete one (and only one) of the 837indicators \texttt{ippmod(i.....)\index{ippmod}} in the subroutine 838\texttt{usppmo}. By default, all the indicators \texttt{ippmod(i.....)} are 839initialised at -1, which means that no specific physics is activated. 840 841\begin{list}{$\bullet$}{} 842 \item Diffusion flame in the framework of ``3 points'' rapid complete 843 chemistry: indicator {\bf \tt ippmod(icod3p\index{icod3p})} 844 \begin{list}{$\rightarrow$}{} 845 \item \texttt{ippmod(icod3p)} = 0 adiabatic conditions 846 \item \texttt{ippmod(icod3p)} = 1 permeatic conditions (enthalpy 847 transport) 848 \item \texttt{ippmod(icod3p)} =-1 module not activated 849 \end{list} 850 \item Eddy Break Up pre-mixed flame: indicator {\bf \tt 851 ippmod(icoebu\index{icoebu})} 852 \begin{list}{$\rightarrow$}{} 853 \item \texttt{ippmod(icoebu\index{icoebu})} = 0 adiabatic 854 conditions at constant richness 855 \item \texttt{ippmod(icoebu)} = 1 permeatic conditions at 856 constant richness 857 \item \texttt{ippmod(icoebu)} = 2 adiabatic conditions at 858 variable richness 859 \item \texttt{ippmod(icoebu)} = 3 permeatic conditions at 860 variable richness 861 \item \texttt{ippmod(icoebu)} =-1 module not activated 862 \end{list} 863 \item Libby-Williams pre-mixed flame: indicator {\bf \tt ippmod(icolwc\index{icolwc})} 864 \begin{list}{$\rightarrow$}{} 865 \item \texttt{ippmod(icolwc)}=0 two peak model with adiabiatic conditions. 866 \item \texttt{ippmod(icolwc)}=1 two peak model with permeatic conditions. 867 \item \texttt{ippmod(icolwc)}=2 three peak model with adiabiatic conditions. 868 \item \texttt{ippmod(icolwc)}=3 three peak model with permeatic conditions. 869 \item \texttt{ippmod(icolwc)}=4 four peak model with adiabiatic conditions. 870 \item \texttt{ippmod(icolwc)}=5 four peak model with permeatic conditions. 871 \item \texttt{ippmod(icolwc)}=-1 module not activated. 872 \end{list} 873 \item Multi-coals and multi-classes pulverised coal combustion: 874 indicator {\bf \tt ippmod(iccoal\index{iccoal})} 875 The number of different coals must be less than or equal to 876 \texttt{ncharm\index{ncharm}} = 3. The number of particle size 877 classes \texttt{nclpch\index{nclpch}(icha)} for the coal 878 \texttt{icha}, must 879 be less than or equal to \texttt{ncpcmx\index{ncpcmx}} = 10. 880 \begin{list}{$\rightarrow$}{} 881 \item \texttt{ippmod(iccoal)} = 0 imbalance between the 882 temperature of the continuous and the solid phases 883 \item \texttt{ippmod(iccoal)} = 1 otherwise 884 \item \texttt{ippmod(iccoal)} =-1 module not activated 885 \end{list} 886 887 \item Multi-classes pulverised heavy fuel combustion: 888 indicator {\bf \tt ippmod(icfuel\index{icfuel})} 889 \begin{list}{$\rightarrow$}{} 890 \item \texttt{ippmod(icfuel)} = 0 module activated 891 \item \texttt{ippmod(icfuel)} =-1 module not activated 892 \end{list} 893 894 \item Lagrangian modelling of multi-coals and 895 multi-classes pulverised coal combustion: 896 indicator {\bf \tt ippmod(icpl3c\index{icpl3c})} 897 The number of different coals must be less than or equal to 898 \texttt{ncharm\index{ncharm}} = 3. The number of particle size 899 classes \texttt{nclpch\index{nclpch}(icha)} for the coal 900 \texttt{icha}, must be less than or equal to 901 \texttt{ncpcmx\index{ncpcmx}} = 10. 902 \begin{list}{$\rightarrow$}{} 903 \item \texttt{ippmod(icpl3c)} = 1 coupling with the Lagrangian 904 module, with transport of $H_2$ 905 \item \texttt{ippmod(icpl3c)} =-1 module not activated 906 \end{list} 907 \item Electric arcs module (Joule effect and Laplace forces): 908 indicator {\bf \tt ippmod(ielarc\index{ielarc})} 909 \begin{list}{$\rightarrow$}{} 910 \item \texttt{ippmod(ielarc)} = 1 determination of the magnetic field by 911 means of the Ampere's theorem (not available) 912 \item \texttt{ippmod(ielarc)} = 2 determination of the magnetic 913 field by means of the vector potential 914 \item \texttt{ippmod(ielarc)} =-1 module not activated 915 \end{list} 916 \item Joule effect module (Laplace forces not taken into account): 917 indicator {\bf \tt ippmod(ieljou\index{ieljou})} 918 \begin{list}{$\rightarrow$}{} 919 \item \texttt{ippmod(ieljou)} = 1 use of a real potential 920 \item \texttt{ippmod(ieljou)} = 2 use of a complex potential 921 \item \texttt{ippmod(ieljou)} = 3 use of real potential and specific boundary conditions for transformers. 922 \item \texttt{ippmod(ieljou)} = 4 use of complex potential and specific boundary conditions for transformers. 923 \item \texttt{ippmod(ieljou)} =-1 module not activated 924 \end{list} 925 \item Compressible module: indicator {\bf \tt 926 ippmod(icompf\index{icompf})} 927 \begin{list}{$\rightarrow$}{} 928 \item \texttt{ippmod(icompf)} = 0 module activated 929 \item \texttt{ippmod(icompf)} =-1 module not activated 930 \end{list} 931 \item Atmospheric flow module: indicator {\bf \tt 932 ippmod(iatmos\index{icompf})} 933 \begin{list}{$\rightarrow$}{} 934 \item \texttt{ippmod(iatmos)} =-1 module not activated 935 \item \texttt{ippmod(iatmos)} = 0 standard modelling 936 \item \texttt{ippmod(iatmos)} = 1 dry atmosphere 937 \item \texttt{ippmod(iatmos)} = 2 humid atmosphere 938 \end{list} 939% \item cooling towers module: indicator {\bf \tt 940% ippmod(iaeros\index{icompf})} 941% \begin{list}{$\rightarrow$}{} 942% \item \texttt{ippmod(iaeros\index)} =-1 module not activated 943% \item \texttt{ippmod(iaeros\index)} = 0 no model (NOT functional) 944% \item \texttt{ippmod(iaeros\index)} = 1 Poppe's model 945% \item \texttt{ippmod(iaeros\index)} = 2 Merkel's model 946% \end{list} 947\end{list} 948 949{\em WARNING: Only one specific physics module can be activated at the 950same time.} 951 952In the framework of the gas combustion modelling, the user may impose 953his own enthalpy-temperature tabulation (conversion law). He needs then 954to give the 955value zero to the indicator \texttt{indjon\index{indjon}} (the default value 956being 1). For more details, the user may refer to the following note 957(thermochemical files). 958 959\minititre{Note: the thermo-chemical files} 960The user must not forget to place in the directory DATA the 961thermochemical file \texttt{dp\_C3P}, \texttt{dp\_C3PSJ} or 962\texttt{dp\_ELE} (depending on the specific physics module he activated) 963Some example files are placed in the directory \texttt{DATA/REFERENCE} at the creation of the 964study case. Their content is described below. 965 966\begin{list}{$\bullet$}{} 967 \item Example of file for the gas combustion: 968 \begin{list}{$\rightarrow$}{} 969 \item if the enthalpy-temperature conversion data base 970 JANAF is used: \texttt{dp\_C3P} (see 971 array \ref{tab:dpC3P}). 972 973\begin{table}[htbp] 974\begin{center} 975\small{ 976\begin{tabular}{|c|c|c|c|} \hline 977 Lines &Examples of values & Variables & Observations \\ \hline 978 1 & 5 & \texttt{ngaze\index{ngaze}} & Number of current species \\ \hline 979 2 & 10 & \texttt{npo\index{npo}} & Number of points for the \\ 980 & & & enthalpy-temperature table \\ \hline 981 3 & 300. & \texttt{tmin\index{tmin}} & Lower temperature limit \\ 982 & & & for the table \\ \hline 983 4 & 3000. & \texttt{tmax\index{tmax}} & Upper temperature limi t \\ 984 & & & for the tabulation \\ \hline 985 5 & & & Empty line \\ \hline 986 6 & CH4 O2 CO2 H2O N2 & \texttt{nomcoe\index{nomcoe}}(\texttt{ngaze}) & List of the current species \\ \hline 987 7 &.35 .35 .35 .35 .35& \texttt{kabse\index{kabse}}(\texttt{ngaze}) & Absorption coefficient \\ 988 & & & of the current species \\ \hline 989 8 & 4 & \texttt{nato\index{nato}} & Number of elemental species \\ \hline 990 9 &.012 1 0 1 0 0& \texttt{wmolat\index{wmolat}}(\texttt{nato}), & Molar mass of the elemental \\ 991 10 &.001 4 0 0 2 0& & species (first column) \\ 992 11 &.016 0 2 2 1 0&\texttt{atgaze\index{atgaze}}(\texttt{ngaze},\texttt{nato})& Composition of the current species \\ 993 12 &.014 0 0 0 0 2& & as a function of the elemental species \\ 994 & & & (\texttt{ngaze} following columns) \\ \hline 995 13 & 3 & \texttt{ngazg\index{ngazg}} & Number of global species \\ 996 & & & Here, \texttt{ngazg} = 3 (Fuel, Oxidiser and Products) \\ \hline 997 14 & 1. 0. 0. 0. 0. & & Composition of the global species as a \\ 998 15 & 0. 1. 0. 0. 3.76 &\texttt{compog\index{compog}}(\texttt{ngaze},\texttt{ngazg})& function of the current species of line 6 \\ 999 16 & 0. 0. 1. 2. 7.52 & & In the order: Fuel (line 15), \\ 1000 & & & Oxidiser (line 16) and Product (line 17) \\ \hline 1001 17 & 1 & \texttt{nrgaz\index{nrgaz}} & Number of global reactions \\ 1002 & & & Here \texttt{nrgaz} = 1 (always equal to 1 \\ 1003 & & & in this version) \\ \hline 1004 18 & & \texttt{igfuel\index{igfuel}}(\texttt{nrgaz}), & Numbers of the global species concerned by \\ 1005 & 1 2 -1 -9.52 10.52& \texttt{igoxy\index{igoxy}}(\texttt{nrgaz}), & the stoichiometric ratio \\ 1006 & & & (first 2 integers) \\ 1007 & &\texttt{stoeg\index{stoeg}}(\texttt{ngazg},\texttt{nrgaz})& Stoichiometry in global species reaction. \\ 1008 & & & Negative for the reactants (here \\ 1009 & & & ``Fuel'' and ``Oxidiser'') and positive for \\ 1010 & & & the products (here ``Products'') \\ \hline 1011\end{tabular} 1012} 1013\caption{Example of file for the gas combustion when JANAF is used: \texttt{dp\_C3P}}\label{tab:dpC3P} 1014\end{center} 1015\end{table} 1016 1017 \item if the user provides his own enthalpy-temperature tabulation 1018 (there must be three chemical species and only 1019 one reaction): \texttt{dp\_C3PSJ} (see 1020 array \ref{tab:dpC3PSJ}). This file replaces \texttt{dp\_C3P}. 1021 1022\begin{table}[htbp] 1023\begin{center} 1024\small{ 1025\begin{tabular}{|c|c|c|c|} \hline 1026 Lines & Examples of values & Variables & Observations \\ \hline 1027 1 & 6 & \texttt{npo} & Number of tabulation points \\ \hline 1028 2 & 50. -0.32E+07 -0.22E+06 -0.13E+08& & \\ 1029 3 & 250. -0.68E+06 -0.44E+05 -0.13E+08&\texttt{th\index{th}}(\texttt{npo}), & Temperature(first column), \\ 1030 4 & 450. 0.21E+07 0.14E+06 -0.13E+08& \texttt{ehgazg\index{ehgazg}}(1,\texttt{npo}),& mass enthalpies of fuel, oxidiser \\ 1031 5 & 650. 0.50E+07 0.33E+06 -0.12E+08& \texttt{ehgazg}(2,\texttt{npo}), & and products (columns 2,3 and 4) \\ 1032 6 & 850. 0.80E+07 0.54E+06 -0.12E+08& \texttt{ehgazg}(3,\texttt{npo}) & from line 2 to line \texttt{npo}+1 \\ 1033 7 &1050. 0.11E+08 0.76E+06 -0.11E+08& & \\ \hline 1034 8 & .00219 .1387 .159 &\texttt{wmolg(1)\index{wmolg}}, & Molar masses of fuel, \\ 1035 & & \texttt{wmolg(2)},& oxidiser \\ 1036 & & \texttt{wmolg(3)} & and products \\ \hline 1037 9 & .11111 & \texttt{fs(1)\index{fs(1)}} & Mixing rate at the stoichiometry \\ 1038 & & & (relating to Fuel and Oxidiser) \\ \hline 1039 10 & 0.4 0.5 0.87 &\texttt{ckabsg\index{ckabsg}(1)}, & Absorption coefficients of the fuel, \\ 1040 & & \texttt{ckabsg(2)}, & oxidiser \\ 1041 & & \texttt{ckabsg(3)} & and products \\ \hline 1042 11 & 1. 2. & \texttt{xco2\index{xco2}}, \texttt{xh2o\index{xh2o}}& Molar coefficients of $CO_2$ \\ 1043 & & & and $H_2O$ in the products \\ 1044 & & & (using Modak radiation) \\ \hline 1045\end{tabular} 1046} 1047\caption{Example of file for the gas combustion when the user provides 1048 his own enthalpy-temperature table 1049 (there must be three species and only one 1050 reaction): \texttt{dp\_C3PSJ} (this file replaces 1051 \texttt{dp\_C3P})}\label{tab:dpC3PSJ} 1052\end{center} 1053\end{table} 1054 \end{list} 1055 1056 \item Example of file for the electric arcs: \texttt{dp\_ELE} (see 1057 array \ref{tab:dpELE}). 1058 1059\begin{table}[htbp] 1060\begin{center} 1061\small{ 1062\begin{tabular}{|c|l|c|c|} \hline 1063 Lines & Examples of values & Variables & Observations \\ \hline 1064 1 &\# Free format ASCII file ... & & Free comment \\ \hline 1065 2 &\# Comment lines ... & & Free comment \\ \hline 1066 3 &\# ... & & Free comment \\ \hline 1067 4 &\# Argon propoerties ... & & Free comment \\ \hline 1068 5 &\# ... & & Free comment \\ \hline 1069 6 &\# No of NGAZG and No ... & & Free comment \\ \hline 1070 7 &\# NGAZG NPO ... & & Free comment \\ \hline 1071 8 & 1 238 & \texttt{ngazg\index{ngazg}} & Number of species \\ 1072 & & \texttt{npo\index{npo}} & Number of given temperature points for \\ 1073 & & & the tabulated physical properties \\ 1074 & & & (\texttt{npo} $\leqslant$ \texttt{npot} set in \texttt{ppthch}) \\ 1075 & & & So there will be \texttt{ngazg} blocks of \texttt{npo} lines each \\ \hline 1076 9 &\# ... & & Free comment \\ \hline 1077 14 & 0 & \texttt{ixkabe\index{ixkabe}} & Radiation options for \texttt{xkabe\index{xkabe}} \\ \hline 1078 15 &\# ... & & Free comment \\ \hline 1079 16 &\# Propreties ... & & Free comment \\ \hline 1080 17 &\# ~~~~T~~~~~~~~~~~H ... & & Free comment \\ \hline 1081 18 &\# Temperature Enthalpy ... & & Free comment \\ \hline 1082 19 &\# ... & & Free comment \\ \hline 1083 20 &\# ~~~~~K~~~~~~~~~J/kg ... & & Free comment \\ \hline 1084 21 &\# ... & & Free comment \\ \hline 1085 22 & ~~~300.~~~~~~14000. ... & & In line tabulation of the physical properties \\ 1086 & & & as a function of the temperature in Kelvin \\ 1087 & & & for each of the \texttt{ngazg} species \\ 1088 & & \texttt{h} & Enthalpy in J/kg \\ 1089 & & \texttt{roel} & Density in kg/m3 \\ 1090 & & \texttt{cpel} & Specific heat in J/(kg K) \\ 1091 & & \texttt{sigel} & Electric conductivity in Ohm/m \\ 1092 & & \texttt{visel} & Dynamic viscosity in kg/(m s) \\ 1093 & & \texttt{xlabel} & Thermal conductivity in W/(m K) \\ 1094 & & \texttt{xkabel\index{xkabel}} & Absorption coefficient (radiation) \\ \hline 1095\end{tabular} 1096} 1097\caption{Example of file for the electric arcs module: 1098 \texttt{dp\_ELE}}\label{tab:dpELE} 1099\end{center} 1100\end{table} 1101 1102\end{list} 1103 1104\clearpage 1105 1106%================================== 1107\subsection[Pulverised coal and gas combustion module] 1108{Pulverised coal and gas combustion module (needs update)} 1109%================================== 1110%================================== 1111\subsubsubsection{Initialisation of the variables}\label{sec:Ini-coal} 1112%================================== 1113For coal combustion, it is possible to initialise the specific variables in the Graphical User Interface (GUI) or in the subroutine 1114\texttt{cs\_user\_initialization}. In the GUI, when a coal combustion physics is selected in the item ``Calculation features'' under the heading 1115``Thermophysical models'', an additional item appears: ``Pulverized coal combustion''. In this item the user can define coal types, their composition, the oxidant and reactions parameters, see \figurename~\ref{fig:Ini-coal1} to \figurename~\ref{fig:Ini-coal5}. 1116 1117\begin{figure}[!ht] 1118\begin{center} 1119\includegraphics[width=12cm]{gui_coal_classes} 1120\caption{Thermophysical models - Pulverized coal combustion, coal classes} 1121\label{fig:Ini-coal1} 1122\end{center} 1123\end{figure} 1124 1125\begin{figure}[!ht] 1126\begin{center} 1127\includegraphics[width=11cm]{gui_coal_composition} 1128\caption{Pulverized coal combustion, coal composition} 1129\label{fig:Ini-coal3} 1130\end{center} 1131\end{figure} 1132 1133\begin{figure}[!ht] 1134\begin{center} 1135\includegraphics[width=11cm]{gui_coal_reaction} 1136\caption{Pulverized coal combustion, reaction parameters} 1137\label{fig:Ini-coal4} 1138\end{center} 1139\end{figure} 1140 1141\begin{figure}[!ht] 1142\begin{center} 1143\includegraphics[width=11cm]{gui_coal_oxydant} 1144\caption{Pulverized coal combustion, oxydant} 1145\label{fig:Ini-coal5} 1146\end{center} 1147\end{figure} 1148 1149If the user deals with gas combustion or if he (or she) does not want to use the 1150GUI for coal combustion, the subroutine \texttt{cs\_user\_initialization} must be used (only during the calculation initialisation).\\ 1151In this section, ``specific physics'' will refer to gas combustion or 1152to pulverised coal combustion. 1153 1154These subroutines allow the user to initialise some variables specific 1155to the specific physics activated {\em via} \texttt{usppmo}. As usual, 1156the user may have access to several geometric variables to discriminate 1157between different initialisation zones if needed. 1158 1159It should be recalled again that the user can access the array of values of the 1160variables as described in the \doxygenfile{field.html}{the \texttt{doxygen} 1161documentation dedicated to the fields management}. In the following description, 1162only variables indices \texttt{ivar} are given, but field indices can be 1163retrieved easily by using \texttt{ivarfl(ivar)}. 1164 1165{\em WARNING: in the case of a specific physics modelling, all the 1166variables will be initialised here, even the potential user scalars: {\em 1167\texttt{cs\_user\_initialization}} is no longer used.} 1168 1169 1170\begin{list}{$\bullet$}{} 1171 \item in the case of the EBU pre-mixed flame module, the user can 1172 initialise in every cell \texttt{iel}: the mixing rate 1173 \texttt{isca(ifm)} in variable richness, the 1174 fresh gas mass fraction \\ 1175 \texttt{isca(iygfm\index{iygfm})} 1176 and the mixture enthalpy \texttt{isca(iscalt)} in 1177 permeatic conditions 1178 1179 \item in the case of the rapid complete chemistry diffusion flame 1180 module, the user can initialise in every cell \texttt{iel}: the 1181 mixing rate \texttt{isca(ifm\index{ifm})}, its variance 1182 \texttt{isca(ifp2m\index{ifp2m})} and the mixture mass 1183 enthalpy \texttt{isca(iscalt)} in permeatic conditions 1184 1185 \item in the case of the pulverised coal combustion module, the 1186 user can initialise in every cell \texttt{iel}: 1187 \begin{list}{$\rightarrow$}{} 1188 \item the transport variables related to the solid phase 1189 \begin{list}{}{} 1190 \item \texttt{isca(ixch\index{ixch}(icla))} 1191 the reactive coal mass fraction related to the 1192 class \texttt{icla} (\texttt{icla} from 1 to 1193 \texttt{nclacp} which is the total number of 1194 classes, {\em i.e.} for all the coal type) 1195 \item \texttt{isca(ixck(\index{ixck}icla))} 1196 the coke mass fraction related to the class 1197 \texttt{icla} 1198 \item \texttt{isca(inp\index{inp}(icla))} the 1199 number of particles related to class 1200 \texttt{icla} per kg of air-coal mixture 1201 \item \texttt{isca(ih2\index{ih2}(icla))} the 1202 mass enthalpy related to the class 1203 \texttt{icla} in permeatic conditions 1204 \end{list} 1205 \item \texttt{isca(iscalt)} the mixture enthalpy 1206 \item the transport variables related to the gas phase 1207 \begin{list}{}{} 1208 \item 1209 \texttt{isca(if1m\index{if1m}(icha))} the 1210 mean value of the tracer 1 representing 1211 the light volatile matters released by 1212 the coal \texttt{icha} 1213 \item 1214 \texttt{isca(if2m\index{if2m}(icha))} the 1215 mean value of the tracer 2 representing 1216 the heavy volatile matters released by 1217 the coal \texttt{icha} 1218 \item \texttt{isca(if3m\index{if3m})} 1219 the mean value of the tracer 3 1220 representing the carbon released 1221 as CO during coke burnout 1222 \item \texttt{isca(if4p2m\index{if4p2m})} the 1223 variance associated with the tracer 4 1224 representing the air (the mean value of this 1225 tracer is not transported, it can be deduced 1226 directly from the three others) 1227 \item \texttt{isca(ifp3m\index{ifp3m})} the 1228 variance associated with the tracer 3 1229 \end{list} 1230 \end{list} 1231\end{list} 1232 1233%================================== 1234\subsubsection{Boundary conditions}\label{sec:coal-cl} 1235%================================== 1236In this section, ``specific physics'' refers to gas combustion or 1237to pulverised coal combustion.\\ 1238For coal combustion, it is possible to manage the boundary conditions in the Graphical User Interface (GUI). When the coal combustion physics is selected in the heading ``Thermophysical models'', specific boundary conditions are activated for inlets, see \figurename~\ref{fig:cond_lim-coal}. The user fills for each type of coal previously defined (see \S~\ref{sec:Ini-coal}) the initial temperature and initial composition of the inlet flow, as well as the mass flow rate. 1239 1240\begin{figure}[!ht] 1241\begin{center} 1242\includegraphics[width=13cm]{gui_coal_bc} 1243\caption{Boundary conditions for the combustion of coal} 1244\label{fig:cond_lim-coal} 1245\end{center} 1246\end{figure} 1247 1248For gas combustion or if the GUI is not used for coal combustion, the use of 1249\texttt{cs\_user\_boundary\_conditions} (called at every time step) is as 1250mandatory as \texttt{cs\_user\_parameters.f90} and \texttt{usppmo} to run a calculation involving specific physics. The way of using them is the same as using 1251 in the framework of standard calculations, that is, run several loops on the boundary faces lists (cf. \S\ref{sec:fvm_selector}) 1252marked out by their colors, groups, or geometrical criterion, where 1253the type of face, the type of boundary condition for each variable and 1254eventually the value of each variable are defined. 1255 1256{\em WARNING: In the case of a specific physics modelling, all the 1257boundary conditions for every variable must be defined here, even for 1258the eventual user scalars: {\em \texttt{cs\_user\_boundary\_conditions}} is not used at all.}\\ 1259 1260In the case of a specific physics modelling, a zone number \texttt{izone} 1261\footnote{\texttt{izone} must be less than the maximum number of boundary 1262zone allowable by the code, \texttt{nozppm}. This is fixed at 2000 in 1263 \texttt{pppvar};not to be modified} (for 1264instance the color \texttt{icoul}) is associated with every boundary face, in 1265order to gather together all the boundary faces of the same type. In 1266comparison to \texttt{cs\_user\_boundary\_conditions}, the main change from the user point of 1267view concerns the faces whose boundary conditions belong to the type 1268\texttt{itypfb=ientre\index{ientre}}: 1269 1270\begin{list}{$\bullet$}{} 1271 \item for the EBU pre-mixed flame module: 1272 \begin{list}{$\rightarrow$}{} 1273 \item the user can choose between the ``burned gas 1274 inlet'' type (marked out by the burned gas indicator 1275 \texttt{ientgb\index{ientgb}(izone\index{izone})}=1) and the 1276 ``fresh gas inlet'' type (marked out by 1277 the fresh gas indicator 1278 \texttt{ientgf\index{ientgf}(izone)}=1) 1279 \item for each inlet type (fresh or burned 1280 gas), a mass flow or a velocity must be imposed: 1281 1282 \begin{list}{-}{} 1283 \item to impose the mass flow, 1284 \begin{list}{-}{} 1285 \item the user gives to 1286 the indicator 1287 \texttt{iqimp\index{iqimp}(izone)} 1288 the value 1, 1289 \item the 1290 mass flow value is set in 1291 \texttt{qimp\index{qimp}(izone)} 1292 (positive value, in $kgs^{-1}$) 1293 \item finally he imposes the 1294 velocity vector direction 1295 by giving the components of 1296 a direction vector in 1297 \texttt{rcodcl\index{rcodcl}(ifac,iu\index{iu})}, \texttt{rcodcl(ifac,iv\index{iv})} and \texttt{rcodcl(ifac,iw\index{iw})} 1298 \end{list} 1299 1300{\em WARNING: 1301\begin{list}{-}{} 1302\item the variable \texttt{qimp(izone)} refers to the mass flow across the whole 1303 zone \texttt{izone} and not across a boundary face (specifically for the axi-symmetric calculations, the inlet surface of the mesh must be broken up) 1304\item the variable \texttt{qimp(izone)} deals with the inflow across the area \texttt{izoz} and only across this zone; it is recommended to pay attention to the boundary conditions. 1305\item the velocity direction vector is neither necessarily normed, nor 1306 necessarily incoming. 1307\end{list}} 1308 1309 \item to impose a velocity, the user 1310 must give to the indicator 1311 \texttt{iqimp(izone)} the value 0 and set 1312 the three velocity components (in 1313 $m.s^{-1}$) in 1314 \texttt{rcodcl(ifac,iu)}, 1315 \texttt{rcodcl(ifac,iv)} and 1316 \texttt{rcodcl(ifac,iw)} 1317 \end{list} 1318 \item finally he specifies for each gas inlet type 1319 the mixing rate \texttt{fment\index{fment}(izone)} and 1320 the temperature \texttt{tkent\index{tkent}(izone)} in Kelvin 1321 \end{list} 1322 1323 \item for the ``3 points'' diffusion flame module: 1324 \begin{list}{$\rightarrow$}{} 1325 \item the user can choose between the ``oxidiser 1326 inlet'' type marked out by 1327 \texttt{ientox\index{ientox}(izone)}=1 and the ``fuel 1328 inlet'' type marked out by 1329 \texttt{ientfu\index{ientfu}(izone)}=1 1330 \item concerning the input mass flow or the input 1331 velocity, the method is the same as for the 1332 EBU pre-mixed flame module 1333 \item finally, the user sets the temperatures 1334 \texttt{tinoxy\index{tinoxy}} for each oxidiser inlet 1335 and \texttt{tinfue\index{tinfue}}, for each fuel inlet 1336 1337{\em Note: In the standard version, only the cases with only one 1338 oxidising inlet type and one fuel inlet type 1339 can be treated. In particular, there must be 1340 only one input temperature for the oxidiser 1341 (\texttt{tinoxy}) and one input temperature for the 1342 fuel (\texttt{tinfuel}).} 1343 \end{list} 1344 1345 \item for the pulverised coal module: 1346 \begin{list}{$\rightarrow$}{} 1347 \item the inlet faces can belong to the ``primary 1348 air and pulverised coal inlet'' type, marked 1349 out by \texttt{ientcp\index{ientcp}(izone)}=1, or to 1350 the ``secondary or tertiary air inlet'' type, 1351 marked out by \texttt{ientat\index{ientat}(izone)}=1 1352 \item in a way which is similar to the process 1353 described in the framework of the EBU module, 1354 the user chooses for every inlet face to 1355 impose the mass flow or not (\texttt{iqimp(izone)}=1 or 1356 0). If the mass flow is imposed, the user 1357 must set the air mass flow value 1358 \texttt{qimpat\index{qimpat}(izone)}, its direction in 1359 \texttt{rcodcl(ifac,iu)}, \texttt{rcodcl(ifac,iv)} 1360 and \\ \texttt{rcodcl(ifac,iw)} and if 1361 \item incoming air temperature \texttt{timpat\index{timpat}(izone)} in 1362 Kelvin. If the velocity is imposed, he must 1363 set \texttt{rcodcl(ifac,iu)}, \\ 1364 \texttt{rcodcl(ifac,iv)} and \texttt{rcodcl(ifac,iw)}. 1365 1366 \item if the inlet belongs to the ``primary air and 1367 pluverised coal'' \texttt{type (ientcp(izone) = 1)} the 1368 user must also define for each coal type \texttt{icha}: 1369 the mass flow 1370 \texttt{qimpcp\index{qimpcp}(izone,icha)}, the 1371 granulometric distribution 1372 \texttt{distch\index{distch}(izone,icha,iclapc)} 1373 related to each class \texttt{iclacp}, and the 1374 injection temperature 1375 \texttt{timpcp\index{timpcp}(izone,icha)} 1376 1377 \end{list} 1378\end{list} 1379 1380%================================== 1381\subsubsection{Initialisation of the options of the variables} 1382%================================== 1383In the case of coal combustion, time averages, chronological records and logss follow-ups can be set in the Graphical User Interface (GUI) or in the subroutines \texttt{cs\_user\_combustion}. In the GUI, under the heading ``Calculation control'', additional variables appear in the list in the items ``Time averages'' and ``Profiles'', as well as in the item Volume solution control'', see \figurename~\ref{fig:t_average-coal} and \figurename~\ref{fig:V_control-coal}. 1384 1385\begin{figure}[!ht] 1386\begin{center} 1387\includegraphics[width=12cm]{gui_coal_time_average} 1388\caption{Calculation control - Time averages} 1389\label{fig:t_average-coal} 1390\end{center} 1391\end{figure} 1392 1393\begin{figure}[!ht] 1394\begin{center} 1395\includegraphics[width=13cm]{gui_coal_solution_control} 1396\caption{Calculation control - Volume solution control} 1397\label{fig:V_control-coal} 1398\end{center} 1399\end{figure} 1400 1401In this section, ``specific physics'' refers to gas combustion or 1402pulverised coal combustion. 1403 1404For gas combustion or if the GUI is not used for coal combustion, the 3 1405subroutines \texttt{cs\_user\_combustion} can be used to complete 1406\texttt{cs\_user\_parameters.f90} for the considered specific physics. These 1407subroutines are called at the calculation start. 1408They allow to: 1409\begin{list}{$\bullet$}{} 1410\item activate, for the variables which are specific to the activated specific 1411 physics module, chronological records at the probes defined in 1412 \texttt{cs\_user\_parameters.f90}.\\ 1413 Concerning the main variables (velocity, pressure, etc ...) the user 1414 must still complete \texttt{cs\_user\_parameters.f90} if he wants to get 1415 chronological records, printings in the log or chronological 1416 outputs. 1417 The variables which can be activated by the user for each specific 1418 physics are listed below. The solved variables (of variable indices 1419 \texttt{ivar}) and the properties of indices \texttt{iprop} (defined at 1420 the cell \texttt{iel} by \texttt{cpro\_prop(iel)} which is obtained 1421 by calling \texttt{field\_get\_val\_s(iprop, cpro\_prop)}) 1422 are listed below: 1423 \begin{list}{$\rightarrow$}{} 1424 \item EBU pre-mixed flame modelling: 1425 \begin{list}{-}{} 1426 \item Solved variables 1427 \begin{list}{\texttt{ivar} = }{} 1428 \item \texttt{isca(iygfm\index{iygfm})} fresh gas mass fraction 1429 \item \texttt{isca(ifm\index{ifm})} mixing rate 1430 \item \texttt{isca(ihm\index{ihm})} enthalpy, if transported 1431 \end{list} 1432 \item Properties \texttt{cpro\_prop(iel)} 1433 \begin{list}{\texttt{iprop} = }{} 1434 \item \texttt{itemp\index{itemp}} temperature 1435 \item \texttt{iym(1)\index{iym(1)}} fuel mass fraction 1436 \item \texttt{iym(2)\index{iym(2)}} oxidiser mass fraction 1437 \item \texttt{iym(3)\index{iym(3)}} product mass fraction 1438 \item \texttt{ickabs\index{ickabs}} absorption 1439 coefficient, when the radiation modelling is 1440 activated 1441 \item \texttt{it3m\index{it3m}} and \texttt{it4m\index{it4m}} 1442 ``$T^3$'' and ``$T^4$'' terms, when the radiation 1443 modelling is activated 1444 \end{list} 1445 \end{list} 1446 \item rapid complete chemistry diffusion flame modelling: 1447 \begin{list}{}{} 1448 \item everything is identical to the ``EBU'' case, except 1449 the fresh gas mass fraction which is replaced by the 1450 variance of the mixing rate 1451 \texttt{ivar=isca(ifp2m\index{ifp2m})} 1452 \end{list} 1453 \item pulverised coal modelling with 3 combustibles: 1454 \begin{list}{}{} 1455 \item {\em variables shared by the two phases}: 1456 \begin{list}{-}{} 1457 \item Solved variables 1458 \begin{list}{\texttt{ivar} = }{} 1459 \item \texttt{isca(ihm\index{ihm})}: gas-coal 1460 mixture enthalpy 1461 \item \texttt{isca(immel\index{immel})}: molar mass 1462 of the gas mixture 1463 \end{list} 1464 \end{list} 1465 \item {\em variables specific to the dispersed phase}: 1466 \begin{list}{-}{} 1467 \item Solved variables 1468 \begin{list}{\texttt{ivar} = }{} 1469 \item \texttt{isca(ixck\index{ixck}(icla))}: coke mass 1470 fraction related to the class \texttt{icla} 1471 \item \texttt{isca(ixch\index{ixch}(icla))}: reactive coal 1472 mass fraction related to the class \texttt{icla} 1473 \item \texttt{isca(inp\index{inp}(icla))}: number of 1474 particles of the class \texttt{icla} per kg of 1475 air-coal mixture 1476 \item \texttt{isca(ih2\index{ih2}(icla))}: mass enthalpy 1477 of the coal of class \texttt{icla}, if we are in 1478 permeatic conditions 1479 \end{list} 1480 \item Properties \texttt{cpro\_prop(iel)} 1481 \begin{list}{\texttt{iprop} = }{} 1482 \item \texttt{immel\index{immel}}: molar mass of the gas 1483 mixture 1484 \item \texttt{itemp2\index{itemp2}(icla)}: temperature of 1485 the particles of the class \texttt{icla} 1486 \item \texttt{irom2\index{irom2}(icla)}: density of 1487 the particles of the class \texttt{icla} 1488 \item \texttt{idiam2\index{idiam2}(icla)}: diameter of the 1489 particles of the class \texttt{icla} 1490 \item \texttt{igmdch\index{igmdch}(icla)}: disappearance 1491 rate of the reactive coal of the class \texttt{icla} 1492 \item \texttt{igmdv1\index{igmdv1}(icla)}: mass transfer 1493 caused by the release of light volatiles 1494 from the class \texttt{icla} 1495 \item \texttt{igmdv2\index{igmdv2}(icla)}: mass transfer 1496 caused by the release of heavy volatiles 1497 from the class \texttt{icla} 1498 \item \texttt{igmhet\index{igmhet}(icla)}: coke 1499 disappearance rate during the coke burnout 1500 of the class \texttt{icla} 1501 \item \texttt{ix2\index{ix2}(icla)}: solid mass fraction 1502 of the class \texttt{icla} 1503 \end{list} 1504 \end{list} 1505 \item {\em variables specific to the continuous phase}: 1506 \begin{list}{-}{} 1507 \item Solved variables 1508 \begin{list}{\texttt{ivar} = }{} 1509 \item \texttt{isca(if1m\index{if1m}(icha))}: mean value of 1510 the tracer 1 representing the light 1511 volatiles released by the coal \texttt{icha} 1512 \item \texttt{isca(if2m\index{if2m}(icha))}: mean value of 1513 the tracer 2 representing the heavy 1514 volatiles released by the coal \texttt{icha} 1515 \item \texttt{isca(if3m)\index{if3m}}: mean value of the 1516 tracer 3 representing the carbon released as 1517 CO during coke burnout 1518 \item \texttt{isca(if4pm\index{if4pm})}: variance of the 1519 tracer 4 representing the air 1520 \item \texttt{isca(if3p2m\index{if3p2m})}: variance of the 1521 tracer 3 1522 \end{list} 1523 \item Properties \texttt{cpro\_prop(iel)} 1524 \begin{list}{\texttt{iprop} = }{} 1525 \item \texttt{itemp\index{itemp1}}: temperature of the 1526 gas mixture 1527 \item \texttt{iym1(1)\index{iym1(1)}}: mass fraction of 1528 $CH_{X1m}$ (light volatiles) in the gas 1529 mixture 1530 \item \texttt{iym1(2)\index{iym1(2)}}: mass fraction of 1531 $CH_{X2m}$ (heavy volatiles) in the gas 1532 mixture 1533 \item \texttt{iym1(3)\index{iym1(3)}}: mass fraction of 1534 CO in the gas mixture 1535 \item \texttt{iym1(4)\index{iym1(4)}}: mass fraction of 1536 $O_2$ in the gas mixture 1537 \item \texttt{iym1(5)\index{iym1(5)}}: mass fraction of 1538 $CO_2$ in the gas mixture 1539 \item \texttt{iym1(6)\index{iym1(6)}}: mass fraction of 1540 $H_2O$ in the gas mixture 1541 \item \texttt{iym1(7)\index{iym1(7)}}: mass fraction of 1542 $N_2$ in the gas mixture 1543 \end{list} 1544 \end{list} 1545 \end{list} 1546 \end{list} 1547 1548 \item set the relaxation coefficient of the density \texttt{srrom}, with \\ 1549$\rho^{n+1}=\texttt{srrom}*\rho^{n}+(1-\texttt{srrom})\rho^{n+1}$\\ 1550(the default value is \texttt{srrom\index{srrom}} = 0.8. At the 1551 beginning of a calculation, a sub-relaxation of 0.95 may reduce 1552 the numerical ``shocks''). 1553 1554 \item set the dynamic viscosity \texttt{diftl0}. By default 1555 \texttt{diftl0\index{diftl0}}= 4.25 $kgm^{-1}s^{-1}$ 1556(the dynamic diffusivity being the ratio between the thermal 1557 conductivity $\lambda$ and the mixture specific heat $C_p$ in the 1558 equation of enthalpy). 1559 1560 \item set the value of the constant \texttt{cebu\index{cebu}} of the Eddy Break 1561 Up model (only in \texttt{cs\_user\_combustion}. By default \texttt{cebu}=2.5) 1562\end{list} 1563 1564%================================== 1565\subsection{Heavy fuel oil combustion module} 1566%================================== 1567%================================== 1568\subsubsection{Initialisation of transported variables} 1569%================================== 1570To initialise or modify (in case of a continuation) values of transported 1571variables and of the time step, the standard subroutine \texttt{cs\_user\_initialization} is used. 1572 1573Physical properties are stored using the \texttt{cs\_field} API (cell center). For instance, to obtain \texttt{rom(iel)}, 1574the mean density (in $kg.m^{-3}$), one must declare a \texttt{ncelet} array \texttt{cpro\_rom} and then call 1575\texttt{call field\_get\_val\_s(icrom, cpro\_rom)}.\\ 1576Physical properties (\texttt{rom, viscl, cp, ...}) are computed in \texttt{ppphyv} and are not to be modified here. 1577 1578The \texttt{cs\_user\_initialization-fuel.f90} example illustrates how the user 1579may initialise quantities related to gaseous species and droplets compositions 1580in addition to the chosen turbulent model. 1581 1582%================================== 1583\subsubsection{Boundary conditions} 1584%================================== 1585Boundary conditions are defined as usual on a per-face basis in 1586\texttt{cs\_user\_boundary\_conditions}. They may be assigned in two ways: 1587\begin{list}{.}{} 1588\item for ``standard'' boundary conditions (inlet, free outlet, wall, symmetry): a code is defined in the array \texttt{itypfb} (of dimensions equal to the number of boundary faces). This code will then be used by a non-user subroutine to assign the conditions. 1589\item for ``non-standard'' conditions: see details given in 1590 \texttt{cs\_user\_boundary\_conditions-fuel.f90} example. 1591\end{list} 1592 1593%================================== 1594\subsection{Radiative thermal transfers in semi-transparent gray media} 1595%================================== 1596%================================== 1597\subsubsection{Initialisation of the radiation main parameters} 1598%================================== 1599 1600The main radiation parameters can be initialise in the Graphical User Interface (GUI) or in the user subroutine \texttt{cs\_user\_radiative\_transfer\_param}. In the GUI, under the heading ``Thermophysical models'', when one of the two thermal radiative transfers models is selected, see \figurename~\ref{fig:0_ray}, additional items appear. The user is asked to choose the number of directions for angular discretisation, to define the absorption coefficient and select if the radiative calculation are restarted or not, 1601see \figurename~\ref{fig:1_ray} and \figurename~\ref{fig:3_ray}. When ``Advanced options'' is selected for both models \figurename~\ref{fig:2_ray} or \figurename~\ref{fig:4_ray} appear, the user must fill the resolution frequency and verbosity levels. In addition, the activation of the radiative transfer leads to the creation of an item ``Surface solution control'' under the heading ``Calculation control'', 1602see \figurename~\ref{fig:5_ray}, where radiative transfer variables can be selected to appear in the output log. 1603 1604\begin{figure}[ht] 1605\begin{center} 1606\includegraphics[width=0.9\textwidth]{gui_rad_transf_do_params} 1607\caption{Radiative transfers - parameters of the DO method} 1608\label{fig:1_ray} 1609\end{center} 1610\end{figure} 1611 1612\begin{figure}[ht] 1613\begin{center} 1614\includegraphics[width=7cm]{gui_rad_transf_do_advanced} 1615\caption{Radiative transfers - advanced parameters of the DO method} 1616\label{fig:2_ray} 1617\end{center} 1618\end{figure} 1619 1620\begin{figure}[ht] 1621\begin{center} 1622\includegraphics[width=0.9\textwidth]{gui_rad_transf_p1_params} 1623\caption{Radiative transfers - parameters of the P-1 model} 1624\label{fig:3_ray} 1625\end{center} 1626\end{figure} 1627 1628\begin{figure}[ht] 1629\begin{center} 1630\includegraphics[width=7cm]{gui_rad_transf_p1_advanced} 1631\caption{Radiative transfers - advanced parameters of the P-1 model} 1632\label{fig:4_ray} 1633\end{center} 1634\end{figure} 1635 1636\begin{figure}[ht] 1637\begin{center} 1638\includegraphics[width=0.9\textwidth]{gui_rad_transf_post_output} 1639\caption{Calculation control - Radiative transfers post-processing output} 1640\label{fig:5_ray} 1641\end{center} 1642\end{figure} 1643 1644If the GUI is not used, \texttt{cs\_user\_radiative\_transfer\_param} is one of the two subroutine which must be completed by the user for all 1645calculations including radiative thermal transfers. It is called only during the calculation initialisation. It is composed of three headings. The first one is dedicated to the activation 1646of the radiation module, only in the case of classic physics. \\ 1647{\em WARNING: when a calculation is ran using a specific physics module, 1648this first heading must not be completed. The radiation module is then 1649activated or not, according to the parameter file related to the considered 1650specific physics.} \\ 1651 1652\noindent 1653In the second heading the basic parameters of the radiation module are indicated.\\ 1654Finally, the third heading deals with the selection of the 1655post-processing graphic outputs. The variables to treat are splitted 1656into two categories: the volumetric variables and those related to the 1657boundary faces.\\ 1658 1659\noindent 1660For more details about the different parameters, the user may refer to the 1661keyword list (\S~\ref{sec:prg_motscles}). 1662 1663 1664%================================== 1665\subsubsection{Radiative transfers boundary conditions} 1666%================================== 1667These informations can be filled by the user through the Graphical User Interface 1668(GUI) or by using the subroutine \texttt{cs\_user\_radiative\_transfer\_bcs.c} (called every time step). If the 1669interface is used, when one of the ``Radiative transfers'' options is selected in 1670\figurename~\ref{fig:gui_thermal_scalar}, it activates specific boundary conditions each time 1671a ``Wall'' is defined, see \figurename~\ref{fig:6_ray}. The user can then choose 1672between 3 cases. The parameters the user must specify are displayed for one of 1673them in \figurename~\ref{fig:7_ray}. 1674 1675\begin{figure}[ht] 1676\begin{center} 1677\includegraphics[width=0.9\textwidth]{gui_rad_transf_wall_model} 1678\caption{Boundary conditions - choice of wall thermal radiative transfers} 1679\label{fig:6_ray} 1680\end{center} 1681\end{figure} 1682 1683\begin{figure}[ht] 1684\begin{center} 1685\includegraphics[width=0.9\textwidth]{gui_rad_transf_wall_params} 1686\caption{Boundary conditions - example of wall thermal radiative transfer} 1687\label{fig:7_ray} 1688\end{center} 1689\end{figure} 1690 1691When the GUI is not used, \texttt{cs\_user\_radiative\_transfer\_bcs.c} is needed for 1692every calculation which includes radiative thermal transfers. It is used to give all the 1693necessary parameters concerning, in the one case, the wall temperature 1694calculation, and in the other, the coupling between the thermal 1695scalar (temperature or enthalpy), and the radiation module at the 1696calculation domain boundaries. It must be noted that the boundary conditions 1697concerning the thermal scalar which may have been defined in the 1698GUI or in subroutine \texttt{cs\_user\_boundary\_conditions} will be modified by 1699the radiation module according to the data given in \texttt{cs\_user\_radiative\_transfer\_bcs} (cf. \S\ref{sec:fvm_selector}).\\ 1700A boundary condition type stored in the array \texttt{isothp} is associated with 1701each boundary face. There are five different types: 1702 1703\begin{list}{$\bullet$}{} 1704 1705\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_GREY}}: grey of black wall, with 1706 temperature defined by main (fluid) boundary conditions, 1707 1708\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_GREY\_EXTERIOR\_T}}: 1709 grey or black wall, calculation of the temperature by means of a flux balance 1710 with an exterior temperature, 1711 1712\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_REFL\_EXTERIOR\_T}}: 1713 for a reflecting wall face, calculation of the 1714 temperature by means of a flux balance. This is equivalent to 1715 using \textbf{CS\_BOUNDARY\_RAD\_WALL\_GRAY\_EXTERIOR\_T}, 1716 with zero emissivity. 1717 1718\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_GRAY\_COND\_FLUX}}: 1719 grey or black wall face to which a conduction 1720 flux is imposed, 1721 1722\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_REFL\_COND\_FLUX}}: 1723 reflecting wall face to which a conduction 1724 flux is imposed, which is equivalent to impose this flux directly 1725 to the fluid. 1726 1727\item \texttt{\textbf{ifinfe}}: for an open boundary (inlet or outlet) or symmetry face, 1728 simulate an infinite extrusion by applying a Neumann condition to 1729 the radiation equations, 1730 1731\end{list} 1732 1733\noindent 1734Depending on the selected boundary condition type at every wall face, 1735the code needs to be given some additional information: 1736 1737\begin{list}{$\bullet$}{} 1738 1739\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_GRAY}}: the array \texttt{epsp} must 1740 be completed with the emissivity value (positive). 1741 1742\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_GREY\_EXTERIOR\_T}}: wall emissivity 1743 (strictly positive, in \texttt{epsp}), thickness (in \texttt{epap}), 1744 thermal conductivity (in \texttt{xlamp}) and external temperature 1745 (in \texttt{textp}) in order to calculate a conduction flux across the wall. 1746 1747\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_REFL\_EXTERIOR\_T}}: 1748 wall thickness (in \texttt{epap}) and thermal 1749 conductivity (in \texttt{xlamp}) and an external temperature 1750 (in \texttt{textp}). 1751 This is equivalent to using \textbf{CS\_BOUNDARY\_RAD\_WALL\_GREY\_EXTERIOR\_T} 1752 with a wall emissivity of 0, and an exchange coefficient-type boundary 1753 condition, where the exchange coefficient is equal to 1754 \texttt{xlamp} / \texttt{textp}. 1755 1756\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_GRAY\_COND\_FLUX}}: 1757 wall emissivity (in \texttt{epsp}) and conduction 1758 flux (in $W/m^2$ whatever the thermal scalar, enthalpy or temperature) in 1759 the array \texttt{rcodcl}. The value of \texttt{rcodcl} is positive when the 1760 conduction flux is directed from the inside of the fluid domain to the 1761 outside (for instance, when the fluid heats the walls). If the 1762 conduction flux is null, the wall is adiabatic. 1763 1764\item \texttt{\textbf{CS\_BOUNDARY\_RAD\_WALL\_REFL\_COND\_FLUX}}: 1765 must be given: the conduction flux (in $W/m^2$ whatever the thermal 1766 scalar) in the array \texttt{rcodcl}. The value of \texttt{rcodcl} is 1767 positive when the conduction flux is directed from the inside of the 1768 fluid domain to the outside (for instance, when the fluid heats the 1769 walls). If the conduction flux is null, the wall is adiabatic. The flux 1770 received by \texttt{rcodcl} is directly imposed as boundary condition for 1771 the fluid. 1772 1773\end{list} 1774 1775%================================== 1776\subsubsection{Absorption coefficient of the medium, boundary conditions 1777 for the luminance and calculation of the net radiative flux} 1778%================================== 1779 1780When the absorption coefficient is not constant, the subroutine 1781\texttt{cs\_user\_rad\_transfer\_absorption} is called instead at each time 1782step. It is composed of three parts. In the first one, the user 1783must provide the absorption coefficient of the medium in the array CK, 1784for each cell of the fluid mesh. By default, the absorption coefficient 1785of the medium is 0, which corresponds to a transparent medium.\\ 1786 1787{\em WARNING: when a specific physics is activated, it is forbidden to 1788give a value to the absorption coefficient in this subroutine. In this 1789case, the coefficient is either calculated automatically, or provided by the user {\em via} a 1790thermo-chemical parameter file (dp\_C3P or dp\_C3PSJ for gas combustion, 1791and dp\_FCP for pulverised coal combustion).}\\ 1792 1793\noindent 1794The two following parts of this subroutine concern a more advanced use 1795of the radiation module. It is about imposing boundary conditions to the 1796equation of radiative transfer and net radiative flux calculation, in 1797coherence with the luminance at the boundary faces, when the user wants 1798to give it a particular value. In most cases, the given examples do not 1799need to be modified. 1800 1801%================================== 1802\subsection{Conjugate heat transfer} 1803%================================== 1804 1805%======================================== 1806\subsubsection{Thermal module in a 1D wall} 1807%======================================== 1808 1809\noindent 1810\textit{subroutine called at every time step} 1811 1812This subroutine takes into account the wall-affected thermal inertia. 1813 Some boundary faces are treated as a solid wall with a given thickness, on 1814 which the code resolves a one-dimensional equation for the heat conduction. 1815 The coupling between the 1D module and the fluid works in a similar way to 1816 the coupling with the \syrthes. By construction, the user is not able to 1817 account for the heat transfer between different parts of the wall. A physical 1818 analysis of each problem, case by case is required in order to evaluate the relevance 1819 of its usage by way of a report of the simple conditions (temperature, zero-flux 1820 ) or a coupling with \syrthes.\\ 1821 1822The use of this code requires that the thermal scalar is 1823defined as (\texttt{iscalt}$>0$). 1824 1825{\em WARNING: The 1D thermal module is developed assuming the thermal scalar 1826 as a temperature. If the thermal scalar is an enthalpy, the code calls the 1827 enthalpy to temperature conversion as defined by the model defaults, 1828 or by the user in \texttt{cs\_user\_physical\_properties} for each 1829 transfer of data between the fluid and the wall in order to convert the 1830 enthalpy to temperature and vice-versa. If the thermal 1831 variable is the total (compressible) energy, the thermal module will not work.} 1832 1833%================================== 1834\subsubsection{Internal Fluid-Thermal coupling} 1835%================================== 1836 1837When at least one volume zone is defined as being solid 1838(see \figurename\ref{fig:internal_coupling_vol_solid}), scalar variables (especially 1839thermal scalar variables) may be solved in a fully coupled manner across the fluid 1840and solid domains. 1841 1842For this purpose, the ``Internal coupling'' should be activated for the desired variables 1843in the matching tab of the ``Coupling parameters'' page, as shown in figure 1844\figurename\ref{fig:internal_coupling_param}). This section should appear when 1845at least one volume zone is defined as solid. 1846 1847\begin{figure}[ht] 1848\begin{center} 1849\includegraphics[width=0.9\textwidth]{gui_internal_coupling_vol_zone} 1850\caption{Solid volume zone definition} 1851\label{fig:internal_coupling_vol_solid} 1852\end{center} 1853\end{figure} 1854 1855\begin{figure}[ht] 1856\begin{center} 1857\includegraphics[width=0.9\textwidth]{gui_internal_coupling} 1858\caption{Conjugate heat transfer: internal coupling} 1859\label{fig:internal_coupling_param} 1860\end{center} 1861\end{figure} 1862 1863%================================== 1864\subsubsection{Fluid-Thermal coupling with \syrthes} 1865%================================== 1866Coupling \CS with \syrthes for conjugate heat transfer can be defined through 1867the GUI or the \texttt{cs\_syrthes\_coupling} user function. 1868To set such a coupling in the GUI, a thermal scalar must be 1869selected first in the item ``Thermal scalar'' under the heading ``Thermophysical models''. 1870At least one wall boundary condition must be set to ``SYRTHES coupling'' type, and 1871the name of the associated \syrthes instance (i.e. base directory name of the associated 1872solid case definition) be set, as shown in, \figurename\ref{fig:syrthes_bc}. 1873The ``Syrthes coupling'' tab will then be available in the ``Coupling parameters'' 1874section (see \figurename\ref{fig:syrthes}), fo further advanced or global settings. 1875The zones where the coupling occurs must be defined and a projection axis can be 1876specified in case of 2D coupling. 1877 1878\begin{figure}[ht] 1879\begin{center} 1880\includegraphics[width=0.9\textwidth]{gui_syrthes_coupling_bc} 1881\caption{Boundary conditions - coupling with \syrthes} 1882\label{fig:syrthes_bc} 1883\end{center} 1884\end{figure} 1885 1886\begin{figure}[ht] 1887\begin{center} 1888\includegraphics[width=0.9\textwidth]{gui_syrthes_coupling} 1889\caption{Coupling parameters - coupling with \syrthes} 1890\label{fig:syrthes} 1891\end{center} 1892\end{figure} 1893 1894If the function \texttt{cs\_user\_syrthes\_coupling} is used, the user must 1895specify the arguments passed to the '\texttt{cs\_syr\_coupling\_define}' function. 1896 These arguments are: 1897\begin{list}{-}{} 1898 \item \texttt{syrthes\_name} is the matching \syrthes application name (useful only when more than one \syrthes and one \CS domain are present), 1899 \item \texttt{boundary\_criteria} is the surface selection criteria, 1900 \item \texttt{volume\_criteria} is the volume selection criteria, 1901 \item \texttt{projection\_axis}: ' ' if the user wishes to use a 3D standard coupling, 1902 or specify '$x$', '$y$', or '$z$' as the projection axis if a 2D coupling with \syrthes is used, 1903 \item \texttt{verbosity} is the verbosity level. 1904 \item \texttt{visualization} is the visualization level. 1905\end{list} 1906Examples are provided in \texttt{cs\_user\_coupling.c}. 1907 1908The user may also define global coupling options relative to the handling of 1909time-stepping, by adapting the example \texttt{cs\_user\_coupling} 1910in the \texttt{cs\_user\_coupling.c} file. In the case of multiple couplings, 1911these options are global to all \syrthes and \CS couplings. 1912 1913%================================== 1914\subsection{Particle-tracking (Lagrangian) Module} 1915%================================== 1916 1917%================================== 1918\subsubsection{General information}\label{sec:over-lag} 1919%================================== 1920 1921\begin{itemize} 1922 1923\item[-] The particle-tracking (or Lagrangian) module enables the simulation of poly-dispersed particulate flows, by calculating the trajectories of individual particles, mainly characterized by their diameter and density (if no heat nor mass transfer between particle and fluid are activated). 1924 1925\item[-] The standard use of the particle-tracking module follows the \textbf{Moments/PDF approach}: the instantaneous properties of the underlying flow needed to calculate the particle motion are reconstructed from the averaged values (obtained by Reynolds-Averaged Navier-Stokes simulation) by using stochastic processes. The statistics of interest are then obtained through Monte-Carlo simulation. 1926 1927\item[-] As a consequence, is is important to emphasize that the most important (and physically meaningful) results of a particle-tracking calculation following the Moments/PDF approach are \mbox{\textbf{statistics}}. Volume and surface statistics, steady or unsteady, can be calculated. Individual particle trajectories (as 1D, \textit{EnSight}-readable cases) and displacements (as \textit{EnSight}-readable animations) can also be provided, but only for illustrative purposes. 1928 1929\end{itemize} 1930 1931%================================== 1932\subsubsection{Activating the particle-tracking module}\label{sec:acti-lag} 1933%================================== 1934 1935The activation of the particle-tracking module is performed either: 1936% 1937\begin{itemize} 1938 \item [$\bullet$] in the Graphical User Interface (GUI): \texttt{Calculation features} $\rightarrow$ \texttt{Thermophysical models} $\rightarrow$ \texttt{Eulerian-Lagrangian multi-phase treatment}~$\rightarrow$ ~\texttt{particles and droplets tracking} 1939 \item [$\bullet$] or in the user function \texttt{cs\_user\_lagr\_model}. 1940\end{itemize} 1941 1942%================================== 1943\subsubsection{Basic guidelines for standard simulations} 1944%================================== 1945 1946Except for cases in which the flow conditions depend on time, it is generally recommended to perform a first Lagrangian calculation whose aim is to reach a steady-state (i.e. to reach a time starting from which the relevant statistics do not depend on time anymore). In a second step, a calculation restart is done to calculate the statistics. When the single-phase flow is steady and the particle volume fraction is low enough to neglect the particles influence on the continuous phase behaviour, it is recommended to perform a Lagrangian calculation on a frozen field.\\ 1947 1948It is then possible to calculate steady-state volumetric statistics and to give a statistical weight higher than 1 to the particles, in order to reduce the number of simulated (``numerical'') particles to treat while keeping the right concentrations. Otherwise, when the continuous phase flow is steady, but the two-coupling coupling must be taken into consideration, it is still possible to activate steady statistics. \\ 1949When the continuous phase flow is unsteady, it is no longer possible to use steady statistics. To have correct statistics at every moment in the whole calculation domain, it is imperative to have an established particle seeding and it is recommended (when it is possible) not to impose statistical weights different from the unity. \\ 1950 1951Finally, when the so-called complete model is used for turbulent dispersion modelling, the user must make sure that the volumetric statistics are directly used for the calculation of the locally undisturbed fluid flow field.\\ 1952 1953When the thermal evolution of the particles is activated, the associated particulate scalars are always the inclusion temperature and the locally undisturbed fluid flow temperature expressed in degrees Celsius, whatever the thermal scalar associated with the continuous phase is ({\em i.e.} temperature or enthalpy). If the 1954thermal scalar associated with the continuous phase is the temperature in Kelvin, the unit is converted automatically into Celsius. If the thermal scalar associated with the continuous phase is the enthalpy, a \texttt{temperature} property or postprocessing 1955field must be defined. In all cases, the thermal backward coupling of the dispersed phase on the continuous phase is adapted to the thermal scalar transported by the fluid. 1956 1957 1958%================================== 1959\subsubsection{Prescribing the main modelling parameters}\label{sec:Ini-lag} 1960%================================== 1961 1962\minititre{Use of the GUI} 1963 1964In the GUI, the selection of the Lagrangian module activates the heading \texttt{Particle and droplets tracking} in the tree menu. The initialization is performed in the three items included in this heading: 1965% 1966\begin{itemize} 1967 \item [$\bullet$] \texttt{Global settings}. The user defines in this item the kind of Euler/Lagrange multi-phase treatment, the main parameters, and the specific physics associated with the particles, see ~\figurename~\ref {fig:Ini-Lag1} to \figurename\ref {fig:Ini-Lag3}. 1968 \item [$\bullet$] \texttt{Statistics}. The user can select the volume and boundary statistics to be post-processed. 1969 \item [$\bullet$] \texttt{Output}. An additional entry in the postprocessing section allows defining the output frequency and post-processing options for particles and selecting the variables that will appear in the log. 1970\end{itemize} 1971 1972\begin{figure}[ht] 1973\begin{center} 1974\includegraphics[width=0.9\textwidth]{gui_lagr_global_settings} 1975\caption{Lagrangian module - View of the \texttt{Global Settings} page} 1976\label{fig:Ini-Lag1} 1977\end{center} 1978\end{figure} 1979 1980\begin{figure}[ht] 1981\begin{center} 1982\includegraphics[width=11cm]{gui_lagr_statistics} 1983\caption{Lagrangian module - statistics} 1984\label{fig:Ini-Lag4} 1985\end{center} 1986\end{figure} 1987 1988\begin{figure}[ht] 1989\begin{center} 1990\includegraphics[width=11cm]{gui_lagr_output} 1991\caption{Lagrangian module - output} 1992\label{fig:Ini-Lag5} 1993\end{center} 1994\end{figure} 1995 1996\minititre{Use of the subroutine \texttt{cs\_user\_lagr\_model}} 1997 1998\noindent 1999When the GUI is not used, \texttt{cs\_user\_lagr\_model} must be completed. This function 2000gathers in different headings all the keywords which are 2001necessary to configure the Lagrangian module. The different headings refer to: 2002\begin{list}{$\bullet$}{} 2003\item the global configuration parameters 2004\item the specific physical models describing the particle behaviour 2005\item the backward coupling (influence of the dispersed phase on the 2006 continuous phase) 2007\item the numerical parameters 2008\item the volumetric statistics 2009\item the boundary statistics 2010\end{list} 2011% 2012\noindent 2013For more details about the different parameters, the user may refer to the 2014keyword list (\S~\ref{sec:prg_motscles_lagr}).%FIXME 2015 2016 2017%================================== 2018\subsubsection{Prescribing particle boundary conditions} 2019%================================== 2020In the framework of the multiphase Lagrangian modelling, the management of the boundary conditions concerns the particle behaviour when there is an interaction between its trajectory and a boundary face. These boundary conditions may be imposed independently of those concerning the Eulerian fluid phase (but they are of course generally consistent). The boundary condition zones are actually redefined by the Lagrangian module (cf. \S\ref{sec:fvm_selector}), and a type of particle behaviour is associated with each one. The boundary conditions related to particles can be defined in the Graphical User Interface (GUI) or in the \texttt{cs\_user\_lagr\_boundary\_conditions.c} file. More advanced user-defined boundary conditions can be prescribed in the \texttt{cs\_user\_lagr\_in function} from \texttt{cs\_user\_lagr\_particle.c}. 2021 2022\minititre{Use of the GUI} 2023 2024 In the GUI, selecting the Lagrangian module in the activates the item \texttt{Particle boundary conditions} under the heading \texttt{Boundary conditions} in the tree menu. Different options are available depending on the type of standard boundary conditions selected (wall, inlet/outlet, etc...), 2025 see \figurename~\ref{fig:CL-Lag}. 2026 2027\begin{figure}[ht] 2028\begin{center} 2029\includegraphics[width=0.9\textwidth]{gui_lagr_bc} 2030\caption{Lagrangian module - boundary conditions} 2031\label{fig:CL-Lag} 2032\end{center} 2033\end{figure} 2034 2035%================================== 2036\subsubsection{Advanced particle-tracking set-up} 2037%================================== 2038 2039In this section, some information is provided for a more advanced numerical set-up of a particle-tracking simulation. 2040 2041\minititre{User-defined stochastic differential equations} 2042%------------------------------------------------- 2043 2044\noindent 2045An adaptation in the \texttt{cs\_user\_lagr\_sde} function is required if 2046supplementary user variables are added to the particle state vector. 2047This function is called at each Lagrangian sub-step. 2048 2049\noindent 2050The integration of the stochastic differential equations associated with 2051supplementary particulate variables is done in this function. \\ 2052When the integration scheme of the stochastic differential equations is 2053a first-order (\texttt{nordre} = 1), this subroutine is called once every 2054Lagrangian iteration, if it is a second-order (\texttt{nordre} = 2), it is called 2055twice. \\ 2056 2057\noindent 2058The solved stochastic differential equations must be written in the 2059form: 2060\begin{displaymath} 2061\frac{d \Phi_p}{dt} \,=\, - \frac{\Phi_p - \Pi}{\tau_\phi} 2062\end{displaymath} 2063where $\Phi_p$ is the I\textit{th} supplementary user variable, 2064$\tau_\phi$ is a quantity homogeneous to a characteristic time, and $\Pi$ is 2065a coefficient which may be expressed as a function of the other 2066particulate variables. \\ 2067In order to do the integration of this equation, the following 2068parameters must be provided: 2069\begin{list}{-}{} 2070\item $\tau_\phi$, equation characteristic time 2071 every particle, 2072\item $\Pi$ , equation coefficient. If the 2073 integration scheme is a first-order, then $\Pi$ is expressed as a 2074 function of the particulate variables at the previous iteration, 2075 stored in the array \texttt{eptpa}. If the chosen scheme is a second-order, 2076 then $\Pi$ is expressed at the first call of the function 2077 (prediction step) as a function of the variables at the 2078 previous iteration, then at the second call 2079 (correction step) as a function of the predicted variables. 2080\end{list} 2081 2082\noindent 2083If necessary, the thermal characteristic time $\tau_c$, whose calculation can be modified by the user in the function 2084\texttt{cs\_user\_lagr\_rt}. 2085 2086\minititre{User-defined particle relaxation time} 2087%------------------------------------------------- 2088 2089\noindent 2090The particle relaxation time may be modified in the \texttt{cs\_user\_lagr\_rt} function according to the chosen formulation of the drag coefficient. The particle relaxation time, modified or not by the user, is available in the array \texttt{taup}. 2091 2092\minititre{User-defined particle thermal characteristic time} 2093%------------------------------------------- 2094 2095\noindent 2096The particle thermal characteristic time may be modified in the \texttt{cs\_user\_lagr\_rt\_t} function according to the chosen correlation for the calculation of the 2097Nusselt number. This function is called at each Lagrangian sub-step. 2098 2099 2100%================================== 2101\subsection{Compressible module} 2102%================================== 2103 2104When the compressible module\footnote{For more details concerning the 2105compressible version, the user may refer to the theory guide \cite{theory} and the document ``Implantation 2106d'un algorithme compressible dans \CS'', Rapport EDF 2003, 2107HI-83/03/016/A, P. Mathon, F. Archambeau et J.-M. H\'erard.} is 2108activated, it is recommended to: 2109\begin{list}{-}{} 2110 \item use the option ``time step variable in time and uniform in 2111 space'' (\texttt{idtvar}=1) with a maximum Courant number of 0.4 2112 (\texttt{coumax}=0.4): these choices must be written in \texttt{cs\_user\_parameters.f90} 2113 or specified with the GUI. 2114 \item keep the convective numerical schemes proposed by default (\textit{i.e.}: upwind scheme). 2115\end{list} 2116With the compressible algorithm, the specific total energy is a new solved variable 2117\texttt{isca(ienerg)}). The temperature variable deduced from the specific 2118total energy variable is \texttt{isca(itempk)} for the compressible module.\\ 2119Initialisation of the options of the variables, boundary conditions, initialisation of the variables and 2120management of variable physical properties can be done with the GUI. We describe below the subroutines 2121the user has to fill in without the GUI. 2122 2123%================================== 2124\subsubsection{ Initialisation of the options of the variables} 2125%================================== 2126\label{prg_uscfx12}% 2127\noindent 2128\textit{Subroutines called at each time step.} 2129 2130When the GUI is not being used, the subroutines \texttt{uscfx1} and \texttt{uscfx2} in \texttt{cs\_user\_parameters.f90} 2131must be completed by the user. 2132 2133\texttt{uscfx1} allows to specify: 2134\begin{list}{-}{} 2135\item \texttt{ieos}: equation of state (only perfect gas with a constant adiabatic coefficient, 2136 \texttt{ieos=1} is available, but the user can complete the subroutine 2137 \texttt{cs\_cf\_thermo}, which is not a user subroutine, to add new equations of state). 2138\item \texttt{call field\_set\_key\_int(ivarfl(isca(itempk)), kivisl, ...)}: molecular thermal conductivity, constant (\texttt{-1}) or variable (\texttt{0}). 2139\item \texttt{iviscv}: volumetric molecular viscosity, constant (\texttt{0}) or variable (\texttt{1}). 2140\end{list} 2141 2142\texttt{uscfx2} allows to specify: 2143\begin{list}{-}{} 2144 \item \texttt{ivivar}: molecular viscosity, constant (\texttt{0}) or variable (\texttt{1}). 2145 \item \texttt{field\_set\_key\_double(icavr(isca(itempk)), kvisl0, ...)}: reference molecular thermal conductivity. 2146 \item \texttt{viscv0}: reference volumetric molecular viscosity. 2147 \item \texttt{xmasmr}: molar mass of the perfect gas (\texttt{ieos=1}). 2148 \item \texttt{icfgrp}: specify if the hydrostatic equilibrium must be accounted for in the 2149 boundary conditions. 2150\end{list} 2151 2152%================================== 2153\subsubsection{Management of the boundary conditions} 2154%================================== 2155 2156\noindent 2157\textit{Subroutine called at each time step.} 2158 2159When running the compressible module without a GUI, the \texttt{cs\_user\_boundary\_conditions} subroutine can be used to define specific boundary conditions 2160(see the \texttt{cs\_user\_boundary\_conditions-compressible} file in the directory \texttt{EXAMPLES} 2161for examples of boundary conditions with the compressible module). 2162 2163With the compressible module, the following types of boundary condition are avaliable: 2164 2165\begin{list}{-}{} 2166 \item Inlet/outlet for which velocity and two thermodynamics variables are known. 2167 \item Subsonic inlet with imposed total pressure and total energy. 2168 \item Subsonic outlet with imposed static pressure. 2169 \item Supersonic outlet. 2170 \item Wall (adiabatic or not). 2171 \item Symmetry. 2172\end{list} 2173 2174It is advised to only use these predefined boundary conditions type for the compressible module. 2175 2176%================================== 2177\subsubsection{Initialisation of the variables} 2178%================================== 2179 2180\noindent 2181\textit{Subroutine called only at the initialisation of the calculation} 2182 2183When the GUI is not used, the subroutine \texttt{cs\_user\_initialization} is used 2184initialize the velocity, turbulence and passive scalars (see 2185the \texttt{cs\_user\_initialization-compressible} file in the directory \texttt{EXAMPLES} 2186for examples of initialisations with the compressible module). Concerning 2187pressure, density, temperature and specific total energy, only 2 variables out 2188of these 4 are independent. The user may then initialise the desired variable pair 2189(apart from temperature-energy) and the two other variables will be 2190calculated automatically by giving the right value to the variable 2191\texttt{ithvar} used for the call to the \texttt{cs\_cf\_thermo} routine. 2192 2193%================================== 2194\subsubsection{Management of variable physical properties} 2195%================================== 2196 2197\noindent 2198\textit{Subroutine called at each time step.} 2199 2200Without the GUI, all of the laws governing the physical properties of the fluid 2201(molecular viscosity, molecular volumetric viscosity, molecular thermal conductivity and 2202molecular diffusivity of the user-defined scalars) can be specified in the subroutine \texttt{usphyv} of 2203the \texttt{cs\_user\_physical\_properties} file, 2204which is then called at each time step. This subroutine replaces and is similar to \texttt{usphyv}. 2205 2206The user should check that the defined laws are valid for 2207the whole variation range of the variables. Moreover, as only the perfect gas with a constant 2208adiabatic coefficient equation of state is available, it is not advised to give a law for the isobaric 2209specific heat without modifying the equation of state in the subroutine \texttt{cs\_cf\_thermo} which is not 2210a user subroutine. 2211 2212%================================== 2213\subsection{Management of the electric arcs module} 2214%================================== 2215 2216%================================== 2217\subsubsection{Activating the electric arcs module}\label{sec:acti-lag} 2218%================================== 2219 2220The electric arcs module is activated either: 2221% 2222\begin{itemize} 2223 \item [$\bullet$] in the Graphical User Interface (GUI): \texttt{Calculation features} $\rightarrow$ \texttt{Electrical models} 2224 \item [$\bullet$] or in the user subroutine \texttt{usppmo}, by setting the \texttt{ielarc} or \texttt{ieljou} parameter to a non-null value. 2225\end{itemize} 2226 2227%================================== 2228\subsubsection{Initialisation of the variables} 2229%================================== 2230 2231\noindent 2232\textit{Subroutine called only at initialisation of the calculation} 2233 2234The subroutine \texttt{cs\_user\_initialization} allows the user to initialise some of the specific physics variables prompted via \texttt{usppmo}. It is called only during the initialisation of the calculation. As usual,the user has access to many geometric variables so 2235 that the zones can be treated separately if needed. 2236 2237The values of potential and its constituents are initialised if required. 2238 2239It should be noted that the enthalpy is relevant. 2240 2241\begin{list}{-}{} 2242\item For the electric arcs module, the enthalpy value is taken from the temperature 2243 of reference \texttt{t0} (given in \texttt{cs\_user\_parameters.c}) 2244 from the temperature-enthalpy tables supplied in the data file \texttt{dp\_ELE.} 2245 The user must not intervene here. 2246\item For the Joule effect module, the value of enthalpy must be specified by the user. 2247 Examples of temperature to enthalpy conversion are given in 2248 \texttt{cs\_user\_physical\_properties.c}). If not defined, a simple default 2249 law is used ($H = C_p T$). 2250\end{list} 2251 2252%================================== 2253\subsubsection{Variable physical properties} 2254%================================== 2255 2256All the laws of the variation of physical data of the fluid are written (when necessary) 2257in the subroutine \texttt{cs\_user\_physical\_properties}. It is called at each time step. 2258 2259{\em WARNING: For the electric module, it is here that all the physical variables are defined 2260 (including the relative cells and the eventual user scalars):} \texttt{cs\_user\_physical\_properties} {\em {is not used.}} 2261 2262The user should ensure that the defined variation laws are valid for the whole range of 2263variables. Particular care should be taken with non-linear laws (for example, a 2264 $3^{rd}$ degree polynomial law giving negative values of density) 2265 2266{\em WARNING: In the electric module, all of the physical properties are considered as variables 2267 and are therefore stored using the \texttt{cs\_field} API. \texttt{cp0}, \texttt{viscls0} and \texttt{viscl0} 2268 are not used} 2269 2270For the Joule effect, the user is required to supply the physical properties in the 2271subroutine. Examples are given which are to be adapted by the user. If the temperature is 2272to be determined to calculate the physical properties, the solved variable, enthalpy must 2273 be deduced. The preferred temperature-enthalpy law should be defined 2274 (a general example is provided in (\texttt{cs\_user\_physical\_properties}), 2275 and can be used for the initialisation of the variables in 2276 (\texttt{cs\_user\_initialization})). 2277 For the electric arcs module, the physical properties are interpolated from the data file 2278 \texttt{dp\_ELE} supplied by the user. Modifications are generally not necessary. 2279 2280%================================== 2281\subsubsection{Boundary conditions} 2282%================================== 2283 2284For the electric module,each boundary face in \texttt{cs\_user\_boundary\_conditions} should be associated with a 2285 \texttt{izone} number \footnote{\texttt{izone} must be less than the maximum 2286 value allowed by the code, \texttt{nozzppm}. This is fixed at 2000 in \texttt 2287{ppvar} and cannot be modified.}(the color \texttt{icoul} for example) in 2288 order to group together all the boundary faces of the same type. In the 2289 \texttt{cs\_user\_boundary\_conditions} report, the main change from the users point of view concerns the 2290 specification of the boundary conditions of the potential, which isn't 2291 implied by default. The Dirichlet and Neumann conditions must be imposed 2292 explicitly using \texttt{icodcl} and \texttt{rcodcl} (as would be done for 2293 the classical scalar). 2294 2295Furthermore, if one wishes to slow down the power dissipation (Joule 2296effect module) or the current (electric arcs module) from the imposed values 2297\texttt{(puismp\index{puismp}} and \texttt{couimp\index{couimp}} respectively), 2298 they can be changed by the potential scalar as shown below: 2299 2300\begin{list}{-}{} 2301\item For the electric arcs, the imposed potential difference can be a fixed variable: 2302 for example, the cathode can be fixed at 0 and the potential at the anode 2303 contains the variable \texttt{dpot}. This variable is initialised in 2304in \texttt{cs\_user\_parameters.90} 2305 by an estimated potential difference. If \texttt{ielcor}=1 (see 2306 \texttt{cs\_user\_parameters.f90}), \texttt{dpot} is updated automatically during the 2307 calculation to obtain the required current. 2308\item For the Joule effect module, \texttt{dpot} is again used with the same 2309 signification as in the electric arcs module. If \texttt{dpot} is not wanted 2310 in the setting of the boundary conditions, the variable \texttt{coejou} can be 2311 used. \texttt{coejou} is the coefficient by which the potential difference is 2312 multiplied to obtain the desired power dissipation. By default this begins at 2313 1 and is updated automatically. If \texttt{ielcor}=1 (see \texttt 2314{cs\_user\_parameters.f90}), multiply the imposed potentials in 2315\texttt{cs\_user\_boundary\_conditions} by \texttt{coejou} at each time step to 2316achieve the desired power dissipation. 2317 \end{list} 2318 2319 {\em WARNING: In the case of alternating current, attention should be paid to the values of potential 2320 imposed at the limits: the variable named "real potential" represents an affective 2321 value if the current is in single phase, and a "real part" if not.} 2322\begin{list}{-}{} 2323\item For the Joule studies, a complex potential is sometimes needed 2324 (\texttt{ippmod(ieljou)}=2): this is the case in particular where the current 2325 has three phases. To have access to the phase of the potential, and not just to its 2326 amplitude, the two variables must be deleted: in \CS, there are two arrays 2327 specified for this role, the real part and the imaginary 2328 part of the potential. For use in the code, these variables are named 2329 ``real potential'' and ``imaginary potential''. For an alternative 2330 sinusoidal potential $Pp$, the maximum value is noted as $Pp_\text{max}$, 2331 the phase is noted as $\phi$, the real potential 2332 and the imaginary potential are respectively $Pp_\text{max}\,cos\phi$ and 2333$Pp_\text{max}\,sin\phi$. 2334\item For the Joule studies in which one does not have access to the phases, the real 2335 potential (imaginary part =0) will suffice (\texttt{ippmod(ieljou)=1}): this is 2336 obviously the case with 2337 continuous current, but also with single phase alternative current. In \CS 2338there is only 1 variable for the potential, called "real potential". Pay attention to 2339 the fact that in alternate current, the "real potential" represents a effective value 2340 of potential , $\frac{1}{\sqrt{2}}\,Pp_\text{max}$ (in continuous current there is no 2341 such ambiguity). 2342\end{list} 2343 2344\minititre{Additions for transformers} 2345 2346The following additional boundary conditions must be defined for tansformers: 2347\begin{list}{$\bullet$}{} 2348\item the intensity at each electrode 2349\item the voltage on each terminal of transformers. To achieve it, the intensity, 2350 the rvoltage at each termin, the Rvoltage, and the total intensity of the 2351transformer are calculated. 2352\end{list} 2353 2354Finally, a test is performed to check if the offset is zero or if a boundary 2355 face is in contact with the ground. 2356 2357%==================================== 2358\subsubsection {Initialisation of the variable options} 2359%================================== 2360\label{prg_cs_user_parameters}% 2361 2362The subroutine \texttt{cs\_user\_parameters} (in \texttt{cs\_user\_parameters.c}) 2363 is called at each time step. It allows: 2364\begin{list}{$\bullet$}{} 2365\item to give the coefficient of relaxation of the density \texttt{srrom}:\\ 2366$\rho^{n+1}=\texttt{srrom}*\rho^{n}+(1-\texttt{srrom})\rho^{n}$\\ 2367(for the electric arcs, the sub-relaxation is taken into account during the 2nd time 2368 step;) 2369 2370\item to indicate if the data will be fixed in the power dissipation or 2371 in the current, done in \texttt{ielcor}. 2372\item target either the current fixed as \texttt{couimp} (electric arcs module) 2373 or the power dissipation \texttt{puism} (Joule module effect). 2374\item to fix the initial value of potential difference \texttt{dpot}, 2375 the for the calculations with a single fixed parameter as \texttt{couimp} 2376 or \texttt{puism}. 2377\item to define type of scaling model for electric arcs \texttt{modrec}. If scaling 2378by a resetting plane is choosen then \texttt{idreca} defines the current density component 2379and \texttt{crit\_reca} the plane used for resetting of electromagnetic variables. 2380\end{list} 2381 2382%================================== 2383\subsection{\CS-\CS coupling} 2384%================================== 2385 2386\noindent 2387\textit{Subroutine called once during the calculation initialisation.} 2388 2389The user function \texttt{cs\_user\_saturne\_coupling} (in 2390 \texttt{cs\_user\_coupling.c} is used to couple \CS with itself. 2391 It is used for turbo-machine applications for instance, the first \CS managing 2392 the fluid around the rotor and the other the fluid around the stator. 2393In the case of a coupling between two \CS instances, first argument \texttt{saturne\_name} 2394 of the function '\texttt{cs\_sat\_coupling\_define}' is ignored. 2395 In case of multiple couplings, a coupling will be matched with available \CS 2396 instances based on that argument, which should match the directory name for the 2397 given coupled domain..\\ 2398The arguments of '\texttt{cs\_sat\_coupling\_define}' are: 2399\begin{list}{-}{} 2400\item \texttt{saturne\_name}: the matching \CS application name, 2401\item \texttt{volume\_sup\_criteria}: the cell selection criteria for support, 2402\item \texttt{boundary\_sup\_criteria}: the boundary face selection criteria for support (not functional), 2403\item \texttt{volume\_cpl\_criteria}: the cell selection criteria for coupled cells, 2404\item \texttt{boundary\_cpl\_criteria}: the boundary face selection criteria for coupled faces, 2405\item \texttt{verbosity}: the verbosity level. 2406\end{list} 2407 2408 2409%================================== 2410\subsection{Fluid-Structure external coupling} 2411%================================== 2412 2413\noindent 2414\textit{Subroutine called only once} 2415 2416The subroutine \texttt{usaste} belongs to the module dedicated to external 2417 Fluid-Structure coupling with \textit{Code\_Aster}. Here one defines the boundary 2418 faces coupled with \textit{Code\_Aster} and the fluid forces components which are 2419 given to structural calculation. When using external coupling with \textit{Code\_Aster}, 2420 structure numbers necessarily need to be negative; the references of coupled faces being 2421 i.e. -1, -2, \emph{etc.} 2422The subroutine performs the following operations: 2423\begin{list}{-}{} 2424 \item '\texttt{getfbr}' is called to get a list of elements matching a 2425geometrical criterion or reference number then a structure number (negative value) is associated 2426 to these elements. 2427 \item the value passed to \texttt{asddlf}, for user-chosen component, for every negative 2428 structure number, defines the movement imposed to the external structure. 2429\end{list} 2430% 2431\clearpage 2432% 2433%================================== 2434\subsection{ALE module} 2435%================================== 2436%================================== 2437\subsubsection{Initialisation of the options} 2438%================================== 2439This initialisation can be performed in the Graphical User Interface (GUI) 2440 or in the subroutines \texttt{usipph} and \texttt{usstr1}. Firstly, 2441 when the ``Mobile mesh'' is selected in GUI under the ``Calculation features'' 2442 heading, additional options are displayed. The user must choose the type of mesh 2443 viscosity and describe its spatial distribution, see \figurename~\ref{fig:Ini-ale}. 2444% 2445\begin{figure}[!ht] 2446\begin{center} 2447\includegraphics[width=12cm]{gui_ale_mei} 2448\caption{Thermophysical models - mobile mesh (ALE method)} 2449\label{fig:Ini-ale} 2450\end{center} 2451\end{figure} 2452% 2453The following paragraphs are relevant if the GUI is not used. 2454 2455\minititre{Subroutine \texttt{usipph}} 2456\noindent 2457\textit{Subroutine called at the beginning.} 2458This subroutine completes \texttt{cs\_user\_parameters.f90}. 2459 2460\texttt{usipph} allows setting options for the ALE module, and in 2461particular to activate the ALE module (\texttt{iale=1}). 2462 2463\minititre{Subroutine \texttt{usstr1}} 2464 2465This subroutine reads in \texttt{cs\_user\_fluid\_structure\_interaction.f90}. It allows to specify the following pieces of information for the structure module: 2466\begin{list}{-}{} 2467 \item the index of the structure, (\texttt{idfstr(ifac)} where \texttt{ifac} is the index of the face). Then the total number of structures \texttt{nbstru} is automatically computed by the code. Be careful, the value must belong to 1, ..., \texttt{nbstru}. 2468 \item the initial value of displacement, velocity and acceleration 2469 (\texttt{xstr0}, \texttt{xstreq} and \texttt{vstr0}). 2470\end{list} 2471 2472Below is a list of the different variables that might be modified: 2473 2474\begin{list}{$\bullet$}{} 2475 2476\item{\texttt{idfstr(ifac)}} \\ 2477{the index of the structure, (\texttt{idfstr(ifac)} where \texttt{ifac} is the index of the face), 0 if the face is not coupled to any structure.} 2478 2479\item{\texttt{xstr0(i,k)}} \\ 2480{initial position of a structure, where \texttt{i} is the dimension of space 2481and \texttt{k} the index of the structure} 2482 2483\item{\texttt{xstreq(i,k)}} \\ 2484{equilibrum position of a structure, where \texttt{i} is the dimension of space 2485and \texttt{k} the index of the structure} 2486 2487\item{\texttt{vstr0(i,k)}} \\ 2488{initial velicity of a structure, where \texttt{i} is the dimension of space 2489and \texttt{k} the index of the structure } 2490\end{list} 2491 2492%================================== 2493\subsubsection{Mesh velocity boundary conditions} 2494%================================== 2495These boundary conditions can be managed through the Graphical User Interface (GUI) 2496 or using the subroutine \texttt{usalcl} (called at each time step). With the GUI, 2497 when the item ``Mobile mesh'' is activated the item ``Fluid structure interaction'' 2498 appears under the heading ``Boundary conditions''. Two types of fluid-structure 2499 coupling are offered. The first one is internal, using a simplified structure 2500 model and the second is external with \textit{Code\_Aster}, see \figurename~ 2501 \ref{fig:CL-ale1} and \figurename~\ref{fig:CL-ale2}. 2502% 2503\begin{figure}[ht] 2504\begin{center} 2505\includegraphics[width=0.9\textwidth]{gui_ale_internal} 2506\caption{Boundary conditions - internal coupling} 2507\label{fig:CL-ale1} 2508\end{center} 2509\end{figure} 2510% 2511\begin{figure}[ht] 2512\begin{center} 2513\includegraphics[width=0.9\textwidth]{gui_ale_external} 2514\caption{Boundary conditions - external coupling} 2515\label{fig:CL-ale2} 2516\end{center} 2517\end{figure} 2518 2519\minititre{Subroutine \texttt{usalcl}} 2520When the GUI is not used, the use of \texttt{usalcl} is mandatory to run 2521a calculation using 2522the ale module just as it is in \texttt{cs\_user\_parameters.f90}. It is used the same way as 2523\texttt{cs\_user\_boundary\_conditions} in the framework of 2524standard calculations, that is to say a loop on the boundary faces 2525marked out by their colour (or more generally by a property of their 2526family), where the type of mesh velocity boundary condition is definied for 2527each variable. 2528 2529The main numerical variables are described below. 2530 2531\variab{ialtyb}{ialtyb(nfabor)}{ia}{In the ale module, the user 2532defines the mesh velocity from the colour of the boundary faces, or 2533more generally from their properties (colours, groups, ...), from the 2534boundary conditions defined in \texttt{cs\_user\_boundary\_conditions}, or even from their 2535coordinates. To do so, the array \texttt{ialtyb(nfabor)} gives for each face 2536\texttt{ifac} the mesh velocity boundary condition types marked out by the key 2537words \texttt{ivimpo\index{ivimpo}}, \texttt{igliss\index{igliss}}, 2538\texttt{ibfixe\index{ibfixe}} or \texttt{ifresf\index{ifresf}}}. 2539 2540\begin{list}{$\bullet$}{} 2541 2542\item If \texttt{ialtyb(ifac) = ivimpo}: imposed velocity. 2543 2544\begin{list}{$\rightarrow$}{} 2545\item In the cases where all the nodes of a face have a imposed displacement, 2546 it is not necessary to fill the tables with mesh velocity boundary conditions 2547 for this face, these will be erased. In the other case, 2548 the value of the Dirichlet must be given in \texttt{rcodcl(ifac,ivar,1)} for 2549 every value of \texttt{ivar} (\texttt{iuma}, \texttt{ivma} and \texttt{iwma}). 2550 The other boxes of \texttt{rcodcl} and \texttt{icodcl} are completed automatically. 2551 2552 The tangential mesh velocity is taken like a tape speed under the 2553 boundary conditions of wall for the fluid, except if wall fluid velocity 2554 was specified by the user in the interface or \texttt{cs\_user\_boundary\_conditions} (in which case 2555 it is this speed which is considered). 2556\end{list} 2557 2558 \item if \texttt{ialtyb(ifac) = ibfixe}: fixed wall 2559\begin{list}{$\rightarrow$}{} 2560 \item the velocity is null. 2561\end{list} 2562 2563 \item if \texttt{ialtyb(ifac) = igliss}: sliding wall 2564\begin{list}{$\rightarrow$}{} 2565\item symmetry boundary condition on the mesh velocity vector, which means a homogeneous Neumann on the tangential mesh velocity and a zero Dirichlet on the normal mesh velocity. 2566\end{list} 2567 2568 \item if \texttt{ialtyb(ifac) = ifresf}: free-surface 2569\begin{list}{$\rightarrow$}{} 2570\item an imposed mesh velocity such that the fluid mass flux is equal to the mesh displacement in order to mimic the free-surface automatically. Note that the boundary condition on the fluid velocity must be set separately (homogeneous Neumann condition for instance). 2571\end{list} 2572 2573\end{list} 2574 2575%================================== 2576\subsubsection{Modification of the mesh viscosity} 2577%================================== 2578 2579The user subroutine \texttt{cs\_user\_physical\_properties} can be used along the ALE (Arbitrary Lagrangian Eulerian 2580 Method) module, and allows modifying the mesh viscosity. 2581 It is called before the time loop, and before reading restart files 2582 (so the mesh is always in its initial position at this stage). 2583The user can modify mesh viscosity values to prevent cells and nodes from huge 2584 displacements in awkward areas, such as boundary layer for example. 2585 2586Note that for more complex settings, the mesh viscosity could be modified in 2587 \texttt{cs\_user\_initialization} or \texttt{cs\_user\_extra\_operations}. 2588 The matching field's name is \texttt{mesh\_viscosity}. 2589 2590%================================== 2591\subsubsection{Fluid - Structure internal coupling}\label{sec:ALE} 2592%================================== 2593 2594In the subroutine \texttt{cs\_user\_fluid\_structure\_interaction} the user provides the parameters of two other subroutines. 2595\texttt{usstr1} is called at the beginning of the calculation. It is used to define 2596 and initialise the internal structures where fluid-Structure coupling occurs. 2597For each boundary face \texttt{ifac}, \texttt{idfstr(ifac)} is the index of the structure 2598 the face belongs to (if \texttt{idfstr(ifac)} = 0, the face \texttt{ifac} doesn't belong 2599 to any structure). When using internal coupling, structure index necessarily must be 2600 strictly positive and smaller than the number of structures. The number of "internal" structures is automatically defined with the maximum 2601 value of the \texttt{idfstr} table, meaning that internal structure numbers must be defined 2602 sequentially with positive values, beginning with integer value '1'. 2603 2604For each internal structure the user can define: 2605\begin{list}{-}{} 2606 \item an initial velocity \texttt{vstr0} 2607 \item an initial displacement \texttt{xstr0} ({\em i.e.} \texttt{xstr0} is the value of the 2608 displacement \texttt{xstr} compared to the initial mesh at time t = 0) 2609 \item a displacement compared to equilibrium \texttt{xstreq} (i.e. \texttt{xstreq} 2610 is the initial displacement of the internal structure compared to its position at 2611 equilibrium; at each time step t and for a displacement \text{xstr}(t), the associated 2612 internal structure will undergo a force $-k*(\text{}t+XSTREQ)$ due to the spring). 2613\end{list} 2614\text{xstr0} and \text{vstr0} are initialised with the value 0. 2615When starting a calculation using ALE, or re-starting a calculation with ALE, based 2616 on a first calculation without ALE, an initial iteration 0 is automatically performed 2617 in order to take initial arrays \text{xstr0}, \text{vstr0} and \text{xstreq} into 2618 account. In any other case, add the following expression '\text{italin=1}' in subroutine 2619 \text{usipsu}, so that the code can deal with the arrays \text{xstr0}, \text{vstr0} and \text{xstreq}. 2620 2621When \texttt{ihistr} is set to 1, the code writes in the output the history of the 2622 displacement, of the structural velocity, of the structural acceleration and of the 2623 fluid force. The value of structural history output step is the same as the one for 2624 standard variables \text{nthist}. 2625 2626The second subroutine, \text{usstr2}, is called at each iteration. One defines in this 2627 subroutine structural parameters (considered as potentially time dependent): {\em i.e.}, 2628 mass m \text{xmstru}, friction coefficients c \text{xcstru}, and stiffness k \text{xkstru}. 2629 \text{forstr} array gives fluid stresses acting on each internal structure. Moreover it is also 2630 possible to take external forces (gravity for example ) into account. 2631\begin{list}{.}{} 2632 \item the \text{xstr} array indicates the displacement of the structure compared to its position in the initial mesh, 2633 \item the \text{xstr0} array gives the displacement of the structures in the initial mesh 2634 compared to structural equilibrium, 2635 \item the \text{vstr} array stands for structural velocity. 2636\end{list} 2637\text{xstr}, \text{xstr0} and \text{vstr} are \text{DATA} tables that can be used to 2638 define the Mass, Friction and Stiffness arays. These are not to be modified. 2639 2640The 3D structural equation that is solved is the following one: 2641\begin{equation}\label{eq:FluidStruct} 2642\displaystyle 2643\tens{m}.\partial_{tt}\vect{x}+\tens{c}.\partial_{t}\vect{x}+\tens{k}.\left(\vect{x}+\vect{x_0}\right)=\vect{f}, 2644\end{equation} 2645where $x$ stands for the structural displacement compared to initial mesh position 2646 \text{xstr}, $x_0$ represents 2647 the displacement of the structure in initial mesh compared to equilibrium. 2648Note that $\tens{m}$,$\tens{c}$, and $\tens{k}$ are 3\text{x}3 matrices. 2649Equation \eqref{eq:FluidStruct} is solved using a Newmark HHT algorithm. 2650Note that the time step used to solve this equation, \text{dtstr}, can be 2651 different from the one of fluid calculations. The user is free to define \text{dtstr} 2652 array. At the beginning of the calculation \text{dtstr} is initialised to the value of 2653 \text{dtcel} (fluid time step). 2654 2655%================================== 2656\subsection{Management of the structure property} 2657%================================== 2658 2659The use of \texttt{usstr2} is mandatory to run a calculation using the ALE 2660 module with a structure module. It is called at each time step. 2661 2662For each structure, the system that will be solved is: 2663 2664\begin{equation} 2665M.x^{''}+C.x^{'}+K.(x-x_{0}) = 0 2666\end{equation} 2667 2668where 2669 2670\begin{list}{-}{} 2671 \item $M$ is the mass structure (\texttt{xmstru}). 2672 \item $C$ is the damping coefficient of the structure (\texttt{xcstru}). 2673 \item $K$ is the spring constant or force constant of the structure (\texttt{xkstru}). 2674 \item $x_{0}$ is the initial position. 2675\end{list} 2676 2677Below is a list of the different variables that might be modified: 2678 2679\begin{list}{$\bullet$}{} 2680 2681\item{\texttt{xmstru(i,j,k})} \\ 2682{mass matrix of the structure, where \texttt{i},\texttt{j} is 2683the array of mass structure and \texttt{k} the index of the structure.} 2684 2685\item{\texttt{xcstru(i,j,k})}\\ 2686{damping matrix coefficient of the structure, where \texttt{i},\texttt{j} is the array of 2687damping coefficient and \texttt{k} the index of the structure.} 2688 2689\item{\texttt{xkstru(i,j,k)}}\\ 2690{spring matrix constant of the structure, where \texttt{i},\texttt{j} is the array of spring 2691constant and \texttt{k} the index of the structure.} 2692 2693\item{\texttt{forstr(i,k)}}\\ 2694{force vector of the structure, where \texttt{i} is the force vector and 2695\texttt{k} the index of the structure.} 2696\end{list} 2697 2698%=========================================================== 2699\subsection{Management of the atmospheric module} 2700%=========================================================== 2701 2702This section describes how to set a calculation using the atmospheric module 2703of \CS. Each paragraph describes a step of the data setting process. 2704 2705\subsubsection{Directory structure}\label{sec:atmo_struct} 2706% 2707The flowchart (\figurename~\ref{fig:organisation}) recalls the directory 2708structure of a study generated by \CS (see also \ref{sec:prg_architecture}). 2709When using the atmospheric module, the structure is identical but a file called 2710\texttt{meteo} may be added to the data settings in order to provide vertical 2711profiles of the main variables. This file should be put in the \texttt{DATA} 2712directory. For more details about the \texttt{meteo} file, see \S~ 2713~\ref{sec:atmo_BCs}). 2714% 2715\begin{figure}[ht] 2716 \centerline{\includegraphics[width=0.6\textwidth]{gui_atmospheric_user_s_guide_v91.png}} 2717 \caption{Organization of a study (specific files of atmospheric version in bold type)} 2718 \label{fig:organisation} 2719\end{figure} 2720% 2721\subsubsection{The atmospheric mesh features} 2722% 2723An atmospheric mesh has the following specific features: 2724% 2725\begin{itemize} 2726\item The boundary located at the top of the domain should be a plane. 2727So, horizontal wind speed at a given altitude can be prescribed at the top 2728face as an inlet boundary. 2729\item Cells may have very different sizes, from very small (near ground or 2730buildings) to very large (near the top of domain or far from zone of interest). 2731\item Vertical resolution: from tiny cells (e.g. $\Delta $\upshape z = 1 m) near 2732the ground to a few hundreds of meters at the top. 2733\item Horizontal resolution: from a few meters to hundreds of meters. 2734\item The length ratio between two adjacent cells (in each direction) should 2735preferably be between $0.7$ and $1.3$. 2736\item The z axis represents the vertical axis. 2737\end{itemize} 2738% 2739A topography map can be used to generate a mesh. In this case, the preprocessor 2740 mode is particularly useful to check the quality of the mesh (run type Mesh 2741quality criteria). 2742% 2743\subsubsection{Atmospheric flow model and steady/unsteady algorithm} 2744% 2745The Graphical User Interface (GUI) may be used to enable the atmospheric flow 2746module and set up the following calculation parameters in the 2747\texttt{Thermophysical models}-\texttt{Calculation features} page 2748(see \figurename~\ref{fig:steady}): 2749% 2750\subsubsubsection{The atmospheric flow model} 2751% 2752The user can choose one of the following atmospheric flow models: 2753% 2754\begin{itemize} 2755\item \textbf{Constant density}: To simulate neutral atmosphere. 2756\item \textbf{Dry atmosphere}: To simulate dry, thermally-stratified 2757atmospheric flows (enables \texttt{Potential temperature} as thermal model). 2758\item \textbf{Humid atmosphere}: To simulate thermally stratified atmospheric 2759flows (air-water mixture) with phase changes (enables \texttt{Liquid potential 2760temperature} as thermal model). The model is described in 2761Bouzereau \cite{bouzereau}. 2762\end{itemize} 2763% 2764\subsubsubsection{The time algorithm} 2765% 2766\begin{itemize} 2767\item Steady flow algorithm: is the one usually set. It sets a time step 2768variable in space and time. It has to be selected if constant boundary 2769conditions are used. 2770\item Unsteady flow algorithm has to be selected for time varying boundary 2771conditions (the time step can then be variable in time or constant). 2772\end{itemize} 2773% 2774\begin{figure}[ht] 2775\centerline{\includegraphics[width=0.7\textwidth]{gui_atmospheric_user_s_guide_v92.png}} 2776\caption{Selection of atmospheric model} 2777\label{fig:steady} 2778\end{figure} 2779% 2780\begin{figure}[ht] 2781\centerline{\includegraphics[width=0.9\textwidth]{gui_atmospheric_user_s_guide_v93.png}} 2782\caption{Selection of steady/unsteady flow algorithm} 2783\label{fig:global} 2784\end{figure} 2785% 2786Table \tablename~\ref{tab:param_list} can help to choose the right parameters 2787depending on the type of atmospheric flow. 2788% 2789\subsubsubsection{Warnings} 2790The following points have to be considered when setting the parameters 2791described above: 2792\begin{itemize} 2793\item The potential temperature thermal model and the liquid potential 2794temperature one (see the paragraph ``Atmospheric main variables'' for the 2795definition) requires that the vertical component of the gravity is set to 2796$g_z=-9.81 m.s^{-2}$ ($g_x=g_y=0 m.s^{-2}$), 2797otherwise pressure and density won't be correctly computed. 2798\item As well, the use of scalar with drift for atmospheric dispersion requires 2799the gravity to be set to $g_z=-9.81$ ($g_x=g_y=0 m.s^{-2}$), even if the density 2800is constant. 2801\end{itemize} 2802% 2803\subsubsection{Physical properties} 2804% 2805The specific heat value has to be set to the atmospheric value 2806$C_{p}=1005 J/kg/K$. 2807% 2808\begin{table}[ht]\label{tab:param_list} 2809\begin{center} 2810\begin{tabular}{|p{80pt}|p{70pt}|p{70pt}|p{80pt}|p{100pt}|} 2811\hline 2812\textbf{Parameters} & \textbf{Constant density} & \textbf{Dry atmosphere} & \textbf{Humid atmosphere} & \textbf{Explanation} \\ 2813\hline 2814pressure boundary condition & Neumann first order & Extrapolation & Extrapolation & 2815In case of \textbf{Extrapolation}, the pressure gradient is assumed (and set) constant, whereas in case of \textbf{Neumann first order}, 2816the pressure gradient is assumed (and set) to zero. \\ 2817\hline 2818Improved pressure interpolation in stratified flows & no & yes & yes & If yes, exact balance between the hydrostatic part of the pressure gradient and the gravity term $\rho g$ is numerically ensured. \\ 2819\hline 2820Gravity (gravity is assumed aligned with the z-axis) & $g_z=0$ or $g_z=-9.81 m.s^{-2}$ (the latter is useful for scalar with drift) & $g_z=-9.81 m.s^{-2}$ & $g_z=-9.81 m.s^{-2}$ & \\ 2821\hline 2822Thermal variable & no & potential temperature & liquid potential temperature & \\ 2823\hline 2824Others variables & no & no & total water content, droplets number & \\ 2825\hline 2826\end{tabular}\label{tab1} 2827\caption[List of parameters]{List of parameters} 2828\end{center} 2829\end{table} 2830% 2831\subsubsection{Boundary and initial conditions}\label{sec:atmo_BCs} 2832% 2833The \texttt{meteo} file can be used to define initial conditions for the 2834different fields and to set up the inlet boundary conditions. For the velocity 2835field, \CS can automatically detect if the boundary is an inlet boundary or an 2836outflow boundary, according to the wind speed components given in the 2837\texttt{meteo} file with respect to the boundary face orientation. This is often 2838used for the lateral boundaries of the atmospheric domain, especially if the 2839profile is evolving in time. In the case of inlet flow, the data given in the 2840\texttt{meteo} file will be used as the input data (Dirichlet boundary condition) 2841for velocity, temperature, humidity and turbulent variables. In the case of 2842outflow, a Neumann boundary condition is automatically imposed (except for the 2843pressure). The unit of temperature in the \texttt{meteo} file is the degree 2844Celsius whereas the unit in the GUI is the kelvin. 2845 2846To be taken into account, the \texttt{meteo} file has to be selected in the GUI 2847(\texttt{Atmospheric flows} page, see \figurename~\ref{fig:meteo}) and the check 2848box on the side ticked. This file gives the profiles of prognostic atmospheric 2849variables containing one or a list of time stamps. The file has to be put in the 2850\texttt{DATA} directory. 2851% 2852\begin{figure}[htbp] 2853\centerline{\includegraphics[width=0.9\textwidth]{gui_atmo_read.png}} 2854\caption{Selection of the \texttt{meteo} file} 2855\label{fig:meteo} 2856\end{figure} 2857% 2858An example of file \texttt{meteo} is given in the directory 2859\texttt{DATA/REFERENCE/}. The file format has to be strictly respected. 2860The horizontal coordinates are not used at the present time (except when 2861boundary conditions are based on several meteorological vertical profiles) 2862and the vertical profiles are defined with the altitude above sea level. The 2863highest altitude of the profile should be above the top of the simulation domain 2864and the lowest altitude of the profile should be below or equal to the lowest 2865level of the simulation domain. The line at the end of the \texttt{meteo} file 2866should not be empty. 2867 2868If the boundary conditions are variable in time, the vertical profiles for 2869the different time stamps have to be written sequentially in the \texttt{meteo} 2870file. 2871 2872You can also set the profiles of atmospheric variables directly in the GUI. 2873The following boundary conditions can be selected in the GUI: 2874% 2875\begin{itemize} 2876\item Inlet/Outlet is automatically calculated for lateral boundaries 2877(e.g. North, West\textellipsis ) of the computational domain 2878(see \figurename~\ref{fig:inlet}). 2879\item Inlet for the top of the domain (see \figurename~\ref{fig:top}). 2880\item Rough wall for building walls (see \figurename~\ref{fig:walls}) or for 2881the ground (see \figurename~\ref{fig:ground}). 2882The user has to enter the roughness length. In case of variable roughness 2883length, the user has to provide the land use data and the association 2884between the roughness length values and land use categories. 2885\end{itemize} 2886% 2887\begin{figure}[htbp] 2888\centerline{\includegraphics[width=0.9\textwidth]{gui_atmospheric_user_s_guide_v95.png}} 2889\caption{Selection of automatic inlet/ outlet for boundary conditions } 2890\label{fig:inlet} 2891\end{figure} 2892% 2893\begin{figure}[htbp] 2894\centerline{\includegraphics[width=0.9\textwidth]{gui_atmospheric_user_s_guide_v96.png}} 2895\caption{Selection of the boundary condition for the top of the domain } 2896\label{fig:top} 2897\end{figure} 2898% 2899\begin{figure}[htbp] 2900\centerline{\includegraphics[width=0.9\textwidth]{gui_atmospheric_user_s_guide_v97.png}} 2901\caption{Selection of the boundary condition for building walls} 2902\label{fig:walls} 2903\end{figure} 2904% 2905\begin{figure}[htbp] 2906\centerline{\includegraphics[width=0.9\textwidth]{gui_atmospheric_user_s_guide_v98.png}} 2907\caption{Selection of the boundary condition for the ground} 2908\label{fig:ground} 2909\end{figure} 2910 2911\paragraph{Remark:} If a meteorological file is given, it is used by default to 2912initialize the variables. If a meteorological file is not given, the user can 2913use the standard \CS initial and boundary conditions set up but has to be aware 2914that even small inconsistencies can create very large buoyancy forces and 2915spurious circulations. 2916 2917\subsubsubsection{Boundary conditions based on several meteorological vertical profiles} 2918 2919In some cases, especially when outputs of a mesoscale model are used, you 2920need to build input boundary conditions from several meteorological vertical 2921wind profiles. Cressman interpolation is then used to create the boundary 2922conditions. The following files need to be put in the \texttt{DATA} directory: 2923\begin{itemize} 2924\item All \texttt{meteo} files giving the different vertical profiles of 2925prognostic variables (wind, temperature, turbulent kinetic energy and 2926dissipation). 2927\item A file called \texttt{imbrication\_files\_list.txt} which is a list 2928of the \texttt{meteo} files used. 2929\item A separate \texttt{meteo} file which is used for the initial conditions 2930and to impose inlet boundary conditions for the variables for which Cressman 2931interpolation is not used (for example: temperature, turbulent kinetic energy). 2932This file must follow the rules indicated previously. 2933\end{itemize} 2934The following files should be put in the SRC directory: 2935\begin{itemize} 2936\item The user source file \texttt{cs\_user\_parameters.f90}. In this file, set 2937the \texttt{cressman\_} flag of each variable, for which the Cressman 2938interpolation should be enabled, to \texttt{.true.}. 2939\end{itemize} 2940% 2941\subsubsection{User subroutines} 2942% 2943The user subroutines are used when the graphical user interface is not 2944sufficient to set up the calculation. We give some examples of user file for 2945atmospheric application: 2946\begin{itemize} 2947\item \texttt{cs\_user\_source\_terms.f90}: to add a source term in the 2948prognostic equations for forest canopy modelling, wind turbine wake modelling... 2949See the associated \doxygenfile{cs_user_source_terms.html}{\texttt{doxygen} documentation 2950for examples of use of \texttt{cs\_user\_source\_terms.f90}}. 2951\item \texttt{cs\_user\_parameters.f90}: to activate the Cressman interpolation. 2952For example, it is used to impose inhomogeneous boundary conditions. See the associated 2953\doxygenfile{f_parameters.html}{\texttt{doxygen} documentation for examples of use of 2954\texttt{cs\_user\_parameters.f90}}. 2955\item \texttt{cs\_user\_extra\_operations-extract.f90}: to generate vertical 2956profiles for post processing. See the associated 2957\doxygenfile{cs_user_extra_operations_examples.html}{\texttt{doxygen} documentation 2958for examples of use of \texttt{cs\_user\_extra\_operations.f90}}. 2959\item \texttt{cs\_user\_boundary\_conditions-atmospheric.f90}: show how to set 2960up the boundary conditions and to put a heterogeneous roughness length... 2961See the associated \doxygenfile{cs_user_boundary_conditions_examples.html}{\texttt{doxygen} 2962documentation for examples of use of \texttt{cs\_user\_boundary\_conditions.f90}}. 2963\end{itemize} 2964% 2965\paragraph{Remark:} 2966If the computation is set without the GUI, other user subroutines such as the 2967following have to be used: 2968\begin{itemize} 2969\item \texttt{cs\_user\_initialization-atmospheric.f90}: allows to initialize 2970or modify (in case of a restarted calculation) the calculation variables and 2971the values of the time step. See the associated 2972\doxygenfile{user_initialization.html}{\texttt{doxygen} documentation for 2973examples of use of \texttt{cs\_user\_initialization.f90}}. 2974 2975\item \texttt{cs\_user\_boundary\_conditions-atmospheric.f90}: allows to define 2976all the boundary conditions. For each type of boundary condition, faces should 2977be grouped as physical zones characterized by an arbitrary number \texttt{izone} 2978chosen by the user. If a boundary condition is retrieved from a meteorological 2979profile, the variable \texttt{iprofm(izone)} of the zone has to be set to 1. 2980The vertical profiles of atmospheric variables can be described in this file. 2981\end{itemize} 2982% 2983Examples are available in the directory \texttt{SRC/EXAMPLE}. 2984% 2985\subsubsection{Physical models} 2986% 2987\subsubsubsection{Atmospheric dispersion of pollutants} 2988% 2989To simulate the atmospheric dispersion of pollutant, one first need to define 2990the source(s) term(s). That is to say the location i.e. the list of cells or 2991boundary faces, the total air flow, the emitted mass fraction of pollutant, 2992the emission temperature and the speed with the associated turbulent parameters. 2993The mass fraction of pollutant is simulated through a user added scalar that 2994could be a `scalar with drift' if wanted (aerosols for example). 2995 2996The simulations can be done using 3 different methods: 2997\begin{enumerate} 2998\item Using a mass source term, that is added in the Navier-Stokes 2999equations using the \texttt{cs\_user\_mass\_source\_terms.f90} user subroutine. 3000 3001\item Prescribing a boundary condition code ``total imposed mass flux`` for 3002some boundary faces using the \texttt{cs\_user\_boundary\_conditions.f90} user 3003subroutine. 3004 3005\item Using a scalar source term. In this case, the air inflow is not taken 3006into account. The user has to add an explicit part to the equations 3007for the scalar through the \texttt{cs\_user\_source\_terms.f90} file. This is 3008done by selecting the cells and adding the source term \texttt{crvexp (cells)} 3009which equals to the air flux multiplied by the mass fraction, while the 3010implicit part \texttt{crvimp} is set to zero. 3011\end{enumerate} 3012 3013The first method is recommended, but one must take care that each source 3014influences the dispersion of the others, which is physically realistic. So 3015if the impact of several sources has to be analyzed independently it has first 3016to be verified that these influences are negligible or as many simulations 3017as there are sources have to be run. 3018 3019With the second method, the same problem of sources interactions appears, and 3020moreover standard Dirichlet conditions should not be used (use 3021\texttt{itypfb=i\_convective\_inlet} and \texttt{icodcl=13} instead) as 3022the exact emission rate cannot be prescribed because the diffusive part 3023(usually negligible) cannot be quantified. Additionally, it requires that 3024the boundary faces of the emission are explicitly represented in the mesh. 3025 3026Finally the third method does not take into account the jet effect of the 3027emission and so must be used only if it is sure that the emission does not 3028modify the flow. 3029 3030Whatever solution is chosen, the mass conservation should be verified by using 3031for example the 3032\texttt{cs\_user\_extra\_operations-scalar\_balance\_by\_zone.f90} file. 3033 3034\subsubsubsection{Soil/atmosphere interaction model} 3035 3036This model is based on the force restore model (Deardorff \cite{deardorff}). 3037It takes into account heat and humidity exchanges between the ground and the 3038atmosphere at daily scale and the time evolution of ground surface temperature 3039and humidity. Surface temperature is calculated with a prognostic equation 3040whereas a 2-layers model is used to compute surface humidity. 3041 3042The parameter \texttt{iatsoil} in the file \texttt{atini0.f90} needs to be equal to one to 3043activate the model. Then, the source file \texttt{solvar.f90} is used. 3044 3045Three variables need to be initialized in the file \texttt{atini0.f90}: deep soil 3046temperature, surface temperature and humidity. 3047 3048The user needs to give the values of the model constants in the file 3049\texttt{solcat.f90}: roughness length, albedo, emissivity... 3050 3051In case of a 3D simulation domain, land use data has to be provided for the domain. 3052Values of model constants for the land use categories have also to be 3053provided. 3054 3055\subsubsubsection{Radiative model (1D)} 3056 3057The 1D-radiative model calculates the radiative exchange between different 3058atmospheric layers and the surface radiative fluxes. 3059 3060The radiative exchange is computed separately for two wave lengths intervals 3061 3062\begin{itemize} 3063\item Calculation in the infrared spectral domain (file \texttt{rayir.f90}) 3064\item Calculation in the spectral range of solar radiation (file 3065\texttt{rayso.f90}) 3066\end{itemize} 3067This 1D-radiative model is needed if the soil/atmosphere interaction model 3068is activated. 3069 3070This model is activated if the parameter \texttt{iatra1} is equal to one in the 3071file \texttt{cs\_users\_parameters.f90}. 3072 3073\subsubsection{Atmospheric main variables} 3074 3075For more details on the topic of atmospheric boundary layers, see Stull 3076\cite{stull}. 3077% 3078\begin{itemize} 3079\item Definition of the potential temperature: 3080\[ 3081\theta =T\left(\frac{P}{P_{r}}\right)^{-\frac{R_{d}}{C_{p}}} 3082\] 3083\item Definition of liquid potential temperature: 3084\[ 3085\theta_{l} = \theta \left( 1-\frac{L}{C_{p}T} q_{l} \right) 3086\] 3087\item Definition of virtual temperature: 3088\[ 3089T_{v} = \left(1+0.61q\right)T 3090\] 3091\item Gas law: 3092\[ 3093P = \rho \frac{R}{M_{d}}\left(1+0,61q\right)T 3094\] 3095with $R=R_{d} M_{d}$. 3096\item Hydrostatic state: 3097\[ 3098\frac{\partial P}{\partial z} = -\rho g 3099\] 3100\end{itemize} 3101% 3102\begin{table}[htbp] 3103\begin{center} 3104\begin{tabular}{|l|l|l|l|} 3105\hline 3106\textbf{Constant name} & \textbf{Symbol} & \textbf{Values} & \textbf{Unit} \\ 3107\hline 3108Gravity acceleration at sea level & $g$ & $9.81$ & $m.s^{-2}$ \\ 3109\hline 3110Effective Molecular Mass for dry air & $M_{d}$ & $28.97$ & $kg.kmol^{-1}$ \\ 3111\hline 3112Standard reference pressure & $P_{r}$ & $10^{5}$ & $Pa$ \\ 3113\hline 3114Universal gas constant & $R$ & $8.3143$ & $J.K^{-1}.mol$ \\ 3115\hline 3116Gas constant for dry air & $R_{d}$ & $287$ & $J.kg^{-1}.K^{-1}$ \\ 3117\hline 3118\end{tabular}\label{tab2} 3119\end{center} 3120\caption{Constant name} 3121\end{table} 3122 3123\begin{table}[htbp] 3124\begin{center} 3125\begin{tabular}{|l|p{113pt}|} 3126\hline 3127\textbf{Variable name} &\textbf{Symbol} \\ 3128\hline 3129Specific heat capacity of dry air & $C_{p}$ \\ 3130\hline 3131Atmospheric pressure & $P$ \\ 3132\hline 3133Specific humidity & $q$ \\ 3134\hline 3135Specific content for liquid water & $q_{l}$ \\ 3136\hline 3137Temperature & $T$ \\ 3138\hline 3139Virtual temperature & $T_{v}$ \\ 3140\hline 3141Potential temperature & $\theta$ \\ 3142\hline 3143Liquid potential temperature & $\theta_{l}$ \\ 3144\hline 3145Latent heat of vaporization & $L$ \\ 3146\hline 3147Density & $\rho $ \\ 3148\hline 3149Altitude & $z$ \\ 3150\hline 3151\end{tabular}\label{tab3} 3152\caption{Variable name} 3153\end{center} 3154\end{table} 3155% 3156\subsubsection{Recommendations}\label{subsubsec:recommendations} 3157% 3158This part is a list of recommendations for atmospheric numerical simulations. 3159% 3160\begin{itemize} 3161\item Enough probes at different vertical levels in the domain should be used 3162to check the convergence of the calculation. 3163\item An inflow boundary condition at the top level of the domain should be set 3164(symmetry and automatic inlet/outlet are not appropriate). 3165\item A Courant number too small or too big has to be avoided (see \CS Best 3166Practice Guidelines). That is the reason why the option \texttt{variable time 3167step in space and in time} is recommended for steady simulations when there are 3168large differences of cell size inside the domain (which is generally the case 3169for atmospheric simulations). With this option, it can be necessary to change 3170the reference time step and the time step maximal increase (by default, the time 3171step increase rate is $10\%$). 3172\end{itemize} 3173% 3174In some cases, results can be improved with the following modifications: 3175% 3176\begin{itemize} 3177\item In some case, the turbulent eddy viscosity can drop to unrealistically low 3178values (especially with $k-\varepsilon$ model in stable atmospheric condition). 3179In those cases, it is suggested to put an artificial molecular viscosity around 3180$0.1 m^{2}.s^{-1}$. 3181\item If the main direction of wind is parallel to the boundary of your computing 3182domain, try to set symmetry boundary conditions for the lateral boundaries to 3183avoid inflow and outflow on the same boundary zone (side of your domain). 3184Another possibility is to use a cylindrical mesh. 3185\item To avoid inflow and outflow on the same boundary zone (side of your domain), 3186avoid the case of vertical profile in the input data \texttt{meteo} file with 3187changes of the sign of velocity of wind ($V_x$ or/and $V_y$). 3188\end{itemize} 3189 3190%================================== 3191%\subsection{Cooling tower modelling} 3192%================================== 3193 3194%================================== 3195%\subsubsection{Parameters} 3196%================================== 3197 3198%\noindent 3199%\textit{Subroutine called only during calculation initialisation? OR AT EACH ITERATION?.} 3200 3201%The subroutine \texttt{uscti1} contains calculation parameters such as: 3202%\begin{list}{-}{} 3203% \item temperature parameters, 3204% \item the number of exchange zones at various locations, 3205% \item the air properties. 3206%\end{list} 3207 3208%================================== 3209%\subsubsection{Initialisation of the variables} 3210%================================== 3211 3212%With the additional variables introduced with the air-cooling module, the 3213%following quantities may be initialized by \texttt{cs\_user\_initialization} in 3214%user selected zones in addition to the usual quantities: 3215%\begin{list}{-}{} 3216% \item the air temperature, variable of index \texttt{isca(ihumid)}, 3217% \item the air humidity, variable of index \texttt{isca(itemp4)}, 3218%\end{list} 3219%where \texttt{iel} can be an element found in a list returned by the routine 3220% '\texttt{getcel}'. 3221 3222%================================== 3223%\subsubsection{Definition of the exchange zones} 3224%================================== 3225 3226%The subroutine \texttt{usctdz} is used to define the exchange zones of a cooling 3227% tower. The user provides the following parameters: 3228%\begin{list}{-}{} 3229% \item \texttt{imzech}: its value is related to the model used: 3230% \begin{list}{$\rightarrow$}{} 3231% \item 0: no model is used, 3232% \item 1: Merkel model is used, 3233% \item 2: Poppe model is used, 3234% \end{list} 3235% \item 10 exchange zone parameters. 3236%\end{list} 3237%These arguments are passed to the subroutine '\texttt{defct}' along with a 3238%geometrical selection criterion. 3239 3240%================================== 3241\subsection{Turbomachinery computations} 3242%================================== 3243 3244\subsubsection{Introduction} 3245 3246Two classical models are available in \CS for rotor/stator 3247interactions modelling in turbomachinery computations: the steady 3248approach which is based on the so-called \emph{Frozen Rotor} modelling 3249and the \emph{transient rotor/stator} approach which is based on a 3250sliding mesh technique. 3251 3252\underline{\emph{Warning:}} This section describes these functionalities based on 3253a single \CS computation. An alternative rotor/stator coupling based 3254on coupling of boundary conditions is also possible (and only briefly 3255described in this section) but it is not recommended. 3256 3257\subsubsection{Meshing reccomendations}\label{tbm_mesh} 3258 3259\paragraph{Periodicity}\label{tbm_mesh_perio} 3260 3261The rotational periodicity treatment is possible only in \emph{Frozen 3262Rotor}. However, the interface plane between rotor and stator 3263must match in the azimutal $\theta$ direction: 3264$$\theta_{\text{min}}^{\text{rotor}}(z)=\theta_{\text{min}}^{\text{stator}}(z),\quad\theta_{\text{max}}^{\text{rotor}}(z)=\theta_{\text{max}}^{\text{stator}}(z)$$ 3265for all $z$ through the rotation axis direction. 3266 3267\paragraph{Rotor/stator interface}\label{tbm_mesh_interface} 3268\begin{itemize} 3269\item \emph{Unsteady rotor/stator}: in the input mesh(es), the 3270 interface between rotor and stator domains has to be composed of 3271 \underline{boundary faces}. Then the interface boundary faces are joined 3272 during the computation and become internal faces, as is usual for 3273 mesh joining in the preprocessing stage. A simple way to ensure 3274 joining is not done prematurely is to provide 3275 \underline{separated meshes} for each rotor or stator domain. 3276\item \emph{Frozen Rotor}: the interface can be composed of boundary 3277 faces (in which case the interface boundary faces are joined at 3278 the beginning of the computation) or of internal faces. 3279\end{itemize} 3280 3281\paragraph{Meshing of the interface region}\label{tbm_user_bpg_mesh} 3282 3283As mentioned above, when a rotor/stator interface boundary exists (in 3284particular for the \emph{unsteady rotor/stator} model), boundary faces 3285are joined by the solver during the computation, based on the current 3286rotor position. It is thus important to be aware that the success of 3287a joining operation is strongly dependant on the 3288\underline{quality of the mesh at the interface}. More precisely, 3289the refinement must be as similar as possible at both sides of the 3290interface. Moreover, it is reminded that the tolerance parameter of 3291a joining is a fraction of the shortest edge linked with a vertex of 3292a joined face. Consequently, cells with high aspect ratios where the 3293refinement in the azimutal $\theta$ direction is much coarser than 3294those in one of the two others can also lead to a joining failure. 3295In particular, the user should be careful to avoid elongated 3296viscous layer type cells in curved areas such as a rotor-stator interface. 3297 3298If the meshes at both sides of the interface are very different 3299such that the joining fails, advanced joining parameters are 3300available. However, modifying the mesh is more likely to 3301succeed. The introduction of a somekind of buffer cells layer on 3302both sides of the interface should be very valuable. Ideally, each 3303of the two layers should have the same refinement and a constant 3304azimutal step (this latter recommandation is relevant only for 3305\emph{unsteady rotor/stator} model). 3306 3307\paragraph{Alternative rotor/stator coupling}\label{tbm_user_bpg_cpl} 3308 3309If the meshes at both sides of the interface are very different and 3310can not be modified, a fallback solution is to use the rotor/stator model 3311based on the boundary conditions coupling. 3312 3313\underline{\emph{Warning:}} Contrarily to the mesh joining approach, the 3314boundary conditions coupling approach is not fully conservative. 3315 3316\subsubsection{Turbomachinery dedicated postprocessing functions}\label{tbm_user_post} 3317 3318Useful postprocessing functions relative to the machinery 3319characteristics are available: postprocessing of the couple on the 3320rotor walls and postprocessing of the head generated by the machinery. 3321 3322\subsubsection{Data setting, keywords and examples}\label{tbm_date_setting} 3323Data setting, keywords and examples for turbomachinery computations 3324(mesh joining or boundary conditions coupling), are provided in 3325\doxygenfile{turbomachinery.html}{the dedicated \texttt{doxygen} 3326documentation}. 3327 3328%================================== 3329\subsection{Cavitation module} 3330%================================== 3331 3332The cavitation module is based on an homogeneous mixture model. The 3333physical properties (density and dynamic viscosity) of the mixture 3334depends on a resolved void fraction and constant reference properties 3335of the liquid phase and the gas phase. 3336 3337For a description of the user management of the cavitation module, 3338please refer to \doxygenfile{cavit.html}{the dedicated 3339 \texttt{doxygen} documentation}. 3340
