I am working on a problem involving the calculation of sound pressure using the Helmholtz equation and have complex Normal velocity components (X, Y, Z) for each node in my mesh. According to the Bempp documentation, scalar function spaces are recommended for solving Laplace or Helmholtz problems, while vector function spaces are used for Maxwell’s equations.
Given that velocity is inherently a vector, I am seeking advice on how to properly apply these complex vector components within a scalar function space to solve the Helmholtz equation in Bempp for acoustic problems.
I tried formatting this data into a numpy array with a shape of (number of nodes, 3) to apply it like this:
u_total = bempp.api.GridFunction(spaceU, coefficients=complex_velocity_numpy_array)
I’m encountering an issue where the shapes of the left-hand side (LHS) and the right-hand side (RHS) do not match when solving linear system of equation in GMRES
Could anyone provide guidance on how to approach this? Are there specific techniques within Bempp that would facilitate the integration of vector data into a scalar function space for solving Helmholtz problems?
The Helmholtz equation in Bempp is a scalar partial differential equation which can be used to model harmonic acoustic pressure fields. Notice that elastodynamic simulations are currently not supported in Bempp. That is, one can only model longitudinal sound waves (p-waves) and no shear components (transverse s-waves). Hence, the acoustic pressure is a scalar field.
In the case of boundary integral equations, the right-hand-side and solution vectors consist of Dirichlet and Neumann traces of the incident and total pressure fields, respectively. The Dirichlet trace is the field evaluated at the surface and the Neumann trace is the normal derivative at the surface. In other words, the Neumann trace (\nabla p \cdot \mathbf{n}) is related to the normal particle velocity and is a scalar function.
The grid functions in Bempp can be complex-valued.
The Bempp documentation on the grid and grid functions explains how to retrieve the unit normal directions at the surface mesh.
Thank you for answering @Elwin
I am working on a comparative analysis of sound pressure calculations around a vibrating hollow steel cube (25mm length, 1mm thickness) using Abaqus and BEMPP. The measurements are taken at the nodes of a sphere with a 200mm radius centered on the cube. The cube is excited at a steady state of 500 -1000 Hz, and the setup is intended to simulate free-field conditions as shown in the image.
In Abaqus, I’ve modeled free-field conditions by implementing infinite elements at the outer boundary of the sphere, with the region between the sphere and the cube filled with acoustic tetrahedral elements. In BEMPP, I used the null field approach as explained in the loudspeaker tutorial, in which I calculated the p_total on the cube first and then on the plotting grid which is the outer sphere by using these equations.
It’s intriguing that the deviations increase with the number of faces loaded with normal velocities. I am seeking insights into potential reasons for these deviations. Could it be related to how the different methods (FEM in Abaqus and BEM in BEMPP) handle the boundary interactions or something else?
in the equation for calculating pressure on the outer sphere, should I use helmholtz_potential.double_layer or helmholtz_far.double_layer?
Any guidance or suggestions on how to further investigate and mitigate these discrepancies would be very valuable.
It is difficult to pinpoint the reason for the differences between Abaqus and Bempp.
Here, I guess you’ll need bempp.api.operators.potential.helmholtz.double_layer to calculate the pressure field. The function bempp.api.operators.far_field.helmholtz.double_layer uses the far field approximation and calculates \lim_{r \to \infty} r |p(r,\phi,\theta)| for a given angle (\phi,\theta).
Please be aware that if you use piecewise constant (DP0) spaces, you need to specify the normal velocity in each triangle centre. If you use piecewise linear (P1) spaces, you specify the normal velocity in each vertex of the grid. In that case, you’ll have to integrate over test/basis functions that cover different faces of the cube. Hence, be careful how you handle the edges of the cube.
Thank you very much for taking the time to answer my question.
I used P1 linear space and applied velocity to each vertex. I excluded the vertices at the edges so that it is easier to apply global X, Y, and Z velocities as normal velocities for faces normal to the X, Y, and Z axes respectively. I couldn’t understand the “Integrate over basis function that covers different faces of the cube”. Could you kindly provide further details on this point? I would greatly appreciate guidance on where exactly this integration should be implemented in the below-mentioned code. Thank you!
In Abaqus, I used infinite elements as they model the far-field region. However, in Bempp, the results calculated from bempp.api.operators.potential.helmholtz.double_layer are more closely aligned with Abaqus results than those from bempp.api.operators.far_field.helmholtz.double_layer. Could you provide some insight into why this might be the case?
The P1 elements in Bempp implement continuous piecewise-linear functions that have value one in a specific node, value zero in all other nodes, and are linear in each triangle’s interior. See Function Spaces — The Bempp Book for more information and a picture.
When creating a GridFunction, you create a function f(\mathbf{x}) = \sum_{i=1}^N \alpha_i g_i(\mathbf{x}) where g_i(\mathbf{x}) is a basis function and \alpha_i is its coefficient. In your case, g_i is a P1 function corresponding to grid node i and the \alpha_i values come from the complex_velocity array.
Indeed, at a grid node \mathbf{x}_j, we have f(\mathbf{x}_j) = \alpha_j so you can interpret the coefficients as the function value in the nodes. Inside each triangle, the function is a weighted average of three basis functions.
At the edges of a cube, the normal direction does not exist. At the same time, a basis function g_j for a grid node j at the cube’s edge has support on all surrounding triangles. The normal direction exists on each of these triangles, but is different depending on the cube’s face. When specifying an analytical function for the argument fun, Bempp will take this into account. When specifying the coefficients argument, you will have to specify correct values. One idea is to interpret the normal direction at the cube’s edge as an average of the connecting faces, but make sure to normalize the vector to unit length.
Regarding the potential operators: since you evaluate the field relatively near to the object, you should use operators.potential, which calculates the potential operators in the representation formula. The far-field operators in operators.far_field use the far-field approximation of r \to \infty with r the distance to the global origin. This approximation allows to simplify the calculations by eliminating higher-order terms in r. Hence, the far-field operators are less accurate for the near field. Since you evaluate at 200 mm away from a 25 mm cube, I would use the near field operators.
