(* 3D vector field *)
F={2x+3y+4z,x^2+y^2+z^2,x^3+y^3+z^3}
(* sphere around (x,y,z) with radius r *)
q={x+r Sin[ph]Cos[th],y+r Sin[ph]Sin[th],z+r Cos[ph]};
volume=4Pi r^3/3;

(* flux density *) 
qth=D[q,th];
qph=D[q,ph];
n=Simplify[Cross[qph,qth]]
Foq=Simplify[F/.{x->q[[1]],y->q[[2]],z->q[[3]]}];
flux=Integrate[Integrate[Dot[Foq,n],{th,0,2 Pi}],{ph,0,Pi}]
Limit[flux/volume,r->0]

Div[F,{x,y,z}]