Equation of Motion for a Newtonian fluid with constant rho and mu

Thus del dot v = 0

Section 3.2 in BS&L

> restart;

> with(linalg):

> tau:=(mu,vx,vy,vz)->evalm(-mu*(matrix([[2*diff(vx,x),diff(vx,y)+diff(vy,x),diff(vz,x)+diff(vx,z)],[diff(vx,y)+diff(vy,x),2*diff(vy,y),diff(vy,z)+diff(vz,y)],[diff(vz,x)+diff(vx,z),diff(vy,z)+diff(vz,y),2*diff(vz,z)]]))); Eq. 3.2-11 to 16 with del dot v = 0.

> tau(mu,vx(x,y),vy(x,y),vz)[1,1]; The 1,1 element of tau.

> ddt:=t->vector([diff(t[1,1],x)+diff(t[2,1],y)+diff(t[3,1],z),diff(t[1,2],x)+diff(t[2,2],y)+diff(t[3,2],z),diff(t[1,3],x)+diff(t[2,3],y)+diff(t[3,3],z)]); del dot the tensor t

> ddt(tau(mu,vx(x,y,z,t),vy(x,y,z,t),vz(x,y,z,t)));

> delsqv:=v->evalm(map(diff,v,x$2)+map(diff,v,y$2)+map(diff,v,z$2)); del squared of a vector

> Ds:=(s,v)->diff(s,t)+v[1]*diff(s,x)+v[2]*diff(s,y)+v[3]*diff(s,z);

> DDv:=v->diverge(v,vector([x,y,z]));

> eqcon2a:=(rho,v)->Ds(rho,v)+rho*DDv(v); From p74

> cont:=eqcon2a(rho,vector([vx(x,y,z,t),vy(x,y,z,t),vz(x,y,z,t)]))/rho; For constant rho

> dif20:=evalm(mu*delsqv(vector([vx(x,y,z,t),vy(x,y,z,t),vz(x,y,z,t)]))+ddt(tau(mu,vx(x,y,z,t),vy(x,y,z,t),vz(x,y,z,t)))); The difference between mu*delsq(v) and del dot tau with tau given in eqs. 3.2-11 to 16.

> simplify(%);

> dif20:=simplify(evalm(simplify(%)/mu));

> simplify(evalm(%+grad(cont,vector([x,y,z])))); since cont is always 0 this is legitimate

and we find that del dot tau in 3.2-10 equals mu*del^2(v) for a Newtonian fluid with constant fluid properties.

> Dv:=(vin,v)->evalm(map(diff,vin,t)+v[1]*map(diff,vin,x)+v[2]*map(diff,vin,y)+v[3]*map(diff,vin,z)); divergence of a vector

> dels:=s->grad(s,vector([x,y,z])); del of a scalar

> eqmot:=evalm(rho*Dv(vector([vx(x,y,z,t),vy(x,y,z,t),vz(x,y,z,t)]),vector([vx(x,y,z,t),vy(x,y,z,t),vz(x,y,z,t)]))+dels(p(x,y,z,t))-mu*delsqv(vector([vx(x,y,z,t),vy(x,y,z,t),vz(x,y,z,t)]))-rho*vector([gx(x,y,z,t),gy(x,y,z,t),gz(x,y,z,t)])); full eq. of motion for Newtonian fluid with constant properties

> eqmot2d:=evalm(rho*Dv(vector([vx(x,y),vy(x,y),0]),vector([vx(x,y),vy(x,y),0]))+dels(p(x,y))-mu*delsqv(vector([vx(x,y),vy(x,y),0]))-rho*vector([0,gy,0]));

>