kill(all)$

fluid  : divergence of tensor cylindrical coordinates( = tensor),version 1.0.
peter.vlasschaert@gmail.com,07/06/2019

rem :version 1.0a
===============
-------------------------------------------------------------------
=>  krondelta(s1,s2)                   , rule
=>  kron_tensor(x,m1,m2,m3)  ,  x = er,eθ,ez
=>  fullkron_tensor(m1,m2,m3)
=>  comp_fullkron_tensor(xx,k,m1,m2,m3) 
-------------------------------------------------------------------
xx = ei ,k = i ; i=r,θ,z
m1,m2, m3 , 'see below'.
-------------------------------------------------------------------

rem : before 'part 3: tensor (= first        order tensor ) = vector '
                        'part 4 :tensor (= second order tensor )                 '

load packages

load(lrats)$

load(linearalgebra)$

kronecker delta : operator
---------------------------------------

1e) : kronecker delta operator: 'maxima'
=================================

p1:kron_delta(e1,e1);

p2a:kron_delta(e2,e1);
p2b:kron_delta(e1,e2);
p2c:kron_delta(x,y);

2e) : my kronecker delta 
=====================

krondelta(s1,s2):=block([],if s1=s2 then 1 else 0);

rem  :    s1 = s2   =>  1
               s1 ≠ s2   =>  0

p2a1:krondelta(er,eθ);
p2b1:krondelta(eθ,er);
p2c1:krondelta(ez,ez);
p2d1:krondelta(x,x);

divergence of a vector:operator
------------------------------------------------

depends([f],[r,θ,z]);

gradf : cylindrical coordinates
=========================

gradf(m) := diff(m,r)*er+1/r*diff(m,θ)*eθ+diff(m,z)*ez;

p3a:gradf(f);

p3a1:args(p3a);

p3a2:part(p3a1[1],1);
p3a3:part(p3a1[1],2);

tensor representation : cylindrical coordinates.

p3a4:[e1=er,e2=eθ,e3=ez];
p3a5:[a11=arr,a12=arθ,a13=arz,a21=aθr,a22=aθθ,a23=aθz,a31=azr,a32=azθ,a33=azz];
p3a6:flatten([p3a4,p3a5]);

simp:false;

p3b1: a11.e1.e1+a12.e1.e2+a13.e1.e3;
p3b2: a21.e2.e1+a22.e2.e2+a23.e2.e3;
p3b3: a31.e3.e1+a32.e3.e2+a33.e3.e3;
p3b:subst(p3a6,flatten(p3b1+p3b2+p3b3));

result :
======
.  :  non commutative operator , see , 'p3b'
see: part3.  (why)

non commutative operator :
=======================

p3d1:args(p3b);

same : representation of  'p3d1' of a tensor.
------------------------------------------------------------------

p3d1a:subst(p3a4,[e1,e1,e1,e2,e2,e2,e3,e3,e3]);

p3d1b:subst(p3a4,[e1,e2,e3,e1,e2,e3,e1,e2,e3]);

p31c:subst(p3a5,[a11,a12,a13,a21,a22,a23,a31,a32,a33]);

p31c1:[a[r,r],a[r,θ],a[r,z],a[θ,r],a[θ,θ],a[θ,z],a[z,r],a[z,θ],a[z,z]];

solve : problem non commutative operator
------------------------------------------------------------------

1e) parts : nabla operator.
======================

p4a1:part(p3a1[1],1,1);
p4a2:part(p3a1[2],1);
p4a3:part(p3a1[3],1);

simp:true;

2e) use krondelta
===============

by : ez
----------

p5a1:makelist(krondelta(p4a1,p3d1a[i]),i,1,length(p3d1a));

p5a11:makelist(p5a1[i]*p31c1[i],i,1,length(p3d1a));

p5a11a:makelist(p5a1[i]*p3d1b[i],i,1,length(p3d1a));

p5a11b:makelist(p5a11[i]*p5a11a[i],i,1,length(p3d1a));

by : eθ
----------

p5a2:makelist(krondelta(p4a2,p3d1a[i]),i,1,length(p3d1a));

p5a22:makelist(p5a2[i]*p31c1[i],i,1,length(p3d1a));

p5a22a:makelist(p5a2[i]*p3d1b[i],i,1,length(p3d1a));

p5a22b:makelist(p5a22[i]*p5a22a[i],i,1,length(p3d1a));

by : er    
----------

p5a3:makelist(krondelta(p4a3,p3d1a[i]),i,1,length(p3d1a));

p5a33:makelist(p5a3[i]*p31c1[i],i,1,length(p3d1a));

p5a33a:makelist(p5a3[i]*p3d1b[i],i,1,length(p3d1a));

p5a33b:makelist(p5a33[i]*p5a33a[i],i,1,length(p3d1a));

fa:makelist(krondelta(er,p3d1a[i]),i,1,length(p3d1a));
fb:makelist(fa[i]*p31c1[i],i,1,length(p3d1a));
fc:makelist(fa[i]*p3d1b[i],i,1,length(p3d1a));
fd:makelist(fb[i]*fc[i],i,1,length(p3d1a));

kron_tensor(x,m1,m2,m3):=block([],    
faa:makelist(krondelta(x,m1[i]),i,1,length(m1)),
fbb:makelist(faa[i]*m2[i],i,1,length(m1)),
fcc:makelist(faa[i]*m3[i],i,1,length(m1)),
fdd:makelist(fbb[i]*fcc[i],i,1,length(m1))
);

simp;

kron_tensor(er,p3d1a,p31c1,p3d1b);

combine
========

p6:map("+",p5a11b,p5a22b,p5a33b);

fullkron_tensor(m1,m2,m3):= block([],   
    fe1:kron_tensor(er,m1,m2,m3),
    fe2:kron_tensor(eθ,m1,m2,m3),
    fe3:kron_tensor(ez,m1,m2,m3),
    map("+",fe1,fe2,fe3) );
    

same : result 'p6'

-----------------------------------------------------------------------------------------------------------------------------------

p7:fullkron_tensor(p3d1a,p31c1,p3d1b);

------------------------------------------------------------------------------------------------------------------------------------

p7a:map("=",p31c1,p31c);

p7a1:map("=",p31c,p31c1);

p7b:subst(p7a,p7);

We want to find general equation  : ei  component , i=r,θ,z
=================================================

rule : unit vectors : ei , i=r,θ,z
===============

p8:['diff(er,θ) = eθ,'diff(eθ,θ) = -er];

p8a:depends([er,eθ],θ);

t:3;

differentation with respect to 'r'
==============
ari  = > i=r,θ,z

p9a:sum(p7b[i],i,1,t);

depends([arr,arθ,arz],r);

p9b:diff(p9a,r);

use :rule.
---------------

p9c:subst(p8,p9b);

differentation with respect to 'θ'
==============
aθi  = > i=r,θ,z
rem :  θ no distance => multiplied by '1/r'
see : gradf =>  'above'

p10a:sum(p7b[i],i,t+1,2*t);

depends([aθr,aθθ,aθz],θ);

p10b:1/r*diff(p10a,θ),expand;

use:rule
--------------

p10c:subst(p8,p10b);

differentation with respect to 'z'
==============
azi  = > i=r,θ,z

p11a:sum(p7b[i],i,2*t+1,3*t);

depends([azr,azθ,azz],z);

p11b:diff(p11a,z),expand;

use:rule
--------------

p11c:subst(p8,p11b);

final result : divergence of tensor cylindrical coordinates
--------------------------------------------------------------------------------------
a= tensor , second order.

sum :differentation with respect to ,'r,θ,z'
===================================

p12:p9c+p10c+p11c;

p12a:collectterms(p12,er,eθ,ez);

er :  component
==============
corrected term : arr/r ? Why

p12a1:(arr/r+ratcoeff(p12a,er)),exapand;

p12a11:expand(%);

p12a11a:subst(p7a1,p12a11);

eθ :  component
==============
corrected term : arθ/r ? Why

p12a2:arθ/r+ratcoeff(p12a,eθ),expand;

p12a11a:subst(p7a1,p12a2);

ez :  component
==============
corrected term : arz/r  ? Why

p12a3:arz/r+ratcoeff(p12a,ez),expand;

p12a33a:subst(p7a1,p12a3);

tensor (=second order)  => tensor ( = first order)  = vector .


result: (website)
==============

---------------------------------------------------------------------------------------------------------------
part 3:div(vector,  'first        order')  = 3 components
part 4:div(tensor, 'second  order')  = 9 components  = (3*3)     components
----------------------------------------------------------------------------------------------------------------
general :
            div(tensor, 'thrid        order') = 27 components = (3*3*3) components
----------------------------------------------------------------------------------------------------------------
example :second order :stress tensor used , navier - stokes equation.
-----------------------------------------------------------------------------------------------------------------