Physical domain
We consider a one-dimensional system where masses are connected with springs and dampers.
For simplicity, we assume that the mass \(\dimvar{m}\), the spring constant \(\dimvar{k}\), and the damping coefficient \(\dimvar{\mu}\) take identical values for now.
Due to the force balance for each mass, we obtain a set of equations which govern the dynamics as
\[ \begin{align}\begin{aligned}\tder{
\dimvar{\pos{}{}}
}{
\dimvar{t}
}
=
&
\dimvar{\vel{}{}},\\\dimvar{m}
\tder{
\dimvar{\vel{}{}}
}{
\dimvar{t}
}
=
&
\dimvar{k}
\left\{
\dimvar{\pos{}{}} \left( \dimvar{x} - \Delta \dimvar{x} \right)
-
\dimvar{\pos{}{}} \left( \dimvar{x} \right)
\right\}
+
\dimvar{k}
\left\{
\dimvar{\pos{}{}} \left( \dimvar{x} + \Delta \dimvar{x} \right)
-
\dimvar{\pos{}{}} \left( \dimvar{x} \right)
\right\}\\+
&
\dimvar{\mu}
\left\{
\dimvar{\vel{}{}} \left( \dimvar{x} - \Delta \dimvar{x} \right)
-
\dimvar{\vel{}{}} \left( \dimvar{x} \right)
\right\}
+
\dimvar{\mu}
\left\{
\dimvar{\vel{}{}} \left( \dimvar{x} + \Delta \dimvar{x} \right)
-
\dimvar{\vel{}{}} \left( \dimvar{x} \right)
\right\}\\+
&
\dimvar{m} \dimvar{a},\end{aligned}\end{align} \]
where \(\dimvar{\pos{}{}}\) and \(\dimvar{\vel{}{}}\) are the displacements and velocities, respectively, and \(\dimvar{m} \dimvar{a}\) is the external forcing.
In the limit of \(\Delta \dimvar{x} \rightarrow 0\), by utilizing the Taylor series expansion, the second relation converges to
\[\pder{
\dimvar{\vel{}{}}
}{
\dimvar{t}
}
=
\dimvar{c}^2
\pder{}{\dimvar{x}}
\pder{
\dimvar{\pos{}{}}
}{
\dimvar{x}
}
+
\dimvar{\nu}
\pder{}{\dimvar{x}}
\pder{
\dimvar{\vel{}{}}
}{
\dimvar{x}
}
+
\dimvar{a},\]
where
\[ \begin{align}\begin{aligned}\dimvar{c}
&
\equiv
\sqrt{
\frac{
\dimvar{k} \Delta \dimvar{x}^2
}{
\dimvar{m}
}
},\\\dimvar{\nu}
&
\equiv
\frac{
\dimvar{\mu} \Delta \dimvar{x}^2
}{
\dimvar{m}
}\end{aligned}\end{align} \]
are the propagation speed and the kinematic viscosity having the dimensions of \(\left[ L T^{-1} \right]\) and \(\left[ L^2 T^{-1} \right]\), respectively, which are to be given as material properties.
By using a reference length scale \(\dimvar{L}\) and a time scale \(\dimvar{T}\), we obtain a set of non-dimensional governing equations:
\[ \begin{align}\begin{aligned}\pder{\pos{}{}}{t}
&
=
\vel{}{},\\\pder{\vel{}{}}{t}
&
=
c^2
\pder{}{x}
\pder{\pos{}{}}{x}
+
\nu
\pder{}{x}
\pder{\vel{}{}}{x}
+
a,\end{aligned}\end{align} \]
where the following non-dimensional numbers are involved:
\[ \begin{align}\begin{aligned}t
&
\equiv
\frac{
1
}{
\dimvar{T}
}
\dimvar{t},\\\pos{}{}
&
\equiv
\frac{
1
}{
\dimvar{L}
}
\dimvar{\pos{}{}},\\\vel{}{}
&
\equiv
\frac{
\dimvar{T}
}{
\dimvar{L}
}
\dimvar{\vel{}{}},\\c
&
\equiv
\frac{
\dimvar{T}
}{
\dimvar{L}
}
\dimvar{c},\\\nu
&
\equiv
\frac{
\dimvar{T}
}{
\dimvar{L}^2
}
\dimvar{\nu},\\a
&
\equiv
\frac{
\dimvar{T}^2
}{
\dimvar{L}
}
\dimvar{a}.\end{aligned}\end{align} \]
For multi-dimensional domains, we have a straightforward extension:
\[ \begin{align}\begin{aligned}\pder{\pos{}{}}{t}
&
=
\vel{}{},\\\pder{\vel{}{}}{t}
&
=
c^2
\pder{}{x_i}
\pder{\pos{}{}}{x_i}
+
\nu
\pder{}{x_i}
\pder{\vel{}{}}{x_i}
+
a,\end{aligned}\end{align} \]
where the summation rule is adopted for brevity.
To see the energy balance, we multiply the equation of motion by \(\rho \vel{}{}\) and volume-integrate it:
\[\myint{0}{L}{
\rho
\vel{}{}
\pder{\vel{}{}}{t}
}{x}
-
\myint{0}{L}{
\rho
\vel{}{}
c^2
\pder{}{x}
\pder{\pos{}{}}{x}
}{x}
-
\myint{0}{L}{
\rho
\vel{}{}
\nu
\pder{}{x}
\pder{\vel{}{}}{x}
}{x}
-
\myint{0}{L}{
\rho
\vel{}{}
a
}{x}
=
0,\]
where \(L\) is the length of the domain.
With some algebra, we see that the first term is equal to
\[\pder{}{t}
\left(
\myint{0}{L}{
\frac{1}{2}
\rho
\vel{}{}^2
}{x}
\right)
=
\tder{K}{t},\]
which describes the change of the total kinetic energy \(K\) in time.
To reformulate the second term, we adopt the integration by parts:
\[-
\myint{0}{L}{
\rho
c^2
\vel{}{}
\pder{}{x}
\pder{\pos{}{}}{x}
}{x}
=
-
\left[
\rho
c^2
\vel{}{}
\pder{\pos{}{}}{x}
\right]_{0}^{L}
+
\myint{0}{L}{
\rho
c^2
\pder{\vel{}{}}{x}
\pder{\pos{}{}}{x}
}{x},\]
where the former term denotes the boundary contribution, while the latter term is attributed to the change of the total potential energy \(U\) in time as it can be written as
\[\myint{0}{L}{
\rho
c^2
\pder{}{t}
\left(
\pder{\pos{}{}}{x}
\right)
\pder{\pos{}{}}{x}
}{x}
=
\pder{}{t}
\left(
\myint{0}{L}{
\frac{1}{2}
\rho
c^2
\pder{\pos{}{}}{x}
\pder{\pos{}{}}{x}
}{x}
\right)
=
\tder{U}{t}.\]
Similarly, the third term can be reformulated by adopting the integration by parts:
\[-
\myint{0}{L}{
\rho
\nu
\vel{}{}
\pder{}{x}
\pder{\vel{}{}}{x}
}{x}
=
-
\left[
\rho
\nu
\vel{}{}
\pder{\vel{}{}}{x}
\right]_{0}^{L}
+
\myint{0}{L}{
\rho
\nu
\pder{\vel{}{}}{x}
\pder{\vel{}{}}{x}
}{x},\]
where the former term is again the boundary contribution whereas the latter denotes the energy dissipation.
The fourth term simply denotes the energy input due to the external factor, and we remain it as it is.
In summary, we obtain
\[\tder{E}{t}
=
\left[
\rho
c^2
\vel{}{}
\pder{\pos{}{}}{x}
\right]_{0}^{L}
+
\left[
\rho
\nu
\vel{}{}
\pder{\vel{}{}}{x}
\right]_{0}^{L}
-
\myint{0}{L}{
\rho
\nu
\pder{\vel{}{}}{x}
\pder{\vel{}{}}{x}
}{x}
+
\myint{0}{L}{
\rho
\vel{}{}
a
}{x},\]
where \(E \equiv K + U\) is the total energy of the system.
Frequency domain
In this project, we assume the domain is periodic in all directions for simplicity.
This choice allows us to expand the fields using trigonometric functions:
\[ \begin{align}\begin{aligned}\pos{}{} \left( x, t \right)
&
=
\mysum{l}{-\infty}{\infty}
\kpos{l}{} \left( t \right)
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right),\\\vel{}{} \left( x, t \right)
&
=
\mysum{l}{-\infty}{\infty}
\kvel{l}{} \left( t \right)
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right).\end{aligned}\end{align} \]
Although \(\kpos{l}{}\) and \(\kvel{k}{}\) are both function in time, they are not explicitly indicated later for brevity.
Assigning them to the governing equations yields
\[ \begin{align}\begin{aligned}\mysum{l}{-\infty}{\infty}
\tder{\kpos{l}{}}{t}
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right)
=
&
\mysum{l}{-\infty}{\infty}
\kvel{l}{}
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right),\\\mysum{l}{-\infty}{\infty}
\tder{\kvel{l}{}}{t}
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right)
=
&
-
\mysum{l}{-\infty}{\infty}
c^2
\left( \wavenumber{2 \pi}{l}{L} \right)^2
\kpos{l}{}
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right)\\&
-
\mysum{l}{-\infty}{\infty}
\nu
\left( \wavenumber{2 \pi}{l}{L} \right)^2
\kvel{l}{}
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right)\\&
+
\mysum{l}{-\infty}{\infty}
A_l
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right).\end{aligned}\end{align} \]
Due to the orthogonality:
\[ \begin{align}\begin{aligned}&
\myint{0}{L}{
q \left( x, t \right)
\exp \left( \wavenumber{- 2 \pi}{k x}{L} I \right)
}{x}\\=
&
\myint{0}{L}{
\mysum{l}{-\infty}{\infty}
Q_l \left( t \right)
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right)
\exp \left( \wavenumber{- 2 \pi}{k x}{L} I \right)
}{x}\\=
&
\myint{0}{L}{
\mysum{l}{-\infty}{\infty}
Q_l \left( t \right)
\exp \left( \wavenumber{2 \pi}{\left( l - k \right) x}{L} I \right)
}{x}\\=
&
L
Q_{k} \left( t \right),\end{aligned}\end{align} \]
we obtain
\[ \begin{align}\begin{aligned}\tder{\kpos{k}{}}{t}
=
&
\kvel{k}{},\\\tder{\kvel{k}{}}{t}
=
&
-
C_k
\kpos{k}{}
-
N_k
\kvel{k}{}
+
A_k,\end{aligned}\end{align} \]
where
\[ \begin{align}\begin{aligned}C_k
&
\equiv
c^2
\left( \wavenumber{2 \pi}{k}{L} \right)^2,\\N_k
&
\equiv
\nu
\left( \wavenumber{2 \pi}{k}{L} \right)^2.\end{aligned}\end{align} \]
For two-dimensional cases, we have
\[ \begin{align}\begin{aligned}\tder{\kpos{k_x,k_y}{}}{t}
=
&
\kvel{k_x,k_y}{},\\\tder{\kvel{k_x,k_y}{}}{t}
=
&
-
C_{k_x,k_y}
\kpos{k_x,k_y}{}
-
N_{k_x,k_y}
\kvel{k_x,k_y}{}
+
A_{k_x,k_y},\end{aligned}\end{align} \]
where the following symbols are defined for notational simplicity:
\[ \begin{align}\begin{aligned}C_{k_x,k_y}
&
\equiv
c^2
\left( \wavenumber{2 \pi}{k_x}{L_x} \right)^2
+
c^2
\left( \wavenumber{2 \pi}{k_y}{L_y} \right)^2,\\N_{k_x,k_y}
&
\equiv
\nu
\left( \wavenumber{2 \pi}{k_x}{L_x} \right)^2
+
\nu
\left( \wavenumber{2 \pi}{k_y}{L_y} \right)^2.\end{aligned}\end{align} \]
Higher dimensional relations can be obtained in the same manner.
The kinetic and potential energies using the variables in frequency domain lead to
\[ \begin{align}\begin{aligned}K
&
\equiv
\myint{0}{L}{
\frac{1}{2}
\rho
\vel{}{}^2
}{x}\\&
=
\myint{0}{L}{
\frac{1}{2}
\rho
\left(
\mysum{k}{-\infty}{\infty}
\kvel{k}{}
\exp \left( \wavenumber{2 \pi}{k x}{L} I \right)
\right)
\left(
\mysum{l}{-\infty}{\infty}
\kvel{l}{}
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right)
\right)
}{x}\\&
=
\mysum{k}{-\infty}{\infty}
\frac{1}{2}
\rho
L
\left( \kvel{k}{} \right)^2,\\U
&
\equiv
\myint{0}{L}{
\frac{1}{2}
\rho
c^2
\pder{\pos{}{}}{x}
\pder{\pos{}{}}{x}
}{x}\\&
=
\myint{0}{L}{
\frac{1}{2}
\rho
c^2
\left(
\mysum{k}{-\infty}{\infty}
\kpos{k}{}
\wavenumber{2 \pi}{k}{L} I
\exp \left( \wavenumber{2 \pi}{k x}{L} I \right)
\right)
\left(
\mysum{l}{-\infty}{\infty}
\kpos{l}{}
\wavenumber{2 \pi}{l}{L} I
\exp \left( \wavenumber{2 \pi}{l x}{L} I \right)
\right)
}{x}\\&
=
\mysum{k}{-\infty}{\infty}
\frac{1}{2}
\rho
L
C_k
\left( \kpos{k}{} \right)^2,\end{aligned}\end{align} \]
respectively.
Higher dimensional cases can be derived in the same manner:
\[ \begin{align}\begin{aligned}K
&
=
\mysum{k_y}{-\infty}{\infty}
\mysum{k_x}{-\infty}{\infty}
L_x
L_y
\frac{1}{2}
\rho
\left( \kvel{k_x, k_y}{} \right)^2,\\U
&
=
\mysum{k_y}{-\infty}{\infty}
\mysum{k_x}{-\infty}{\infty}
L_x
L_y
\frac{1}{2}
\rho
C_{k_x, k_y}
\left( \kpos{k_x, k_y}{} \right)^2.\end{aligned}\end{align} \]