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