Hello everyone,

recently I have been struggling to find a satisfying result to an integral that I want to solve. I will first describe the background so you get an idea of what I am trying to achieve and then go on to explain my attempts, thoughts and troubles in solving it.

I want to find a symbolic expression for the volume of a liquid column in a tube. The column is cylindrical in shape, except for the top surface which I want to model as a biconic function. This is my starting point:

z = (cx x^2 + cy y^2)/(1 + Sqrt[1 - (1 + kx) cx^2 x^2 - (1 + ky) cy^2 y^2])

To model the height of the column I add the variable h. For now I am only interested in "circular" curvatures so I set the conic constants in x and y to zero. To make it more easily understandable, I replace the curvatures with radii. My idea to find the volume of the cylindrical column is to calculate a circular integral so I replace x and y with the corresponding polar coordinates.

zint = z + h /. {kx -> 0, ky -> 0, x -> r Cos[\[Phi]], y -> r Sin[\[Phi]], cx -> 1/Rx, cy -> 1/Ry}

The resulting function models the surface shape exactly how I want. I verified this visually with plots and manually by plugging in the numbers to model a half-sphere and then calculating the integral. The result gives the exact numerical value that you would get if you took the geometrical formula for the volume of a sphere.

So now the interesting part begins. I set up a multiple intgral, first integrating over the radius (from zero to the half-diameter of the tube) and then over the angle (zero to 2 pi).

vol = Integrate[(zint*r), {\[Phi], 0, 2 \[Pi]}, {r, 0, dia/2}]

This integral however, takes ages to compute and did not give me a result even after letting it run for two days. So then I researched how to speed up the calculation of symbolic integrals and found this thread suggesting that I calculate the indefinite integral first and find the values in a separate step. This approach did give me an expression for the indefinite integral and I also received expressions when calculating the limits. When plugging in numbers the volume is incorrect, sadly.

indef = Integrate[(zint*r), \[Phi], r]
tmp = Limit[indef, r -> dia/2]
vol = Limit[tmp, \[Phi] -> 2 \[Pi]]

My understanding is that by using this approach I indirectly skipped the continuity check that Integrate uses and which might be the reason that the definite integral takes so long to compute. So did I make some kind of faulty simplification here? Possibly because of the division by zero when either Rx or Ry are equal to zero?

My question now is, is there a way to find a symbolic expression that is valid for all values of Rx and Ry (negative and positive but real)? If no, how can I split the Integral? I am also open to other suggestions how to obtain the desired function (relation between curvatures, height and volume).

I can also upload the notebook if it helps in any way. Thank you for reading