%
%  symdrag.m     Steve Cox, Feb 10, 2009
%
%  Compute the drag force on a sphere of radius "a" at low Reynold's
%  number using the Stoke's approximation developed in our
%  drag notes. In particular, we compute the integrand
%
%   Pin = p*n - sig*n   where n = -x/|x| is the normal into the sphere
%                       p is the pressure and sig is the viscous
%                       part of the stress, sig = eta*(grad(v)+grad(v)')
%
%  we find that Pin = (3/2)(U*eta/a)*[1 0 0] is constant on the sphere
%  and so its integral over the sphere is simply (4/3)*pi*a^2*Pin
%

syms x1 x2 x3 U a eta

x = [x1; x2; x3];
 
X = x*transpose(x);
 
I=eye(3); 

u = [U; 0; 0];

r = sqrt(x1^2+x2^2+x3^2);
 
v = u - (3/4)*(a/r)*(I+X/r^2)*u - (1/4)*(a^3/r^3)*(I-3*X/r^2)*u;
 
vx1 = diff(v,x1);
vx2 = diff(v,x2);
vx3 = diff(v,x3);             

gradv = [vx1 vx2 vx3];
 
sig = eta*(gradv+transpose(gradv)); 

sigtn = simple(sig*(-x/r));

sigtns = simple(subs(sigtn,a,r))

p = -(3/2)*eta*(1/r^2)*transpose(u)*x;
pn = p*(-x/r)

Pin = simple(pn - sigtns)