[wcpvandervelden]

Dear all,

I am trying to build a solver which solves a non-homogeneous wave equation.

You may ask, why such a solver? It is used for solving a simplified version of the Ffowcs Williams Hawkings differential equation, used in aeroacoustics. To be specific, I would like to use Curle's analogy, and only consider the pressure perturbation at the wall. Therefore I wrote some code which extracts the pressure at the boundary, I multiply it with the normal outward vector, and taking the divergence of the volVectorField (which is zero in the internal field, and has values in the boundary field) to create a volScalarField. Next, I forced this last volScalarField to be zero in the internal field (initial there were some values because of the corners of the geometry)

Now the next step is to solve the wave equation throughout the domain, with the source term as indicated above.

I started by simply using the following relation:

solve((fvm::d2dt2(rhoPrimeCurle) - c0 * c0 * fvm::laplacian(rhoPrimeCurle)== -CurleTerm));

However, due to the fact that the source term only consist of values at the wall, no solution can be obtained. There are 0 iterations required to solve this.

I am thinking this is because of the fact that the internalField = 0. Therefore, I would like to fill the CurleTerm.internalField() somehow by interpolating the values from the boundary to the nearest cells.

Is there anyone who has some experience with this method, or perhaps, is there anyone with a brilliant idea how to omit this. I assume someone tried to implement the FW-H equation before in OF.

Many thanks in advance

[ngj]

Hi

There is a minor detail, which is unclear to me. You say that you generate a surfaceVectorField, where you only populate the boundary values with values different from vector::zero. These values are bc-pressure times the face normal.

Here is the question: How does the divergence of that field end up with a zero internal field? Those cells adjacent to the boundary must have non-zero values, since the boundary faces themselves are non-zero, whereas the internal faces equal vector::zero.

Maybe you should try applying the divergence of your volVectorField in your wave equation; at least this provides you with a non-homogeneous source term.

Kind regards,

Niels

[wcpvandervelden]

Dear Niels,

Thank you very much for your quick reply. I will try to clarify some things now.

First of all, I indeed only copy the values from the wall patch to a new volVectorField like this:

pB.boundaryField()[wallID] = p.boundaryField()[wallID].patchInternalField()*(-mesh.Sf().boundaryField()[wallID]
/mesh.magSf().boundaryField()[wallID]);

where the last term between the ( ) simply indicates the outward normal vector. Now I will have an empty internal field, and only vectors of the boundary patch called "wall" (since wallID finds the patchID "wall").

Next, I calculate the dipole source term for the wave equation:

CurleTerm = rho * fvc::div(pB);

Indeed, as you already mentioned in your post, this results in some nonzero values in this volScalarField "CurleTerm". However, in my opinion these inner values are useless, since they only appear at the corner of meshes, i.e. where the normal outward vector of the adjacent is not simply (1 0 0).

Therefore, afterwards i force the internalField to be zero, i.e.

CurleTerm.internalField() = 0;

I will now end up with a volScalarField, which is completely empty in the internalField but contain values at the boundary patches. This is also what I want, since the source term analytically contains a dirac delta function which is only 1 on the surface.

Equating the source term to the wave equation results in a homogeneous wave equation without source term... Somehow, when solving a PDE, only the values of the internalField are used, and not the values at the boundary patch. Therefore I am now looking to some kind of script which translate the boundary values in the boundary patch, to the adjacent cells in order to solve for the density in the wave equation.

[ngj]

OK. I am not sure that I understand the motivation to clear the divergence field, but never mind. You could e.g. do the following to copy boundary values to internal values:

labelList count( CurleTerm.internalField().size(), 0 );

forAll( mesh.boundary(), patchi )
{
    if ( WALL TYPE)
    {
          const labelList & fc( mesh.boundary()[patchi].faceCells() );
          const scalarField & ctw( CurleTerm.boundaryField()[patchi] );
        
          forAll( ctw, facei )
          {
                CurleTerm[fc[facei]] += ctw[facei];
                count[fc[facei]]++;
          }
    }
}

forAll( count, celli )
{
     if ( count[celli] )
         curleTerm[celli] /= count[celli];
}

This should do what you want. It takes any boundary value and places it in the cell adjacent to the boundary (this is the information contained in faceCells). The counter takes care of corners, e.g. if you have three boundary faces to a corner cell, then you add the values of all boundaries and divide by the number of contributions. I would have liked to avoid the "if", but I could not come up with a clever counting procedure, which avoided dividing by zero in performing the averaging as

curleTerm.internalField /= count;

Best regards,

Niels

[a_habib]

Dear wcpvandervelden,
I saw your work "Numerical investigation of an upper airway in a patient suffering from stridor" and I wounder If you can share your solver, I'm interested in this field.
Regards.