I am trying to solve Laplace equation over a time-dependent region using cylindrical coordinates. My system is independent of the angle coordinate, and so the equation looks like:

The region represents a container filled with a liquid plus a sphere that falls into it. The bottom and the sides of the container remain fixed. I have used the Finite Element Method to define my region (I have only represented half of it since it is independent of the angle):

Needs["NDSolve`FEM`"]

(*constants*)
ballRadiusSquared = 0.5;
rmax = 10.;
zmin = 0.;
zinit = 4.;
pa = 1.;
g = 9.8;
rho = 1.;
tmax = 0.2;

(*time dependent variable*)
zBallCentre[t_] := -0.5*g*t^2 + zinit + Sqrt[ballRadiusSquared];

(*define region coordinates*)
n = 1000;
(*circle coordinates*)
circTbl = Table[{Sqrt[ballRadiusSquared]*Cos[t], Sqrt[ballRadiusSquared]*Sin[t]}, {t, -Pi/2, Pi/2 - 2 Pi/(2 n), 2 Pi/(2 n)}];
Do[circTbl = ReplacePart[circTbl, {{i, 2}} -> circTbl[[i, 2]] + zBallCentre[0.7]], {i, Length[circTbl]}];
Do[If[circTbl[[i, 2]] > zinit, circTbl = Delete[circTbl, i], circTbl = circTbl], {i, Length[circTbl], 1, -1}];

(*surface coorinates*)
 surfaceTbl = Table[{x, zinit}, {x, 0., rmax, 1/n}];
If[Length[circTbl] == n, surfaceTbl = surfaceTbl, 
    Do[If[circTbl[[Length[circTbl], 1]] > surfaceTbl[[1, 1]],surfaceTbl = Delete[surfaceTbl, 1],surfaceTbl = surfaceTbl],
    {i, 1, n}]];

(*put everything together*)
coordTbl = Join[circTbl, surfaceTbl];
coordTbl = Prepend[coordTbl, {0., zmin}];
coordTbl = Append[coordTbl, {rmax, zmin}];
lineTbl = Partition[Range[Length[coordTbl]], 2, 1];
lineTbl = Prepend[lineTbl, {Length[coordTbl], 1}];

(*mesh the region where de PDE will be solved*)
bmesh = ToBoundaryMesh["Coordinates" -> coordTbl, "BoundaryElements" -> {LineElement[lineTbl]}];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 0.1];
mesh["Wireframe"]

With the following output for t=0.7:

I have used many points for the upper part because it is going to change a lot over time.
Now, I want to impose a Neumann boundary condition on the sphere. The Neumann boundary condition will be given by the component of the velocity of the sphere that is perpendicular to its surface at a given point. So, I have to find a way to impose a different Neumann boundary condition for each of the lines that join the points around the sphere (the points are defined in circTbl).
Do you know how to do it?

Note:
The other boundary conditions that I am using are:

DirCond = DirichletCondition[\[Phi][r, z] == 0., r == rmax];
NeumCond = NeumannValue[0., r == 0] + NeumannValue[0., r == rmax] + NeumannValue[0., z == zmin];

And I am solving the Laplace equation in the following way:

solution = NDSolveValue[{D[r*D[\[Phi][r, z], r], r]/r + D[\[Phi][r, z], {z, 2}] == NeumCond, DirCond}, \[Phi], {r, z} \[Element] mesh, "ExtrapolationHandler" -> {Automatic, "WarningMessage" -> False}]