Using BEMPP to Calculate Far-Field Sound Pressure Around a Hollow Cube with Imported Surface Velocities

Hello Everone,

I’m working on a student project where I need to use BEMPP to calculate the sound pressure in the far-field around a hollow cube. The cube is excited by a steady-state load, and I’ve conducted a mode-based steady-state dynamic analysis in Abaqus to obtain the surface velocities of the cube.

I have the normal velocities of all surface elements from Abaqus at various frequencies, stored in an Excel file. I would like to know how to use this velocity data in BEMPP to calculate the sound pressure around the cube using far-field boundary conditions. My goal is to compare the sound pressure results from BEMPP with those obtained from an acoustic cavity mesh analysis performed in Abaqus.

I have reviewed Tutorial 5: Loudspeaker Radiation from the BEMPP tutorials, where the cone velocity is set as a constant (0.001 m/s). However, in my case, I have velocity data that is not constant for the cube’s surface. Could anybody please guide how to incorporate this data into BEMPP?

Additionally, is it possible to import a mesh created in Hypermesh into BEMPP? If so, what format is needed, and are any specific steps required to ensure compatibility?

Thank you for your assistance!

Hi Kevin,

The way I compute frequency responses with bempp-cl is by looping through the frequency array and radiating surfaces of my system. I don’t know if there is a better way to do it (which doesn’t involve loops).

import numpy
import bempp.api
from bempp.api.linalg import gmres
from bempp.api.operators.boundary import helmholtz, sparse


#%% mesh
# grid = however you import your grid

#%% study parameters
frequency = np.logspace(1, 3, 50)
c         = 343                   # speed of sound in air
rho       = 1.22                  # air density
omega     = 2 * np.pi * frequency # angular frequency
k         = omega / c             # wavenumber
surface   = [1, 2, 3, 4, 5, 6]    # your cube surface groups
nSurf     = len(surfaces)
nFreq     = len(frequency)


#%% spaceP + identity definition
spaceP = bempp.api.function_space(grid, "P", 1)
I      = sparse.identity(spaceP, spaceP, spaceP)


#%% initialize array to store mesh pressure and mesh velocity
p_mesh = np.empty([nFreq, nSurf], dtype=object)
u_mesh = np.empty([nFreq, nSurf], dtype=object)

#%% Some loops (maybe this can be optimized/replaced?)
for s in range(nSurf):     # surfaces
    spaceU = bempp.api.function_space(grid, "DP", 0, segments=[surface[s]]])
    coeff_surface = np.ones(spaceU.grid_dof_count)
    for f in range(nFreq): # frequencies
        D      = helmhotz.double_layer(spaceP, spaceP, spaceP, -k[f])
        S      = helmhotz.single_layer(spaceU, spaceP, spaceP, -k[f])
        u_surf = bempp.api.GridFunction(spaceU, coefficients=coeff_surface)
        lhs    = D - 0.5 * I  # change - to + for interior studies
        rhs    = 1j * omega[f] * rho * S * u_surf

        # compute + store results
        p_mesh[f, s], _ = gmres(lhs, rhs)
        u_mesh[f, s] = u_surf


# Apply a similar approach to calculate pressure at evaluation points.

Few points about the code above:

  • it uses “-k” to change the sign of the Green function (-jkr instead of +jkr): I had “wrong” results using the +jkr one – but idk if this is an error in the way I implement computations or not,
  • gmres results are stored in an empty array of objects, so they can easily be used for calculation of potentials (didn’t find another way to “store” GridFunctions),
  • Unit velocity is set on each radiating surface: the velocity transfer functions are applied on final results. Very handy if you want to change the velocities without running the whole BEM computation again.

I have no idea about the Hypermesh thing. I believe bempp-cl uses meshio, so my guess is any format compatible with meshio is also compatible with bempp-cl.

Cheers

1 Like

Dear Tchoum,

I am extremely thankful for your answer. It has provided me with a much better understanding of the process.

I would be grateful if you could guide me on how to apply velocity transfer functions to the final results.

As part of my project, I am starting with a simple cube model to familiarize myself with BEMPP before moving on to a gearbox casing model that contains around 200,000 nodes. For each of these nodes, I have the normal velocities and phase information at different frequencies, stored in a CSV format. I have managed to load this data into a pandas DataFrame in Python, but I’m unsure how to apply the velocities to the corresponding nodes in BEMPP. While I have the node numbers listed in one of the columns, I’m not sure how to map these to the respective nodes within BEMPP.

Your insights on how to approach this would be extremely helpful.

Thank you very much once again!

For importing velocity transfer function, I would use numpy.loadtxt(), as it returns an array that you can easily manipulate, slice, or multiply with other numpy arrays. I’m not familiar with Pandas, so it may also be capable of handling this task as well.

To apply velocity transfer functions onto results:

#%% First, compute potentials and get the results into a manageable numpy array
micPosition = np.array([[xmic1, ymic1, zmic1],
                        [xmic2, ymic2, zmic2],
                        [..., ..., ...]]).T
nMic = np.shape(micPosition)[1]

# to store projected pressure
pressure   = np.zeros([nFreq, nMic, nSurf], dtype=complex)  # store results

# pretty sure you can regroup this with the previous loops
for f in range(nFreq):       # looping through frequencies
    for s in range(nSurf):   # looping through surfaces
        spaceU = bempp.api.function_space(grid, "DP", 0, segments=[surface[s]]])
        pressure[f, :, s] = np.reshape(
            helmholtz_potential.double_layer(spaceP, micPosition, -k[f]) *
            p_mesh[f, s] - 1j * omega[f] * rho *
            helmholtz_potential.single_layer(spaceU, micPosition, -k[f]) *
            u_mesh[f, s], nMic)


#%% Apply velocity to results (complex valued)
pressure_p = np.copy(pressure)  # to keep initial results from changing

# velocity of first surface
pressure_p[:, :, 0] *= your_velocity_function  # 1D (velocity) array might need to be tiled or repeated
                                               # into a 2D array

# ...do the same for all surfaces

# sum the total contribution from all faces -> output array of shape (nFreq, nMic)
pressure_s = np.sum(pressure_p, axis=2)

However, this example assumes that each radiating surface has only one uniform velocity defined across its boundary: treating it as an infinitely rigid piston with no modal behavior.

In your case, where you want to apply independent velocities to each node, this approach won’t work. Instead of using a single velocity for an entire surface, you’ll need to apply velocity independently at each degree of freedom (DOF). In that case, you can replace the previous “coeff_surface” argument to allow for independent velocities at each DOF. It would be tricky to explain fully here. But depending on your mesh/application, you can define a function_space that already incorporates your Abaqus nodes as DOFs, which will allow you to apply the velocities directly in coeff_surface – as long as your grid indexing corresponds to the DOF indexing.

This whole thing can be a bit difficult to understand, but I hope it helps!

1 Like

Dear Tchoum,

Again, thank you so much for taking the time to answer my questions.

As you said my grid indexing matched to the DOF indexing. So I applied the velocities directly using the command: u_total = bempp.api.GridFunction(spaceU, coefficients=complex_velocity),
where complex_velocity is a 1D numpy array that contains the velocity (complex value) of all the nodes. The code ran without any errors for a test run for a single frequency. I loaded the velocity as a grid function before solving the linear system of equations.

It would be really helpful if you could answer these questions:

  1. Since I needed to assign velocities to each node, I had to use the function space ‘P’ of continuous polynomials (spaceU = bempp.api.function_space(grid, “P”, 1) ) for velocities. But the velocities are not continuous across nodes. The code ran without any error but I am not sure if this is the correct procedure. Could you please share your insights regarding this?

  2. I did not understand why did you create a loop through the radiating surfaces in your code? Once we create a SpaceP, SpaceU and double_layer, single_layer using these spaces, we can solve this equation to get the pressure values :
    p_total, info = gmres(double_layer - 0.5 * identity,
    1j * omega * rho_0 * single_layer * u_total, tol=1E-5).

Would this approach give the same answer? I couldn’t create a loop over the surface because my mesh created in hypermesh doesn’t have domain_index info to separate the surfaces, but I can still use it in BEMPP because meshio supports hypermesh files.

Hi Kevin,

Regarding the continuous polynomials function space: To be honest, I am not sure to fully understand the mathematics behind BEM. But in my experience, while results where not exactly the same, both continuous and discontinuous polynomials gave quite similar results when studying membrane deformation of a loudspeaker (based on measurement from a laser vibrometer), and both were close to actual anechoic measurements. If you take a look at the loudspeaker tutorial, and compare P1 with DP0, you get (in real/imag):

P1_DP0

If you look at the dB difference (outside the speaker boundaries), it is minimal:
P1_DP0_dB_diff

Now, is it the right approach? It should be worth trying the discontinuous polynomial as well. In that case, you “just” need to interpolate your measurements or simulation data to the element centroids.

Regarding radiating-surface loop: As I understand your system, you won’t need it, I was taking a general example. For instance, on a loudspeaker study with a membrane (see next fig, in red) and port (blue) — both have a their own velocity, so I use two spaceU and u_total to separate each contribution:

ported_box

Here you can see the pressure radiated by both driver and port:
spk_radiation

It really is the only reason why I am looping through surfaces.

1 Like

Thank you so much for your help!
I used a continuous field for velocity. I calculated SPL(DB) results from both Abaqus and BEMPP for a frequency of 500hz on a sphere of radius 150mm from the center of the sound source. The Y-axis of the plot is SPL(DB) and X-axis is the node index of the sphere. It’s not a great match. I’ll keep trying different models.

Update: For the cube, I considered the surface velocities of just one face for both Abaqus (FEM) and BEMPP and got comparable results

That’s a great starting point! The full cube should work properly (and maybe try adding faces gradually). If it doesn’t, you might be having an issue with the Green’s function. In that case, try switching -k to +k (or vice versa, depending on what you’re currently using) in your boundary and potential calculations. If it fixes the issue, apply the same stuff for your gearbox study.

1 Like

Hi Tchoum,

Thank you so much for helping me.

When we apply velocities to the scalar continuous polynomial function space (P) and then load the nodal velocities by the command below:

spaceU = bempp.api.function_space(grid, “P”, 1)

c = np.array([…]) # These are the coefficients
grid_fun = GridFunction(spaceU, coefficients=c)

Do you know in which direction the velocity will be applied to the nodes? Its a bit confusing since the function space is scalar and nodal velocities in acoustic analysis is a vector in the normal direction.

Hi Kevin, sorry for the late reply,

From my understanding, the velocity is always normal to the elements of your mesh (I might be wrong here). However, you can use the dot product between an element’s normal and the desired radiation direction to weight its impact on the radiation. This won’t change the radiation direction, just reduce the amplitude of elements that don’t have a normal parallel to your acoustic vector (it must be defined between 0 and 1, e.g. (1, 0.23, 0), (0, 1, 0), etc.). In Tutorial 5 they use a decorator to do more or less the same thing:

@bempp.api.complex_callable
def u_total_callable(x, n, domain_index, result):
    if domain_index == 2:
        result[0] = n[2];
    else:
        result[0] = 0

Here, each element’s contribution is scaled by the z-component of its normal vector.