Lots of functions to fit numerical data (I use FindFormula below). Below I start with setting up the integral function with integration limits of $1\leq z\leq 10$ for starters and vary just $\phi$ generating a table of {$\phi$,integralFunction[r,p,d,$\phi$]}. Then use FindFormula to fit this data. 

Then I would try to vary two variables like $\phi$ and $r$, generate a table of {$\phi$,r,integralFunction[r,p,d,$\phi$]} and then try to fit this 3D data. Then try to generate a table of another variable and so on and also increase the integration limits. Kinda brute force but it's something you may wish to work with.
 
Here is a start just varying $\phi$ from 0.01 to 0.99. The red fitted function fits the points nicely. 

 (*
 define functions
 *)
 wFun[z_, r_, p_, d_] := (
 r (d + (r + z) Sech[(z + p)/d]^2 - d Tanh[(z + p)/d]) (z + 
 r Tanh[(z + p)/d]))/(d (r + z)^3 Sech[p/d]^2);
 qFun[z_, r_, \[Phi]_] := -2 Exp[(
 z \[Phi] (6 r^2 (\[Phi] - 1)^2 - 3 r z (\[Phi]^2 + \[Phi] - 2) + 
 2 z^2 (\[Phi]^2 + \[Phi] + 1)))/(2 r^2 (\[Phi] - 1)^3)];
 (*
 define integral function
 *)
 myIntFun[r_?NumericQ, p_?NumericQ, d_?NumericQ, \[Phi]_?NumericQ] := 
 NIntegrate[wFun[z, r, p, d] qFun[z, r, \[Phi]], {z, 0, 10}];
 
 (*
 for now, create table varying just phi from 0.01 to 0.99 and fit a \
 formula to this 1D data.
 *)
 phiTable = Table[
 {phi, myIntFun[1, 1, 2, phi]},
 {phi, 0.01, 0.99, 0.01}];
 
 (*
 find a formula for data
 *)
 
 theF[x_] = FindFormula[phiTable, x]
 (*
 superimpose ListPlot of points with fitted function
 *)
 lp = ListPlot[phiTable, PlotStyle -> {Black, PointSize[0.01]}]
 p1 = Plot[theF[x], {x, 0.01, 0.99}, PlotStyle -> Red]
 Show[{lp, p1}]

[![enter image description here][1]][1]


 [1]: https://i.sstatic.net/7GhvB.jpg