In this document we discuss the solution of time-harmonic acoustic fluid-structure interaction problems. We start by reviewing the relevant theory and then present the solution of a simple model problem – the sound radiation from an oscillating circular cylinder that is coated with a compressible elastic layer.
This problem combines the problems discussed in the tutorials illustrating
The figure below shows a sketch of a representative model problem: a circular cylinder is immersed in an inviscid compressible fluid and performs a prescribed two-dimensional harmonic oscillation of radian frequency . The cylinder is coated with a compressible elastic layer. We wish to compute the displacement field in the elastic coating (assumed to be described by the equations of time-harmonic linear elasticity) and the pressure distribution in the fluid (governed by the Helmholtz equation). The two sets of equations interact at the interface between fluid and solid: the fluid pressure exerts a traction onto the elastic layer, while the motion of the elastic layer drives the fluid motion via the non-penetration condition.
We describe the behaviour of the fluid in terms of the displacement field, , of the fluid particles. As usual we employ index notation and the summation convention, and use asterisks to distinguish dimensional quantities from their non-dimensional equivalents. The fluid is inviscid and compressible, with a bulk modulus , such that the acoustic pressure is given by We assume that the fluid motion is irrotational and can be described by a displacement potential , such that We consider steady-state time-harmonic oscillations and write the displacement potential and the pressure as and , respectively, where denotes the real part. For small disturbances, the linearised Euler equation reveals that the time-harmonic pressure is related to the displacement potential via where is the ambient fluid density. We non-dimensionalise all lengths on a problem-specific lengthscale (e.g. the outer radius of the coating layer) such that and . The non-dimensional displacement potential is then governed by the Helmholtz equation
where the square of the non-dimensional wavenumber,
represents the ratio of the typical inertial fluid pressure induced by the wall oscillation to the `stiffness' of the fluid.
We model the coating layer as a linearly elastic solid, described in terms of a displacement field , with stress tensor
where and are the material's Young's modulus and Poisson's ratio, respectively. As before, we assume a time-harmonic solution with frequency so that , and we non-dimensionalise the displacements on and the stress on Young's modulus, , so that and . The deformation of the elastic coating is then governed by the time-harmonic Navier-Lame equations
which depend (implicitly) on Poisson's ratio , and on the (square of the) non-dimensional wavenumber
where is the solid density. The parameter represents the ratio of the typical inertial solid pressure induced by the wall oscillation to the stiffness of the elastic coating. We note that for a `light' coating we have .
The inner surface of the elastic coating, , is subject to the prescribed displacement imposed by the oscillating cylinder. For instance, if the inner cylinder performs axisymmetric oscillations of non-dimensional amplitude , we have
where is the unit vector in the radial direction. The fluid-loaded surface of the elastic coating, , is subject to the fluid pressure. The non-dimensional traction exerted by the fluid onto the solid (on the solid stress scale) is therefore given by
where the are the components of the outer unit normal on the solid boundary and
is the final non-dimensional parameter in the problem. It represents the ratio of the typical inertial fluid pressure induced by the wall oscillation to the stiffness of the elastic coating. The parameter therefore provides a measure of the strength of the fluid-structure interaction (FSI) in the sense that for the elastic coating does not `feel' the presence of the fluid.
The fluid is forced by the normal displacement of the solid. Imposing the non-penetration condition on yields a Neumann condition for the displacement potential,
Finally, the displacement potential for the fluid must satisfy the Sommerfeld radiation condition
which ensures that the oscillating cylinder does not generate any incoming waves.
The implementation of the coupled problem follows the usual procedure for multi-domain problems in
oomph-lib. We discretise the constituent single-physics problems using the existing single-physics elements, here
for the discretisation of the PDEs (1) and (2), respectively. The displacement boundary condition (3) on the inner surface of the elastic coating is imposed as usual by pinning the relevant degrees of freedom, exactly as in a single-physics solid mechanics problem. Similarly, the Sommerfeld radiation condition (6) on the outer boundary of the fluid domain can be imposed by any of the methods available for the solution of the single-physics Helmholtz equation, such as approximate/absorbing boundary conditions (ABCs) or a Dirichlet-to-Neumann mapping.
The boundary conditions (4) and (5) at the fluid-solid interface are traction boundary conditions for the solid, and Neumann boundary conditions for the Helmholtz equation, respectively. In a single-physics problem we would impose such boundary conditions by attaching suitable
FaceElements to the appropriate boundaries of the "bulk" elements, as shown in the sketch below:
TimeHarmonicLinearElasticityTractionElements could be used to impose a (given) traction, , onto the solid;
HelmholtzFluxElements could be used to impose a (given) normal derivative, , on the displacement potential. Both and would usually be specified in a user-defined namespace and accessed via function pointers as indicated in the right half of the sketch.
In the coupled problem, illustrated in the left half of the next sketch, the traction acting on the solid becomes a function of the displacement potential via the boundary condition (4), while the normal derivative of the displacement potential is given in terms of the solid displacement via equation (5). Note that corresponding points on the FSI boundary are identified by matching values of the boundary coordinate which is assumed to be consistent between the two domains.
The implementation of this interaction in the discretised problem is illustrated in the right half of the sketch: We replace the single-physics
HelmholtzFluxFromNormalDisplacementBCElements, and the
TimeHarmonicLinElastLoadedByHelmholtzPressureBCElements. (Yes, we like to be verbose...). Both of these
FaceElements are derived from the
ElementWithExternalElement base class and can therefore store a pointer to an "external" element that provides the information required to impose the appropriate boundary condition. Thus, the
HelmholtzFluxFromNormalDisplacementBCElements store pointers to the "adjacent" time-harmonic linear elasticity elements (from which they obtain the boundary displacement required for the imposition of (5)), while the
TimeHarmonicLinElastLoadedByHelmholtzPressureBCElements store pointers to the "adjacent" Helmholtz elements that provide the value of the displacement potential required for the evaluation of (4).
The identification of the "adjacent" bulk elements can be performed using the
Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh(...) helper function. We note that, as suggested by the sketch above, this function does not require to the two adjacent meshes to have a consistent discretisation – the identification of adjacent elements is based entirely on the (assumed to be consistent) boundary coordinate in the two meshes. We refer to another tutorial for a discussion of how to set up (or change) the parametrisation of mesh boundaries by boundary coordinates.
The animation below shows the deformation of the elastic coating if a non-axisymmetric displacement
(for ) is imposed on the inner boundary of the coating .
Here is a plot of the corresponding pressure field:
Finally, we provide some validation of the computational results by comparing the non-dimensional time-average radiated power
against the analytical solution for axisymmetric forcing ( ) for the parameter values , , and a non-dimensional coating thickness of see Heil, M., Kharrat, T., Cotterill, P.A. & Abrahams, I.D. (2012) Quasi-resonances in sound-insulating coatings. Journal of Sound and Vibration 331 4774-4784 for details. In the computations the integral in (8) is evaluated along the outer boundary of the computational domain.
As usual we define the problem parameters in a namespace.
We wish to perform parameter studies in which we vary the FSI parameter . To make this physically meaningful, we interpret as a measure of the stiffness of the elastic coating (so that an increase in corresponds to a reduction in the layer's elastic modulus ). In that case, the frequency parameter in the time-harmonic linear elasticity equations becomes a dependent parameter and is given in terms of the density ratio and by . We therefore provide a helper function to update the dependent parameter following any change in the independent parameters.
We force the system by imposing a prescribed displacement on the inner surface of the elastic coating and allow this to vary in the azimuthal direction with wavenumber :
The rest of the namespace contains lengthy expressions for various exact solutions and is omitted here.
The driver code is very straightforward. We parse the command line to determine the parameters for the parameter study and build the problem object, using refineable nine-noded quadrilateral elements for the solution of the time-harmonic elasticity and Helmholtz equations.
We then solve the problem for various values of , updating the dependent variables after every increment.
Problem class is templated by the types of the "bulk" elements used to discretise the time-harmonic linear elasticity and Helmholtz equations, respectively. It contains the usual member functions to detach and attach
FaceElements from the bulk meshes before and after any mesh adaptation, respectively.
The private member data includes storage for the various meshes and objects that are used for outputting the results.
We start by building the meshes for the elasticity and Helmholtz equations. Both domains are complete annular regions, so the annular mesh (which is built from a rectangular quad mesh) is periodic.
The solid mesh occupies the region between and where is the thickness of the elastic coating:
The Helmholtz mesh occupies the region between and where is the outer radius of the computational domain where we will apply the Sommerfeld radiation condition. Note that the two meshes are not matching – both meshes have 3 element layers in the radial direction but 10 and 11 in the azimuthal direction, respectively. This is done mainly to illustrate our claim that the multi-domain setup functions can operate with non-matching meshes.
Both bulk meshes are adaptive so we create error estimators for them:
Next we create the mesh that will store the
FaceElements that will apply the Sommerfeld radiation condition, using the specified number of Fourier terms in the Dirichlet-to-Neumann mapping; see the Helmholtz tutorial for details.
Next we pass the problem parameters to the bulk elements. The elasticity elements require a pointer to the elasticity tensor and the frequency parameter :
The Helmholtz elements need a pointer to the (square of the) wavenumber, :
It is always a good idea to check the enumeration of the mesh boundaries to facilitate the application of boundary conditions:
Next we create the meshes containing the various
FaceElements used to apply to the FSI traction boundary condition (4), the FSI flux boundary condition (5) for the Helmholtz equation, and the Sommerfeld radiation condition (6), respectively, using helper functions discussed below.
We add the various sub-meshes to the problem and build the global mesh
The solid displacements are prescribed on the inner boundary (boundary 0) of the solid mesh so we pin all four values (representing the real and imaginary parts of the displacements in the and directions, respectively) and assign the boundary values using the function
Global_Parameters::solid_boundary_displacement(...). (The enumeration of the unknowns is discussed in another tutorial.)
Finally, we set up the fluid-structure interaction, assign the equation numbers, define the output directory and open a trace file to record the radiated power as a function of the FSI parameter .
The mesh adaptation is driven by the error estimates for the bulk elements. The various
FaceElements must therefore be removed from the global mesh before the adaptation takes place. We do this by calling the helper function
delete_face_elements(...) (discussed below) for the three face meshes, before rebuilding the Problem's global mesh.
After the (bulk-)mesh has been adapted, the various
FaceElements must be re-attached. We then (re-)setup the fluid-structure interaction and rebuild the global mesh.
The helper function
delete_face_elements() is used to delete all
FaceElements in a given surface mesh before the mesh adaptation.
create_fsi_traction_elements() creates the
FaceElements required to apply the FSI traction boundary condition (4) on the outer boundary (boundary 2) of the solid mesh:
To function properly, the elements need to know the number of the bulk mesh boundary they are attached to (this allows them to determine the boundary coordinate required to set up the fluid-structure interaction; see Implementation ), and the FSI parameter .
[Note: We omit the listings of the functions
create_helmholtz_DtN_elements() which create the
FaceElements required to apply the FSI flux boundary condition (5) on the inner boundary (boundary 0), and the Sommerfeld radiation condition (6) on the outer boundary (boundary 2) of the Helmholtz mesh because they are very similar. Feel free to inspect the source code.]
The setup of the fluid-structure interaction requires the identification of the "bulk" Helmholtz elements that are adjacent to (the Gauss points of) the
FaceElements that impose the FSI traction boundary condition (4), in terms of the displacement potential computed by these "bulk" elements. This can be done using the helper function
Multi_domain_functions::setup_bulk_elements_adjacent_to_face_mesh(...) which is templated by the type of the "bulk" element and its spatial dimension, and takes as arguments:
Nearly a one-liner:
Exactly the same method can be used for the identification of the "bulk" elasticity elements that are adjacent to (the Gauss points of) the
FaceElements that impose the FSI flux boundary condition (5), using the displacement computed by these "bulk" elements:
The post-processing function
doc_solution(...) computes and outputs the total radiated power, and plots the computed and exact solutions (real and imaginary parts) for all fields.
–el_multiplier(which controls the number of elements in the mesh) and suppressing any automatic (un)refinement by setting the maximum number of adaptations to zero using the
–max_adaptcommand line flag.
RefineableQHelmholtzElement<2,2>to the bi-cubic
A pdf version of this document is available.