Finding an approximate solution to this integral: $alpha(phi,r,p,d)=int_0^infty w(z,r,p,d)Q(z,r,phi)dz$

**** See update as per comments below ****

Lots of functions to fit numerical data (I use FindFormula below). Below I start with setting up the integral function with integration limits of $1leq zleq 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}]

Update:

As per comments below, FindFormula only works with one variable. So the following uses FindFit and applies data of myIntFun[r,1,2,$phi$] to: $$ a+bphi+cphi^2+dr+er^2+frphi $$ in the range of $0.1leq rleq 0.9$ and $0.1leq phileq 0.9$ and compares the fitted formula to a ListPointPlot3D of the data points:

phiTable = Table[
   {phi, r, myIntFun[r, 1, 2, phi]},
   {phi, 0.1, 0.9, 0.05}, {r, 0.1, 0.9, 0.05}];


myParms = 
  FindFit[Flatten[phiTable, 1], 
   a + b x + c x^2 + d y + e y^2 + f x y, {a, b, c, d, e, f}, {x, y}];

myFit[x_, y_] = (a + b x + c x^2 + d y + e y^2 + f x y) /. myParms;

my2DPlot = 
  Plot3D[myFit[x, y], {x, 0, 1}, {y, 0, 1}, 
   PlotStyle -> {Opacity[0.2], Blue}];

lp = ListPointPlot3D[phiTable, BoxRatios -> {1, 1, 1}];
Show[{lp, my2DPlot}]