boundary integral

Main Forum Old-Forum - do not post here
1

Hello everyone,

I have a very simple question. I am doing flux calculations and want to calculate line integral on the “outlet” boundary (the right boundary) of a square mesh. I know my solution i.e u_h is correct because it converges but it seems like I am doing something wrong with flux calculation since it is not giving me the results I expect.

The analytical flux I expect is grad(u_analy)*n and I choose u_analytical to be (x^2+y^2,0) in two dimensions. Since I am in right boundary of a square mesh my n is (1,0). So I expect flux analytical to be (2x , 0). The square has dimensions of [-1,1]^2 so the analytical result is 0. But mu numerical result is far from that… I am doing the boundary integration below…

[code] VH1= H1(mesh, order=2, dirichlet=“wall|inlet”)
gradu00 = GridFunction(VH1)
gradu01 = GridFunction(VH1)
gradu10 = GridFunction(VH1)
gradu11 = GridFunction(VH1)

gradu00.Set(grad(u_h.components[0])[0])
gradu01.Set(grad(u_h.components[0])[1])
gradu10.Set(grad(u_h.components[1])[0])
gradu11.Set(grad(u_h.components[1])[1])
     
et = GridFunction(VH1)
et.Set(CoefficientFunction(1.0), definedon=mesh.Boundaries("outlet"))
flux_x = Integrate(et*( gradu00*n[0] + gradu01*n[1] ), mesh, BND)
flux_y = Integrate(et*( gradu10*n[0] + gradu11*n[1] ), mesh, BND)
print ("flux_x:", flux_x)
print ("flux_y:", flux_y)

[/code]

Am I doing it wrong ?

2

Hello noname,

on the right boundary I expect x to be constant 1 and therefore the analytical flux to be 4, not 0.

What numerical value do you obtain?

Best,
Michael

3

Hi Michael,

I get almost 4 numerically… but I really do not follow why you expect 4 analytically …

4

With the same logic, on the circular boundary defined by level set function x^2+y^2-0.5^2=0 (using the same analytical solution) i calculate n to be <2x/sqrt(4x^2+4y^2) , 2y/sqrt(4x^2+4y^2)> and I calculate grad(u_analy)=<2x , 2y>
that makes analytical flux around the circle which is grad(u_analy).n to be
<4x^2/sqrt(4x^2+4y^2) + 4y^2/sqrt(4x^2+4y^2)> = 1. Am I also wrong here ?

5

  1. Let the domain [-1,1]^2 and u = (x^2+y^2,0). Then grad(u)= (2x, 2y; 0, 0). The right boundary is given by the set { (1,y) : -1<=y<=1} and the according normal vector is n=(1,0). Therefore, grad(u)*n = (2x, 0) on the right boundary. Further x=1 on the right boundary => grad(u)*n = (2, 0) on the right boundary. Integrating this flux over the right boundary (it has length 2) yields int(2, y=-1…1) = 4.
    This is what I also observe numerically.

  2. Let the domain given by the level-set function x^2+y^2-0.5^2=0. The normal vector is then given by n= (2x,2y) (simplifying your expression by using the level-set equation). Thus, grad(u)n = (4x^2 + 4y^2, 0) = (1,0). Integrating over the boundary of the circle yields int(1) = 2r*pi= pi.
    The numerical result is also near pi.

Attached you find the according python-code for both examples.

https://ngsolve.org/media/kunena/attachments/889/bnd_int_circ.py
https://ngsolve.org/media/kunena/attachments/889/bnd_int_rect.py

Best,
Michael

Attachment: bnd_int_circ.py

Attachment: bnd_int_rect.py

6

Thank you Michael, now it is very clear to me…