# laplacian

I am trying to implement the Boudary Layer Impedance model to approximate thermo-viscous losses at solid boundaries for pressure acoustics frequency domain problems, such as described by Bossart et al. (2003) and also by Kampinga (2010).

The idea is to use an impedance boundary condition at solid walls which is first estimated and then updated based on the current solution iteratively until convergence.

In order to calculate the value of the impedance at the boundaries at each iteration, I have to calculate the tangential wavenumber, which is obtained from the tangential laplacian of the pressure at the boundaries.

First question, I tried to calculate the laplacian on the boundaries using:

sol = pde.getSolution()
g_bnd = grad(sol, where = FunctionOnBoundary(domain))
d_bnd = div(g_bnd)

and I get :

Then, assuming I found a solution to that, this is not the value I am looking for because I need the divergence of the components of the gradient parallel to the boundary only. Do I have access to tangential vectors at boundaries?

Many thanks,

Antoine

Rerences
--------------

Bossart, R., Joly, N., & Bruneau, M. (2003). Hybrid numerical and analytical solutions for acoustic boundary problems in thermo-viscous fluids. Journal of Sound and Vibration, 263(1), 69-84.

Kampinga, R. (2010). Viscothermal acoustics using finite elements. Ph.D. thesis. University of Twente, Enschede, The Netherlands.

## Question information

Language:
English Edit question
Status:
For:
esys-escript Edit question
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Lutz Gross (l-gross) said on 2011-12-15: #1

In general escript-finley does not allow 2nd order derivatives due its internal representation of
gradients. you could calcute s=Laplace(sol) by solving the PDE

which makes X= - grad(sol) and y= grad(sol)_i * n_i in the escript PDE . This then can be
interpolated to FunctionOnBoundary(domain)

The calculation of the tangential Laplace will not work with this approach.

 Revision history for this message Antoine Lefebvre (antoine-lefebvre2) said on 2011-12-15: #2

Thanks Lutz,

I found out a way to reformulate the problem using only the gradient and the normal vectors on the boundary.

After performing the calculations, I have a Data object Y, with values for all points on the boundaries but I only need to update the value of the points tagged with a specific value, so I do:

for i in xrange( BC.getNumberOfDataPoints() ):
if BC.getTagNumber(i) == 2:
Yi = Y.getTupleForDataPoint(i)
tmp = wavenumber*Z0*Yi
BC.setValueOfDataPoint(i, [[-tmp,-tmp],[tmp,-tmp]])

pde.setValue(d=BC)

Is this the proper way to update the BC point by point? For a large model, that loop can become expansive, no?

 Revision history for this message Joel Fenwick (j-fenwick1) said on 2011-12-16: #3

Yes the loop would be expensive and it also uses
setValueOfDataPoint()
which we recommend avoiding _if possible_.

I haven't actually tried this but perhaps something like the following would do what you want:

# create a function for all points in Y
tmp=Y*wavenumber*Z0
update=tmp*[[-1, -1], [1, -1]]

# create a function which is 1 for tag 2 and zero for everything else

#selectively copy values into BC where the mask is >0

 Revision history for this message Antoine Lefebvre (antoine-lefebvre2) said on 2011-12-16: #4

Thanks Joel,
The copyWithMask method appears to work properly.

I am still having some issues with that boundary condition I am trying to implement.

I am currently doing:
g_bnd = grad(sol, where = FunctionOnBoundary(domain))

to get the gradient on the boundaries but the component of the gradient normal to the boundary is always 0. Is this the intended behavior?

In my case, I need the normal component. Should I calculate the gradient on Function(domain) instead and use a Projector to get the values on the Face elements?

Thanks again for the help.

 Revision history for this message Antoine Lefebvre (antoine-lefebvre2) said on 2011-12-23: #5

I am trying to calculate the laplacian as suggested by Lutz in answer #1 above but I am not getting correct results. Here is my script with a simple 2d test case. The equation is:

laplacian(p) + k**2*p = 0

with n_i*grad(p)_i = 0 on all boundaries but one where it is set to 1.

-----

cs = 343. # speed of sound (m/sec )
rho = 1.25 # density (kg/m**3)
frequency = 100. # Hz

L = 0.3

wavenumber = 2*pi*frequency/cs

domain = Rectangle(l0=L,l1=0.02,n0=30,n1=10, order=2)

pde = LinearPDE(domain, numEquations=1)

A = pde.createCoefficient("A")
for i in range(domain.getDim()):
A[i,i]=-1

pde.setValue(A=A, D=wavenumber**2)

# normal acceleration BC on input plane
anBC = pde.createCoefficient("y")
anBC.setTaggedValue(1, 1)
pde.setValue(y=anBC)

sol = pde.getSolution()

laplacian = LinearPDE(domain, numEquations=1)

l = laplacian.getSolution()

 Revision history for this message Launchpad Janitor (janitor) said on 2012-01-08: #6

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

 Revision history for this message Lutz Gross (l-gross) said on 2012-01-17: #7

Will take a look at this problem in the next few days.

 Revision history for this message Antoine Lefebvre (antoine-lefebvre2) said on 2012-01-17: #8

Ok. Thanks.

 Revision history for this message Lutz Gross (l-gross) said on 2012-01-20: #9

What is the result for sol you are expecting?

 Revision history for this message Antoine Lefebvre (antoine-lefebvre2) said on 2012-01-20: #10

Do you refer to my original question or to the comment #5 above?

The equation I solved is simply laplacian(p)+k**2*p=0, so when calculating the gradient I should obtain k**2*p as the answer. This is not very useful but I should be able to make it work. This is not what I currently get.

Then, if it works, I'll try to reformulate the method that you proposed in #1 and do a system of three equations to retrieve the 3 spacial components of the vectorial gradient of grad(sol), from which I should be able to calculate the laplacian parallel to the surface and finally go on with implementing the boundary layer impedance model.

 Revision history for this message Lutz Gross (l-gross) said on 2012-01-23: #11

You need to use rich face element when calling "grad(sol,where=FunctionOnBoundary(domain)" and using the component normal to the boundary, see http://esys.esscc.uq.edu.au/esys13/nightly/user/html/node86.html .
Otherwise only the tangential components are calculated properly while the normal component is set to zero - which in you case means that y=0 in the laplacian PDE.

Unfortunately, Rectangle does not use rich face elements by default. To use rich face element you need to call

Rectangle(l0=L,l1=0.02,n0=60,n1=20, order=2, useElementsOnFace=True)

Hope this helps.