Warping function for 3D beams

Dear All,

I am modeling 3D beams with arbitrary cross sections and for that purpose I need to calculate the warping function numerically. The warping function ψ(y,z) is defined as a solution of the following Laplace equation with Neumann boundary conditions:

∂2ψ/∂y2 + ∂2ψ/∂z2 = 0 in Ω,
∂ψ/∂n = z*ny - y*nz on Γ,

where Ω is the cross sectional area, and Γ is its contour. This system can be solved by different numerical methods, for example, finite element method, boundary element method, finite differences...
Numerically I obtain that

∫ψdΩ = 0,

for non-symmetrical cross sections and it is true for the warping function obtained analytically for beams with rectangular cross sections, elliptical or triangular.
I wonder if someone can help me to prove that ∫ψdΩ = 0 for any cross section? Furthermore,

∫y*ψdΩ = 0 and
∫z*ψdΩ = 0.

And it should be possible to prove also these integrals.
For the center of the coordinate system is chosen the twist center (or shear center) of the cross section and the warping function is zero at that point.

Thank you in advance!

Regards,
Stan