This program was contributed by Andreas Hegendörfer <[email protected]>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
Pictures from this code gallery program
Annotated version of Readme.md
Overview
The Sound Transmission Loss (STL) is a measure of the sound insulation of components and is therefore an important acoustic quantity. Here, the STL of a concrete wall is calculated by means of vibroacoustic FEM simulations. In the frequency domain, the elastodynamic equations are solved in the structural domain, while the Helmholtz equation is solved in the acoustic domain, with full coupling at the interface. The deal.II implementation is parallelized using MPI and MUMPS, employs the hp finite element method, and is mainly based on step-62 step-46, and step-27.
Elastodynamic equations
The behavior of the structure \( \Omega_S \) is modelled by means of the linear elastodynamic equations, which are given in strong form as
\[
\rho_S \ddot{u}_i-\left( c_{ijkl} \epsilon_{kl} \right)_{,j} = f_i \qquad \text{in } \Omega_S \\
\qquad \qquad \qquad u_i = \bar{u}_i \qquad \text{on } \partial \Omega_{SD} \\
\qquad \qquad \left( c_{ijkl} \epsilon_{kl} \right) n_j= \bar{t}_i \qquad \text{on } \partial \Omega_{SN} \\
\]
Here, \(u_i \) denotes the displacement, \( \rho \) is the mass density, \( f_i \) is the body force density, \( c_{ijkl} \) is the stiffness tensor, \( n_i \) is the outwards pointing normal unit vector, \( \bar{t}_i \) is the prescribed traction and \( \bar{u}_i \) is the prescribed displacement and the double dot represents the second derivative with respect to time. The linear strain tensor \( \epsilon_{ij} \) is defined as
\[ \epsilon_{ij} = \frac{1}{2}\left( u_{i,j} + u_{j,i} \right) \]
Multiplying by the test function \( \phi_i \), integration over the domain, applying integration by parts and the divergence theorem and introducing the boundary conditions results in the weak form
\begin{align*}
\int_{\Omega_S} \phi_j \rho_S \ddot{u}_j dV + \int_{\Omega_S} \phi_{j,i} c_{ijkl} \epsilon_{kl} dV - \int_{\partial \Omega_{SN} } \phi_j \bar{t}_j = 0
\end{align*}
The time domain equation above is transferred into the frequency domain by performing a Fourier transform with regard to the time variable. Then, the weak form becomes
\begin{align*}
-\int_{\Omega_S} \phi_j \omega^2 \rho_S {u}_j dV + \left[ 1 + i\mathcal{I} \right] \int_{\Omega_S} \phi_{j,i} c_{ijkl} \epsilon_{kl} dV - \int_{\partial \Omega_{SN} } \phi_j \bar{t}_j = 0
\end{align*}
whereby \( \omega \) is the angular frequency. Here, the isotropic loss factor \( \mathcal{I} \) is introduced to account for structural damping and \( i \) is the imaginary unit defined as \( i^2=-1 \).
Acoustic equations
The acoustic equation considered is the wave equation for pressure, given in strong form as
\[
\frac{1}{c^2} \ddot{p} - p_{,ii} = 0 \qquad \text{in } \Omega_A \\
\qquad \qquad p = \bar{p} \qquad \text{on } \partial \Omega_{AD} \\
\qquad \qquad p_{,i} n_i = -\rho_A \dot{\bar{v}}_n \qquad \text{on } \partial \Omega_{AN} \\
\]
Here, \( c \) is the sound speed and \( p \) is the acoustic pressure. As boundary conditions the prescribed normal velocity \( \bar{v}_n \) and the prescribed acoustic pressure \(\bar p \) are considered. Similar as for the elastodynamic equations, the weak form results as
\begin{align*}
\int_{\Omega_A} \eta \frac{1}{c^2} \ddot{p} dV + \int_{\Omega_A} \eta_{,i} p_{,i} dV - \int_{\partial \Omega_{AN}} \eta \rho_A \dot{\bar{v}}_n dA = 0
\end{align*}
whereby \( \eta \) denotes the test function. Performing a Fourier transform with regard to the time variable results in the frequency domain weak form, given with the wave number \( k = \omega / c \) as
\begin{align*}
-\int_{\Omega_A} \eta k^2 {p} dV + \int_{\Omega_A} \eta_{,i} p_{,i} dV - \int_{\partial \Omega_{AN}} \eta \rho_A i \omega {\bar{v}}_n dA = 0
\end{align*}
A perfectly matched layer (PML) is used to simulate an open boundary. A PML is a complex coordinate stretch in the form of
\begin{align*}
x \to x + \frac{i}{\omega} \int^{x} \sigma(x') dx'
\end{align*}
where \( \sigma \) is the PML absorption function. It vanishes outside of the PML, where no absorption is desired but is turned on inside the PML resulting in an absorption of outgoing waves while minimizing reflections back into the domain. Considering the PML, the weak form becomes
\[
-\int_{\Omega_A} \eta k^2 {p} J dV + \int_{\Omega_A} \eta_{,i} p_{,i} \Lambda_i^2 J dV - \int_{\partial \Omega_{AN}} \eta \rho_A i \omega {\bar{v}}_n \frac{J}{\Lambda_n} dA = 0 \\
J = \prod_{i=1}^{\mathrm{dim}} s_i \\
\Lambda_i = \frac{1}{s_i} \\
s_i = 1+\frac{i}{\omega} \sigma_i
\]
In the equations above the transformation of both volume and surface elements is taken into account.
Coupled Mechanic-Acoustic equations
For the structure-acoustic coupling, the normal vector \( n_{i} \) is chosen to conform with the normal vector of the mechanical domain given as
\[
n_{i} = n_{i_S} = -n_{i_A}
\]
Here, the normal vector \( n_{i_S} \) referes to the mechanical domain, while \( n_{i_A} \) denotes the normal vector of the acoustic domain.
The coupling between the mechanic and the acoustic domain takes place at the interface \( \partial \Omega_{SA} \) and is modelled via boundary conditions, specified as
\[
-p n_{i} = \bar{t}_i \qquad \text{on } \partial \Omega_{SA} \\
p_{,i} n_i = \omega^2 \rho_A u_i n_i \qquad \text{on } \partial \Omega_{SA}
\]
Then, the coupled system becomes
\[
-\int_{\Omega_S} \phi_j \omega^2 \rho_S {u}_j dV + \left[ 1 + i\mathcal{I} \right] \int_{\Omega_S} \phi_{j,i} c_{ijkl} \epsilon_{kl} dV - \int_{\partial \Omega_{SN} } \phi_j \bar{t}_j \\
- \int_{\Omega_A} \eta k^2 {p} J dV + \int_{\Omega_A} \eta_{,i} p_{,i} \Lambda_i^2 J dV - \int_{\partial \Omega_{AN}} \eta \rho_A i \omega {\bar{v}}_n \frac{J}{\Lambda_n} dA \\
+ \int_{\partial \Omega_{SA}} \phi_i p n_{i} dA + \int_{\partial \Omega_{SA}} \eta \rho_A \omega^2 u_i n_{i} dA = 0
\]
The coupled system above is implemented in the code.
Problem description and modeling
The STL of a concrete wall is calculated. Conceptually, several methods exist for calculating the STL, including modeling the source and receiver rooms. Here, however, in the interest of computational efficiency the source room is not modeled and the wall is directly excited by a prescribed diffuse sound field with sound power \( P_s\) introducing structural vibrations. These vibrations radiate airborne sound with sound power \( P_r\) into the receiver side, which is considered anechoic and modeled using PMLs, as illustrated in the figure below.
The same problem and modeling setup as in [1] and [2] are considered here.
The geometric dimension of the wall are: height \(H = 4.37\,{m}\), width \(W = 2.84\,\text{m}\) and thickness \(T = 0.203\,\text{m}\). The density of the concrete is \( \rho_S = 2275\,\frac{kg}{m^3} \), its Young’s modulus is \( E = 31.6\,GPa \) and the Poisson’s ratio is \(\nu = 0.2\). As damping, an isotropic loss factor of \(\mathcal{I} = 0.01\) is assumed. The outer boundary of the wall is fixed, and the surrounding wall is considered sound-hard and therefore does not influence the STL. The density of the air is \( \rho_A = 1.204\,\frac{kg}{m^3} \) and the speed of sound in air is \( c = 343\,\frac{m}{s} \).
The diffuse sound pressure field \( p_{diffuse,wall} \) acting on the source side of the wall surface must be admissible, that is, it must satisfy the wave equation. Assuming \( N \) uncorrelated plane waves traveling in \( x\)-direction, an admissible diffuse sound pressure field \( p_{diffuse,wall} \) on the wall surface can be obtained as
\[
p_{diffuse,wall} = \underbrace{ \frac{A}{\sqrt{2N}} \sum^{N}_{n=1} e^{-i\left( k_{n,x}x + k_{n,y}y + k_{n,z}z \right)} e^{i \Phi_n}}_{p_{incident}} + \underbrace{\frac{A}{\sqrt{2N}} \sum^{N}_{n=1} e^{-i\left( -k_{n,x}x + k_{n,y}y + k_{n,z}z \right)} e^{i \Phi_n}}_{p_{reflection}} \\
k_{n,x} = \text{cos}(\Theta_n) \\
k_{n,y} = \text{sin}(\Theta_n) \text{cos}(\phi_n) \\
k_{n,z} = \text{sin}(\Theta_n) \text{sin}(\phi_n)
\]
Here, \( A \) is the amplitude of the plane waves, the polar angles \(\Theta_n\), \(\phi_n\) and the phase angle \(\Phi_n\) are uniformly distributed random variables in the interval \( \left[0,2\pi \right] \), the polar angle \( \Theta_n \) follows from \( \Theta_n = \text{acos}(q_n)\), whereby \( q_n\) is a uniformly distributed variable in the interval \( \left[0,1 \right] \). On the wall surface the diffuse sound pressure field \( p_{diffuse,wall} \) is a superposition of the incoming waves \( p_{incident}\) and their reflections \( p_{reflection}\). In the equations above, a perfectly reflecting wall is assumed, so no reduction of amplitude or phase shift occurs. To overcome this simplification, the source room would need to be modeled, which is not done here.
The STL is defined as
\[
STL = 10 \text{log}_{10} \left(\frac{P_s}{P_r} \right)
\]
and therefore requires the calculation of the sound powers on the source \( P_s\) and receiver \( P_r\) sides.
The sound power on the source side \( P_s\) can be calculated as
\[
P_s = \frac{1}{2}\int_{\partial A_{wall,source} } \text{Re}\left[ p_{incident} v_{i}^{*} n_i \right] dA
\]
Here, the asterix ( \( {}^*\)) indicates complex conjugation, \( A_{wall,source} \) is the surface of the wall on the source side and \( v_i \) is the sound particle velocity, which can be calculated as
\[
v_{i} = \frac{-1}{i\omega \rho_a}p_{incident,i}
\]
The sound power on the receiver side \( P_r\) can be calculated as
\[
P_r = \frac{1}{2}\int_{\partial A_{wall,receiver} } \text{Re}\left[ p \left(i\omega u_i\right)^{*} n_i \right] dA
\]
Here, \( A_{wall,receiver} \) is the surface of the wall on the receiver side.
Results
In the figure below, the STL obtained in the present work is compared with measurement data from [3], with results from a commercial FEM software reported in [1] for a single random seed in the diffuse sound field, and with the corresponding average STL over multiple random seeds as presented in [2]. Overall, the STL calculations are sensitive to the choice of random seeds for the diffuse sound field. The STL ranges between 30 and 60 dB, which corresponds to the incident sound power being between \(10^3\) and \(10^6\) times greater than the transmitted sound power, depending on frequency.
In addition, the figure below shows the magnitude of the solution for the first resonance frequency at around 112 Hz. On the left, the magnitude of the displacement \(u_i\) is presented, while the magnitude of the sound pressure \(p\) directly behind the concrete wall is illustrated on the right side. As enforced by the PML, the sound pressure \(p\) vanishes at the boundaries of the domain.
Ideas for extensions
Implementation of a QuadratureCache class to store calculated data for single cells to be reused for different frequencies, similar to that in step-62. Application of an iterative solver. Consideration of non-matching grids between structural and acoustic domains.
References
[1] COMSOL Multiphysics 6.4, "Sound Transmission Loss Through a Concrete Wall", available at https://www.comsol.com/model/sound-transmission-loss-through-a-concrete-wall-73371.
[2] M.H. Jensen, "Modeling Sound Transmission Loss Through a Concrete Wall", COMSOL Blog, 2020, available at https://www.comsol.com/blogs/modeling-sound-transmission-loss-through-a-concrete-wall.
[3] A. Litvin and H.W. Belliston, “Sound Transmission Loss Through Concrete and Concrete Masonry Walls,” American Concrete Institute, Journal Proceedings, vol. 45, pp. 641–646, 1978.
Compilation and running the code
To run the code, deal.II must be installed together with PETSc (configured with complex number support) and the p4est library.
Similar to the deal.ii examples, run
cmake -DDEAL_II_DIR=/path/to/deal.II .
in this directory to configure the project.
You can switch between debug and release mode by calling either
or
To execute the program in serial run
and for parallel execution (in this case, on j processors) run
mpirun -np j ./parallel_vibroacoustic_solver
Annotated version of parallel_vibroacoustic_solver.cc
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
namespace VibroAcousticProblem
* * * struct InterferenceTaperTransform *
geometric tolerance
constexpr double geom_tol = 1e-10;
By default, boundary indicators are 0
return the isotropic loss factor
return std::complex<double>{1., 0.01};
calculation of the linear strain tensor
const unsigned int shape_func,
const unsigned int q_point)
{
for (
unsigned int i = 0; i < dim; ++i)
for (
unsigned int i = 0; i < dim; ++i)
for (
unsigned int j = i + 1; j < dim; ++j)
Tensor< 1, spacedim > shape_grad_component(const unsigned int i, const unsigned int q_point, const unsigned int component) const
Returning the stiffness tensor of the wall
const double E = 31600. * 1e6;
double lambda = v / (1 - 2 * v) * 1 / (1 + v) * E;
double mu = 0.5 * 1 / (1 + v) * E;
for (
unsigned int i = 0; i < dim; ++i)
for (
unsigned int j = 0; j < dim; ++j)
for (
unsigned int k = 0; k < dim; ++k)
for (
unsigned int l = 0; l < dim; ++l)
stiffness_tensor[i][j][k][l] =
(((i == k) && (j == l) ? mu : 0.0) +
((i == l) && (j == k) ? mu : 0.0) +
((i == j) && (k == l) ? lambda : 0.0));
Returning the density of the wall
Returning the density of air
Returning the speed of sound
Implementation of a perfectly matched layer
class PML :
public Function<dim, std::complex<double>>
explicit PML(
double omega)
: omega(omega)
, pml_coeff(1.e4 / omega)
Assert(dim == 3, ExcNotImplemented());
#define Assert(cond, exc)
x-direction
b_neg[0] = std::numeric_limits<double>::lowest();
y-direction
z-direction
calculates the value of the complex coordinate stretch
Vector<std::complex<double>> &value)
const override;
std::array<double, dim> b_pos, b_neg;
std::array<double, dim> t_pos, t_neg;
Vector<std::complex<double>> &value)
const
for (
unsigned int d = 0; d < dim; ++d)
const double x_prime = p[d] - b_pos[d];
pml_coeff /
std::pow(t_pos[d], pml_coeff_degree);
coeff = a_coeff *
std::pow(x_prime, pml_coeff_degree);
else if (p[d] < b_neg[d])
const double x_prime = b_neg[
d] - p[
d];
pml_coeff /
std::pow(t_neg[d], pml_coeff_degree);
coeff = a_coeff *
std::pow(x_prime, pml_coeff_degree);
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
complex coordinate stretching: s = 1 + i * sigma(x)
value[d] = std::complex<double>(1.0, coeff);
Creation of a diffuse sound field for excitation of the wall on the source side
class DiffuseSoundField :
public Function<dim, std::complex<double>>
DiffuseSoundField(
unsigned int N,
double omega,
MPI_Comm mpi_communicator)
, dist_0_2PI(0.0, 2. *
numbers::PI)
, generator(static_cast<
int>(omega))
Only let one rank 0 create the random variables of the diffuse sound field...
const unsigned int rank =
for (
unsigned int n = 0; n < N; ++n)
phi[n] = dist_0_2PI(generator);
phase[n] = dist_0_2PI(generator);
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
inline ::VectorizedArray< Number, width > acos(const ::VectorizedArray< Number, width > &x)
... and broadcast from rank 0 to all ranks.
T broadcast(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
Generation of kn on all ranks
for (
unsigned int n = 0; n < N; n++)
double scale = (omega / get_sound_speed());
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
Calculation of the total pressure on the source side
std::vector<std::complex<double>> &values,
const unsigned int component = 0)
const override;
Calculation of the gradient of the total pressure on the source side
gradient_list(
const std::vector<
Point<dim>> &points,
std::vector<
Tensor<1, dim, std::complex<double>>> &gradients,
const unsigned int component = 0)
const override;
Calculation of the incident sound pressure on the source side
value_list_incidence(
const std::vector<
Point<dim>> &points,
std::vector<std::complex<double>> &values,
const unsigned int component = 0)
const;
Number of plane waves
For randmoness
std::uniform_real_distribution<double> dist_0_2PI;
std::uniform_real_distribution<double> dist_0_1;
variables of the diffuse field
std::vector<double> phi, Phi, phase, kn_x, kn_y, kn_z;
std::vector<Tensor<1, dim, double>> k_vector;
angular velocity
DiffuseSoundField<dim>::value_list(
const std::vector<
Point<dim>> &points,
std::vector<std::complex<double>> &values,
const unsigned int component)
const
AssertThrow(component == 0, ExcMessage(
"only component 0 is implemented"));
values.resize(points.size());
std::complex<double> j{0., 1.};
for (
unsigned int q = 0; q < points.size(); ++q)
const auto &q_p = points[q];
for (
unsigned int n = 0; n <
N; n++)
#define AssertThrow(cond, exc)
sound waves traveling towards the wall...
double dot_reflection = 0.;
for (
unsigned int d = 0; d < dim; d++)
sound waves traveling towards the wall...
dot_towards += k_vector[n][d] * q_p[d];
... and reflections.
dot_reflection -= k_vector[n][d] * q_p[d];
dot_reflection += k_vector[n][d] * q_p[d];
sound waves traveling towards the wall...
values[q] +=
std::exp(-j * (dot_towards) + j * phase[n]);
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
... and reflections.
values[q] +=
std::exp(-j * (dot_reflection) + j * phase[n]);
values[q] *= 1. / (
std::sqrt(2. *
static_cast<double>(N))) * 1.e6;
DiffuseSoundField<dim>::value_list_incidence(
std::vector<std::complex<double>> &values,
const unsigned int component)
const
AssertThrow(component == 0, ExcMessage(
"only component 0 is implemented"));
std::complex<double> j{0., 1.};
for (
unsigned int q = 0; q < points.size(); ++q)
const auto &q_p = points[q];
for (
unsigned int n = 0; n <
N; n++)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Only sound waves traveling towards the wall.
for (
unsigned int d = 0; d < dim; d++)
Sound waves traveling towards the wall.
dot_towards += k_vector[n][d] * q_p[d];
Only sound waves traveling towards the wall.
values[q] +=
std::exp(-j * (dot_towards) + j * phase[n]);
values[q] *= 1. / (
std::sqrt(2. *
static_cast<double>(N))) * 1.e6;
DiffuseSoundField<dim>::gradient_list(
std::vector<
Tensor<1, dim, std::complex<double>>> &gradients,
const unsigned int component)
const
AssertThrow(component == 0, ExcMessage(
"only component 0 is implemented"));
const std::complex<double> j{0.0, 1.0};
for (
unsigned int q = 0; q < points.size(); ++q)
const auto &q_p = points[q];
for (
unsigned int d = 0;
d < dim; ++
d)
for (
unsigned int n = 0; n <
N; ++n)
for (
unsigned int d = 0;
d < dim; ++
d)
dot += k_vector[n][d] * q_p[d];
const std::complex<double> exp_term =
for (
unsigned int d = 0;
d < dim; ++
d)
gradients[q][d] += (-j * k_vector[n][d]) * exp_term *
(1. /
std::sqrt(2. *
static_cast<double>(N))) *
HarmonicResponse(
double omega);
run(
bool write_output =
false);
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
calculate the magnitude of u and p
calculation of sound power at receiver side for one cell
cell_receiver_sound_power(
calculation of sound power at sender side
calculation of sound power at receiver side
Assemble the air-structure coupling terms This function implements the assembly of the air-structure interface.
assemble_air_structure_interface_term(
FullMatrix<std::complex<double>> &local_interface_matrix,
const QGauss<dim> quadrature_formula_structure, quadrature_formula_air;
const QGauss<dim - 1> face_quadrature_formula_structure,
face_quadrature_formula_air;
DiffuseSoundField<dim> field;
HarmonicResponse<dim>::HarmonicResponse(
double omega)
: mpi_communicator(MPI_COMM_WORLD)
, triangulation(mpi_communicator,
, quadrature_formula_structure(fe_structure.degree + 1)
, quadrature_formula_air(fe_air.degree + 1)
, face_quadrature_formula_structure(fe_structure.degree + 1)
, face_quadrature_formula_air(fe_air.degree + 1)
, dof_handler(triangulation)
, field(1.e3, omega, mpi_communicator)
"HarmonicResponse is only implemented for dim == 3");
fe_collection.push_back(fe_structure);
fe_collection.push_back(fe_air);
q_collection.push_back(quadrature_formula_structure);
q_collection.push_back(quadrature_formula_air);
q_face_collection.push_back(face_quadrature_formula_structure);
q_face_collection.push_back(face_quadrature_formula_air);
HarmonicResponse<dim>::setup_system()
Set material id and active FE indices.
for (
const auto &cell : dof_handler.cell_iterators())
cell->set_material_id(static_cast<unsigned
int>(MaterialID::Concrete));
if ((cell->center()[0] - 203.) > 0.)
cell->set_material_id(
static_cast<unsigned int>(MaterialID::Air));
for (
const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
cell->set_active_fe_index(
static_cast<unsigned
int>(MaterialID::Concrete));
if ((cell->center()[0] - 203.) > 0.)
if (cell->is_locally_owned())
cell->set_active_fe_index(
static_cast<unsigned int>(MaterialID::Air));
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
Definition of FE space
dof_handler.distribute_dofs(fe_collection);
pcout <<
" Number of degrees of freedom = " << dof_handler.n_dofs()
locally_owned_dofs = dof_handler.locally_owned_dofs();
locally_relevant_dofs.clear();
locally_relevant_solution.reinit(locally_owned_dofs,
locally_relevant_magnitude.reinit(locally_owned_dofs, mpi_communicator);
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
set up contraints
constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Set boundary ids
for (
const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (unsigned
int f = 0; f < cell->n_faces(); ++f)
const auto face = cell->face(f);
const auto p = face->center();
static_cast<unsigned int>(MaterialID::Concrete)))
cell->face(f)->set_user_index(
static_cast<unsigned int>(SurfaceID::ReceiverSide));
cell->face(f)->set_user_index(
static_cast<unsigned int>(SurfaceID::SourceSide) &&
static_cast<unsigned int>(MaterialID::Concrete)));
if (cell->face(f)->at_boundary())
if ((p[0] + geom_tol) < 203. && p[0] > geom_tol)
cell->face(f)->set_boundary_id(
static_cast<unsigned int>(SurfaceID::FixedBoundary));
if (
std::abs(p[0] - 1312.) < geom_tol &&
static_cast<unsigned int>(MaterialID::Air))
cell->face(f)->set_boundary_id(
static_cast<unsigned int>(SurfaceID::ZeroPressure));
* * for(const auto &cell :triangulation.active_cell_iterators())
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
The wall is fixed at its outer boundary
static_cast<unsigned int>(SurfaceID::FixedBoundary),
component_mask_displacement);
system_matrix.reinit(locally_owned_dofs,
HarmonicResponse<dim>::assemble_system()
std::vector< bool > component_mask
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
FE values for volume integration
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
Common face quadrature
const QGauss<dim - 1> common_face_quadrature(fe_collection.max_degree() +
FE face values for surface integration
Local element matrix for a cell
Local interface matrix between air and structure DoFs
fe_air.n_dofs_per_cell(), fe_structure.n_dofs_per_cell());
Right-hand side
std::vector<types::global_dof_index> local_dof_indices;
std::vector<types::global_dof_index> neighbor_dof_indices;
for (
const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
Assemble air cells contributions
if (cell->material_id() ==
static_cast<unsigned int>(MaterialID::Air))
hp_fe_values.reinit(cell);
const unsigned int dofs_per_cell = fe_air.n_dofs_per_cell();
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
cell_rhs.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
neighbor_dof_indices.resize(fe_structure.n_dofs_per_cell());
const double sound_speed = get_sound_speed();
const unsigned int n_quadrature_points
const FEValues< dim, spacedim > & get_present_fe_values() const
Wave number k
const std::complex<double> k = (omega / sound_speed);
for (
unsigned int q = 0; q < n_q_points; ++q)
const auto JxW = fe_values.JxW(q);
const Point<dim> &q_point = fe_values.quadrature_point(q);
calucalte lambda and J for PML
pml.vector_value(q_point, lambda);
const std::complex<double> J =
lambda[0] * lambda[1] * lambda[2];
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
const double phi_i = fe_values[pressure].value(i, q);
Gradient of phi_i, which is multiplied by 1/lambda
fe_values[pressure].gradient(i, q);
for (
unsigned int d = 0; d < dim; ++d)
phi_i_i[d] *= 1.0 / lambda[d];
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
fe_values[pressure].value(j, q);
Gradient of phi_j, which is multiplied by 1/lambda
fe_values[pressure].gradient(j, q);
for (
unsigned int d = 0; d < dim; ++d)
phi_j_j[d] *= 1.0 / lambda[d];
cell_matrix[i][j] += phi_i_i * phi_j_j * JxW * J;
cell_matrix[i][j] -= Utilities::fixed_power<2>(k) *
assemble local system to global system.
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix,
Here, the air-structure interface is considered. Similar to step 46, 3 possibilities exist: The neighbor is at the same refinement level and has no children, the neighbor has children and the neighbor is coarser.
for (
const auto f : cell->face_indices())
if (!cell->at_boundary(f))
const auto neighbor = cell->neighbor(f);
if (neighbor->material_id() ==
static_cast<unsigned int>(MaterialID::Concrete))
if (neighbor->center()[0]<203.)
The neighbor is at the same refinement level and has no children.
if ((cell->neighbor(f)->level() == cell->level()) &&
(cell->neighbor(f)->has_children() ==
false))
air_fe_face_values.reinit(cell, f);
elasticity_fe_face_values.reinit(
cell->neighbor_of_neighbor(f));
assemble_air_structure_interface_term(
elasticity_fe_face_values,
cell->neighbor(f)->get_dof_indices(
constraints.distribute_local_to_global(
The neighbor has children.
else if ((cell->neighbor(f)->level() ==
(cell->neighbor(f)->has_children() ==
for (
unsigned int subface = 0;
subface < cell->face(f)->n_children();
air_fe_sub_face_values.reinit(cell,
elasticity_fe_face_values.reinit(
cell->neighbor_child_on_subface(f,
cell->neighbor_of_neighbor(f));
assemble_air_structure_interface_term(
elasticity_fe_face_values,
cell->neighbor_child_on_subface(f, subface)
->get_dof_indices(neighbor_dof_indices);
constraints.distribute_local_to_global(
The neighbor is coarser.
else if (cell->neighbor_is_coarser(f))
air_fe_face_values.reinit(cell, f);
elasticity_fe_sub_face_values.reinit(
cell->neighbor_of_coarser_neighbor(f).first,
cell->neighbor_of_coarser_neighbor(f).second);
assemble_air_structure_interface_term(
elasticity_fe_sub_face_values,
cell->neighbor(f)->get_dof_indices(
constraints.distribute_local_to_global(
Assemble structure cells contributions
else if (cell->is_locally_owned() &&
static_cast<unsigned int>(MaterialID::Concrete))
hp_fe_values.reinit(cell);
const unsigned int dofs_per_cell =
fe_structure.n_dofs_per_cell();
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
cell_rhs.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
neighbor_dof_indices.resize(fe_air.n_dofs_per_cell());
get_stiffness_tensor<dim>();
const double density = get_density_structure();
const std::complex<double> iso_loss = get_iso_loss();
for (
unsigned int q = 0; q < n_q_points; ++q)
const auto JxW = fe_values.JxW(q);
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
fe_values[displacement].value(i, q);
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
fe_values[displacement].value(j, q);
get_strain(fe_values, i, q);
get_strain(fe_values, j, q);
cell_matrix[i][j] += eps_phi_i * stiffness_tensor *
iso_loss * eps_phi_j * JxW;
Utilities::fixed_power<2>(omega) * phi_j * phi_i *
Here the diffuse sound field is considered on the source side of the wall.
for (
const auto face_no : cell->face_indices())
if (cell->face(face_no)->user_index() ==
static_cast<unsigned
int>(SurfaceID::SourceSide))
hp_fe_face_values.reinit(cell, face_no);
const unsigned int n_face_q_points =
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
if (fe_structure.has_support_on_face(i, face_no))
std::vector<std::complex<double>> pressure(
const auto face_quadrature_points =
fe_face_values.get_quadrature_points();
field.value_list(face_quadrature_points,
for (
unsigned int q = 0; q < n_face_q_points;
const auto JxW = fe_face_values.JxW(q);
fe_face_values[displacement].value(i, q);
-fe_face_values.normal_vector(q);
cell_rhs[i] += pressure[q] *
(NormalVector * phi_i) * JxW;
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix,
HarmonicResponse<dim>::assemble_air_structure_interface_term(
FullMatrix<std::complex<double>> &local_interface_matrix,
{
local_interface_matrix = 0;
const auto density_air = get_density_air();
const unsigned int n_face_quadrature_points =
air_fe_face_values.n_quadrature_points;
for (
unsigned int q = 0; q < n_face_quadrature_points; ++q)
-air_fe_face_values.normal_vector(q);
for (
unsigned int i = 0; i < air_fe_face_values.dofs_per_cell; ++i)
const double phi_i = air_fe_face_values[pressure].value(i, q);
j < elasticity_fe_face_values.dofs_per_cell;
elasticity_fe_face_values[displacement].value(j, q);
local_interface_matrix[i][j] +=
density_air * Utilities::fixed_power<2>(omega) * phi_i *
normalVectorStructure * phi_j * air_fe_face_values.JxW(q);
local_interface_matrix[i][j] +=
phi_j * normalVectorStructure * phi_i *
elasticity_fe_face_values.JxW(q);
const FEFaceValues< dim, spacedim > & get_present_fe_values() const
calculation of incident sound power on the source side
HarmonicResponse<dim>::incident_sound_power()
face_quadrature_formula_structure,
std::vector<types::global_dof_index> local_dof_indices;
const double density_air = get_density_air();
for (
const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (const auto face_no : cell->face_indices())
if (cell->face(face_no)->at_boundary() &&
(cell->face(face_no)->user_index() ==
static_cast<unsigned
int>(SurfaceID::SourceSide)) &&
static_cast<unsigned
int>(MaterialID::Concrete)))
structure_fe_face_values.reinit(cell, face_no);
const unsigned int n_face_q_points =
structure_fe_face_values.n_quadrature_points;
std::vector<std::complex<double>> pressure(n_face_q_points);
const auto face_quadrature_points =
structure_fe_face_values.get_quadrature_points();
field.value_list_incidence(face_quadrature_points,
std::vector<Tensor<1, dim, std::complex<double>>>
pressure_gradients(n_face_q_points);
field.gradient_list(face_quadrature_points,
for (
unsigned int q = 0; q < n_face_q_points; ++q)
-structure_fe_face_values.normal_vector(q);
const auto JxW = structure_fe_face_values.JxW(q);
std::complex<double> v_n =
(-1.0 / (omega * std::complex<double>(0., 1.) *
(pressure_gradients[q] * NormalVector);
0.5 * std::real(pressure[q] * std::conj(v_n) * JxW);
Calculating the sound power on the receiver side
HarmonicResponse<dim>::receiver_sound_power()
Common face quadrature
const QGauss<dim - 1> common_face_quadrature(fe_collection.max_degree() +
FE face values for surface integration
for (
const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
for (const auto f : cell->face_indices())
if ((cell->face(f)->user_index() ==
static_cast<unsigned
int>(SurfaceID::ReceiverSide)) &&
static_cast<unsigned
int>(MaterialID::Concrete)))
const auto neighbor = cell->neighbor(f);
The neighbor is at the same refinement level and has no children.
if ((neighbor->level() == cell->level()) &&
(neighbor->has_children() ==
false) &&
(neighbor->material_id() ==
static_cast<unsigned int>(MaterialID::Air)))
elasticity_fe_face_values.reinit(cell, f);
air_fe_face_values.reinit(
neighbor, cell->neighbor_of_neighbor(f));
cell_receiver_sound_power(elasticity_fe_face_values,
The neighbor has children.
else if ((neighbor->level() == cell->level()) &&
(neighbor->has_children() ==
true))
for (
unsigned int subface = 0;
subface < cell->face(f)->n_children();
elasticity_fe_sub_face_values.reinit(cell,
air_fe_face_values.reinit(
cell->neighbor_child_on_subface(f, subface),
cell->neighbor_of_neighbor(f));
sound_power += cell_receiver_sound_power(
elasticity_fe_sub_face_values,
The neighbor is coarser.
else if (cell->neighbor_is_coarser(f))
elasticity_fe_face_values.reinit(cell, f);
air_fe_sub_face_values.reinit(
cell->neighbor_of_coarser_neighbor(f).first,
cell->neighbor_of_coarser_neighbor(f).second);
cell_receiver_sound_power(elasticity_fe_face_values,
HarmonicResponse<dim>::cell_receiver_sound_power(
{
std::vector<types::global_dof_index> local_dof_indices;
const unsigned int n_q_points = air_fe_face_values.n_quadrature_points;
std::vector<std::complex<double>> local_dof_values_pressure(
air_fe_face_values.n_quadrature_points);
air_fe_face_values[pressure].get_function_values(locally_relevant_solution,
local_dof_values_pressure);
std::vector<Tensor<1, dim, std::complex<double>>>
local_dof_values_displacement(n_q_points);
elasticity_fe_face_values[displacement].get_function_values(
locally_relevant_solution, local_dof_values_displacement);
for (
unsigned int q = 0; q < n_q_points; ++q)
air_fe_face_values.normal_vector(q);
const auto JxW = air_fe_face_values.JxW(q);
std::complex<double> normal_velocity =
local_dof_values_displacement[q] * NormalVector *
std::complex<double>(0., 1.) * omega;
std::real(local_dof_values_pressure[q] * std::conj(normal_velocity)) *
HarmonicResponse<dim>::solve()
locally_owned_dofs, mpi_communicator);
solver.solve(system_matrix, completely_distributed_solution, system_rhs);
constraints.distribute(completely_distributed_solution);
locally_relevant_solution = completely_distributed_solution;
HarmonicResponse<dim>::calculate_magnitude()
for (
const auto i : locally_owned_dofs)
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
For postprocessing, the real part represents the magnitude, while the imaginary part vanishes.
locally_relevant_magnitude[i] =
std::abs(value);
HarmonicResponse<dim>::run(
bool write_output)
{
static bool mode_output =
false;
pcout <<
"Debug mode" << std::endl;
pcout <<
"Release mode" << std::endl;
pcout << std::endl <<
"Frequency = " << frequency <<
" Hz" << std::endl;
std::ifstream input_file(
"./tria.inp");
grid_in.read_abaqus(input_file);
triangulation.refine_global(1);
for (
const auto &cell : triangulation.active_cell_iterators())
if (cell->center()[0] > 203.)
triangulation.execute_coarsening_and_refinement();
pcout <<
" setup_system" << std::endl;
pcout <<
" assemble_system" << std::endl;
pcout <<
" solve" << std::endl;
void attach_triangulation(Triangulation< dim, spacedim > &tria)
Calculation of the magnitude of the solution. For postprocessing, the magnitude of the solutions corresponds to the real part and the imaginary part vaishes.
const auto sound_power_source_side_local = incident_sound_power();
const auto sound_power_receiver_side_local = receiver_sound_power();
const auto sound_power_source_side =
::Utilities::MPI::sum(sound_power_source_side_local,
const auto sound_power_receiver_side =
::Utilities::MPI::sum(sound_power_receiver_side_local,
10. * std::log10(sound_power_source_side / sound_power_receiver_side);
pcout <<
" STL = " << stl <<
" dB" << std::endl;
std::ofstream outfile(
"STL_result.txt", std::ios::app);
outfile << frequency <<
" " << stl << std::endl;
std::vector<std::string> solution_names_magnitude(dim,
"magnitude_u");
solution_names_magnitude.push_back(
"magnitude_p");
std::vector<std::string> solution_names(dim,
"u");
solution_names.push_back(
"p");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_out.add_data_vector(locally_relevant_solution,
data_out.add_data_vector(locally_relevant_magnitude,
solution_names_magnitude,
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
@ component_is_part_of_vector
subdomain visualization (unchanged)
for (
unsigned int i = 0; i < subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector(subdomain,
"subdomain");
data_out.build_patches();
ss << std::fixed << std::setprecision(2) << frequency;
std::filesystem::create_directories(
"./output/");
data_out.write_vtu_with_pvtu_record(
"./output/",
"solution_" + ss.str(), 0, mpi_communicator, 2, 8);
main(
int argc,
char *argv[])
{
const unsigned int dim = 3;
f_c = f0 *
std::pow(10.,
static_cast<double>(n) / 60.);
VibroAcousticProblem::HarmonicResponse<dim> elastic_problem(omega);
elastic_problem.run(
true);
catch (std::exception &exc)
<<
"----------------------------------------------------"
std::cerr <<
"Exception on processing: " << std::endl
<< exc.what() << std::endl
<<
"Aborting!" << std::endl
<<
"----------------------------------------------------"
<<
"----------------------------------------------------"
std::cerr <<
"Unknown exception!" << std::endl
<<
"Aborting!" << std::endl
<<
"----------------------------------------------------"
* * int main(int argc, char **argv)
1.9.8