\documentstyle[11pt,a4]{article}
\pagecolor{\rgb{255}{255}{204}}
\setlength{\topmargin}{0cm}
\setlength{\oddsidemargin}{0cm}
\setlength{\evensidemargin}{0cm}
\setlength{\textheight}{22.5cm}
\setlength{\textwidth}{16cm}
\footskip=27pt
\def\baselinestretch{0.950}
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\renewcommand{\thesubsection}{{\em \thesection.\arabic{subsection}}}

\newtheorem{lemma}{\sc \bf Lemma}[section]
\newtheorem{propos}{\sc \bf Proposition}[section]
\newtheorem{theor}{\sc \bf Theorem}[section]
\newtheorem{corr}{\sc \bf Corollary}[section]
\newtheorem{remark}{\sc \bf Remark}[section]
\newtheorem{definition}{\sc \bf Definition}[section]

\newenvironment{lem}{\begin{lemma} \hspace{1mm}}{\end{lemma}}
\newenvironment{prop}{\begin{propos} \hspace{1mm}}{\end{propos}}
\newenvironment{theo}{\begin{theor} \hspace{1mm}}{\end{theor}}
\newenvironment{cor}{\begin{corr} \hspace{1mm}}{\end{corr}}
\newenvironment{rem}{\begin{remark} \em \hspace{1mm}}{\end{remark}}
\newenvironment{defin}{\begin{definition} \r \vs \em }{\end{definition}}


\newcommand{\vb}{\mbox{${\bf v_t\mit}$}}

\newcommand{\divv}{\mbox{${\bf\nabla\mit}\cdot{\bf v_t\mit}$}}
\newcommand{\nb}{\mbox{${\bf\nabla\mit}$}}
\newcommand{\ddt} [1] {{\partial {#1}\over{\partial {t}}}}

%
%  New commands
%
\def\bkR{{\rm I\kern-.17em R}}
%
%  Begin document
%
\pagestyle{empty}
\begin{document}
\thispagestyle{empty}
%
%  Title
%
\title{\bf Subgrid-Scale  Modelling in Acoustic Finite-Difference Time-Domain Simulation }
\author{
     Dick Botteldooren, John De Poorter \\ \ \\
         {\small \it  Department of Information Technology, University of Gent, St.~Pietersnieuwstraat 41, B-9000 Gent, Belgium} }
\date{}
\maketitle
%
%  Abstract
%
\begin{abstract}
\noindent In principle, acoustic problems can be addressed by solving numerically the nonlinear equations that describe the dynamic behavior of a gas, including damping effects. In practise such an approach rarely leads to useful simulation tools. One of the reasons for this are the very different discretization constrains that govern various effects included in the equations. This paper solves the problem by discretizing the acoustic length scale accurately while treating effects that happen on a smaller length-scale as perturbations.
\end{abstract}

{\it Keywords:} finite-difference time-domain (FDTD), subgrid-scale, quasi-stationary solutions, viscosity effects
\medskip

%
%
%

\section{Introduction}


The application field of Finite-Difference Time-Domain (FDTD)  simulation in acoustics is limited when phenomena on different length scales influence propagation. Grid cells with linear dimensions of about one tenth of the wavelength of interest are required when modelling wave propagation. However, objects and phenomena on a smaller scale like small resonators, sharp edges or viscous boundary layers can influence the total wave propagation significantly. Grid refinements are seldomly an efficient solution to deal with this problem because of the drastic increase of computational costs.

A general approach to deal with different scale phenomena in FDTD simulations is described and tested. The approach is based on a coarse FDTD grid of which the dimensions are only determined by the wavelength considerations. The effects of subgrid phenomena are introduced in the FDTD scheme by adding correction terms to the FDTD equations. These correction terms can result from analytical or numerical models of the subgrid-scale effects. In a recent paper this technique was introduced for subgrid-scale objects. In this manuscript this initial work is generalised and extended to viscosity effects. 


\section{Acoustic equations for simulation models}


Fluid dynamic equations that include all effects that may influence sound
propagation in air are well known. Written in a form that is suitable for acoustic
simulations, they read 
\begin{equation} 
\ddt{\rho_t} + \rho_t \divv\ = -\vb \cdot \nb\rho_t,
\label{mass}
\end{equation}

\begin{equation} 
\rho_t \ddt{\vb} + \nb p_t
=
-\rho_t\left({\vb \cdot \nb}\right)\vb
+\left(\mu_B+{\mu\over{3}}\right) \nb \left(\nb \cdot \vb \right) 
+ \mu \nabla^{2} \vb, 
\label{Euler2}
\end{equation}

and

\begin{equation} 
\begin{array}{c}
{\displaystyle
\ddt{p_t} + \gamma p_t \divv
= - \vb \cdot \nb p_t
+ {(\gamma-1)} \mu_B ({\divv})^2
+ 2{(\gamma-1)}\mu( \phi \cdot \nb) \cdot \vb
}
\\
{\displaystyle
+ {(\gamma-1)}{{\kappa}\over{R}} \nabla^{2}{\left({p_t}\over{\rho_t}\right)}
- {{(\gamma-1)}\over{R}}
\sum_i{{{\rho c_{vi}}\over{R\tau_i}}\left\{{ p_t\over\rho_t}
-\int^{t}{{p_t\over\rho_t}\exp\left({{\tau-t}\over{\tau_i}}\right)
{d\tau\over{\tau_i}}}\right\}}. 
}
\end{array}
\label{pressure2}
\end{equation} 
where
\begin{equation}
\phi_{ij}= 
{1 \over 2}\left({{{\partial v_{ti}}\over{\partial x_j}}+{{\partial
v_{tj}}\over{\partial x_i}}}\right) - {1\over{3}}\left(\divv \right) \delta_{ij},
\label{phidef}
\end{equation}
where ${\mu}$ is the shear viscosity, ${\mu_{B}}$ is the bulk viscosity, ${\kappa}$
is the heat conduction constant and ${c_{vi}}$ and
${\tau_{i}}$ describe molecular vibration. The index $t$ refers to global fields,
including acoustic perturbation. It is not possible however to solve these equations
numerically for practical acoustic situations. One of the difficulties involved
concerns the very different length scales that govern different effects. In fluid
dynamics, approximations for subgrid-scale turbulence are sometimes introduced.
Simular ideas are introduced in this work for other kinds of subgrid-scale phenomena.

Many acoustical phenomena however are accurately described by the linearized 
fluid-dynamic equations, retaining as many terms from Eqs. (\ref{mass}), (\ref{Euler2}),
and (\ref{pressure2}) as required. For situations of interest in this work, a viscous
fluid has to be considered. Both the acoustic velocity
${\bf v}$ and the acoustic pressure $p$ are related by the linearized Euler equation 
 
\begin{equation}
\rho_{o}\frac{\partial {\bf v}}{\partial t} = -\nabla p + \mu [\nabla^{2}{\bf v} + \frac{1}{3} \nabla(\nabla\cdot{\bf v})],
\label{Eqstat1}
\end{equation}
and the linearized mass conservation equation combined with the constitutive relation
$p=\rho c^2$
\begin{equation}
\frac{\partial p}{\partial t} = -\rho_{o}c^2 \nabla \cdot {\bf v},
\label{Eqstat2}
\end{equation}
where $\rho_{o}$ is the equilibrium fluid density %, $\mu$ the fluid viscosity 
and $c$ the speed of sound in the fluidum. The combination of Eq.~(\ref{Eqstat1}) and Eq.~(\ref{Eqstat2}) results in the viscous wave equation. It is well known that for most FDTD applications viscous dissipation is only important near hard boundaries in the so-called viscous boundary layer.

The acoustical field is 'quasi-stationary' in regions that are small compared to the wavelength.  This means that although there is a vibrating motion of the fluid, the field corresponds to a steady flow of the incompressible fluid and the wave behaviour vanishes (i.e. $\nabla \cdot {\bf v}$ goes to zero because $c$ goes to infinity). This assumption leads to the Laplace equation for the acoustic pressure
\begin{equation}
\nabla^2 p \approx 0.
\label{Eqstat4}
\end{equation}
The boundary conditions of this region still depend on the wave equations~(\ref{Eqstat1}) and (\ref{Eqstat2}). 

\section{The FDTD equations for subgrid-scale models}

During FDTD calculations the physical space is discretized in grid cells. For  most applications, the dimensions of the grid cells are determined by the smallest wavelength considered. More explicitly, these dimensions should be less than about one-tenth of this wavelength to keep discretisation error small.  However, when sub-wavelength scale objects are important (for example a small orifice) or when viscous boundary layers are present in the calculation space, the grid cell size must be further reduced, resulting in high computational costs. A
numerically more efficient approach is introduced where the subgrid-scale effects are calculated separately. These effects result in correction terms in the FDTD equations used in the coarse grid.   

The FDTD formalism will be investigated in a staggered Cartesian grid, although there are no theoretical obstacles against applying more general grids. The FDTD equations in a staggered grid will first be derived by integrating the analytical wave equations over the Cartesian grid cells. The integral formulation of the mass equation (Eq.~(\ref{Eqstat2})) over a closed cell volume $V_{i}$ of which the surface $S_{i}$ is divided in a set of sub-surfaces $S_{i,j}$ reads
\begin{equation}
\frac{d\langle  p \rangle_{V_{i}}}{d t} = -\frac{\rho_{o}c^2}{V_{i}} \sum_{j}S_{i,j}\langle v_{n} \rangle_{S_{i,j}},
\label{Eqstat5}
\end{equation}
where $v_{n}$ is the velocity component parallel to the outer normal vector of the surfaces $S_{i,j}$, and $\langle \rangle_{V_{i}}$ and $\langle \rangle_{S_{i,j}}$ refers to the spatial averaging operators over the volume and surfaces of the considered FDTD cell. The sub-surface $S_{i,j}$ is the surface in between volume $V_{i}$ and volume $V_{j}$ and is perpendicular to a coordinate axis. 

A second equation describing the time derivative  $\langle v_{n} \rangle_{S_{i,j}}$ can be derived from the Euler equation (Eq.~(\ref{Eqstat1})),
\begin{equation}
\frac{d \langle v_{n} \rangle_{S_{i,j}}}{d t} =  -\frac{1}{\rho_{o}}\langle \frac{\partial p }{\partial n}\rangle_{S_{i,j}} + \frac{\mu}{\rho_{o}} \langle (\nabla^{2}v_{n} + \frac{1}{3} \frac{\partial  }{\partial n} (\nabla\cdot{\bf v})\rangle_{S_{i,j}}.
\label{Eqstat6}
\end{equation}
To obtain a closed system of equations, both the spatially averaged normal derivative of the pressure and the viscosity terms has to be estimated from the mean pressure values over the volumes $V_{i}$. 

In the absence of short scale effects, the spatial variations of the pressure is quasi linear over the grid cells, resulting in a straightforward approximation for $\langle \frac{\partial p }{\partial n}\rangle_{S_{i,j}}$ as a linear function of the mean pressure values over the volumes $V_{i}$ and $V_{j}$.  This way the standard FDTD equations can be obtained. However when sub-wavelength geometries and/or viscous boundary layers are present the linear relation with the volume pressure values is not sufficient and correction terms are necessary.

\section{Conclusions}

The extension of the FDTD equations with correction terms to deal with subgrid-scale phenomena allows the efficient calculation of wave propagation affected by these phenomena. The methodology to determine correction coefficients for viscous boundary layers and subgrid-scale hard objects is derived and tested. Both techniques seem sufficiently accurate when they are compared to theoretical and experimental results. However, when the two techniques are combined to model viscous boundary layers on subgrid-scale hard objects some further approximations have to be made, leading to an reduced accuracy of the viscous dissipation effects. The technique is effectively used in some practical situations.  




\end{document}





