kill(all)$

fluid : stress strain  (= cylindrical coordinates ),version 1.0.
peter.vlasschaert@gmail.com,12/06/2019

load packages
=============

load(lrats)$

differentation  unit vectors: cylindrical
================================

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

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

cylindrical gradient  : gradf
=================

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

stress tensor :σ
==============
-------------------------------------------------------------------------
1e)  p  = pressure                             (matrix :3 by 3 )
2e)   τ  = deviatoric stress tensor (matrix : 3 by 3)
-------------------------------------------------------------------------
rem : p *(matrix : 3  by 3) = p*(unit matrix 3 by 3  )

symbolic : representation
======================
σσ = σ
ττ   =  τ

p2:σσ=-p*Id+ττ;

p2a:σσ=genmatrix(lambda([i,j], σ[i,j]), 3, 3);

p2a1:rhs(p2a);

element :σ ,tensor
----------------------------

p2a2:p2a1[1,1];

p2b:ττ=genmatrix(lambda([i,j], τ[i,j]), 3, 3);

unit matrix  ( = Id)
=========

p2c1:Id=matrix(
 [1,0,0], 
 [0,1,0], 
 [0,0,1]
);

pressure matrix:
==============

p2c2:matrix(
 [p,0,0], 
 [0,p,0], 
 [0,0,p]
);

p3:[p2a,p2b,p2c1];

a) symbolic presentation : σ ( = tensor)
------------------------------------------------------------

p3a:lhs(p2a) = lhs(p2b) - p*lhs(p2c1);

b) matrix symbolic presentation :σ ( = tensor)
-------------------------------------------------------------------

p3b:rhs(p2a) = rhs(p2b) - p*rhs(p2c1);

one element of matrix
---------------------------------

p3c:lhs(p3b)[1,1] =rhs(p3b)[1,1];

devatoric stress tensor ( "incompressible flow",cylindrical coordinates ) .
=============================================================
μ  = constant viscosity
u  = velocity

shear tensor :
============

τ = μ*(gradf(u) + transpose(gradf(u)))

independently of time : gradf
-------------------------------------------

p4: u = ur*er+uθ*eθ+uz*ez;

depends([ur,uθ,uz],[r,θ,z]);

gradf :  ?

p4a:gradc=gradf(rhs(p4));

p4b:expand(p4a);

use rule : differentation unit vectors ' p1'.
===================================

p4c:lhs(p4b) = subst(p1,rhs(p4b));

matrix : representation, 'gradc'
----------------------------------------------
how convert  p4c  => matrix

index notation : 'p4c'
==================

p4c1:[er = e[r],eθ = e[θ], ez = e[z],ur = u[r],uθ = u[θ],uz = u[z]];

p4c2:subst(p4c1,p4c);

p4c3:args(rhs(p4c2));

p4c4:length(p4c3);

How to analyse : 'p4c5'
-----------------------------------

p4c5:p4c3[1];

element analysis ? list element.
------------------------------------------------

p4c6a:part(p4c5,0);

p4c6b:part(p4c5,1);

p4c6c:part(p4c5,1,1);
p4c6d:part(p4c5,1,2);

p4c6e:part(p4c5,2);

use  : pickapart  ?

1e) example
------------------

pickapart(p4c5,3);

2e) example
------------------

pickapart(p4c3[4],3);

freeof(x1,x2,.....,expr)
-------------------------------

freeof(e[θ],p4c5);

op
--------

op(part(p4c5,1));

subst
----

subst(["."= "="],p4c3[4]);

coeff(%,e[θ]);

matrix :
======

p5:makelist(subst(["."= "="],p4c3[i]),i,1,length(p4c3));

if : construction ,check .

p5[11];

use : freeof, 'command'.

freeof(e[θ],lhs(p5[11]));

freeof(e[r],lhs(p5[11]));

find diagonal : gradc
--------------------------------

first: test ? only one value (with this code).
===================================
=>  e[θ] : one missing.

aa(m1,m2):= block([],
 for i: 1 thru length(m1)  
    do(if freeof(m2,lhs(m1[i])) = false and
           freeof(m2,rhs(m1[i])) = false then
        m3: coeff(lhs(m1[i]),m2)
    
        )
);

aa(p5,e[θ]);

m3;

second: test ? 
============
=>  solution:

bb(m1,m2):= block([],
    m4:[],
 for i: 1 thru length(m1)  
    do(  if freeof(m2,lhs(m1[i])) = false and
              freeof(m2,rhs(m1[i])) = false then
        m4: endcons(coeff(lhs(m1[i]),m2),m4)    
        ),
    return(m4)
);

bb(p5,e[θ]);

m4;

diagonal elements : gradc(u).
=========================

p6a:bb(p5,e[r]);
p6b:bb(p5,e[θ]);
p6c:bb(p5,e[z]);

How to find : off diagonal elements:gradc(u).
======================================

cc(m1,m2,m3):= block([],
    m5:[],
 for i: 1 thru length(m1)  
    do(  if freeof(m2,lhs(m1[i])) = false and
              freeof(m3,rhs(m1[i])) = false then
        m5: endcons(coeff(lhs(m1[i]),m2),m5)    
        ),
    return(m5)
);

p7a12:cc(p5,e[r],e[θ]);
p7a13:cc(p5,e[r],e[z]);
p7a21:cc(p5,e[θ],e[r]);
p7a23:cc(p5,e[θ],e[z]);
p7a31:cc(p5,e[z],e[r]);
p7a32:cc(p5,e[z],e[θ]);

matrix : gradc
============

p8:gradc=genmatrix(lambda([i,j], mm[i,j]), 3, 3);

mm[1,1]:sum(p6a[i],i,1,length(p6a));
mm[2,2]:sum(p6b[i],i,1,length(p6b));
mm[3,3]:sum(p6c[i],i,1,length(p6c));

input :diagonal elements
-------------------------------------

p8a:ev(p8,nouns);

input : off diagonal elements
---------------------------------------------

mm[1,2]:sum(p7a12[i],i,1,length(p7a12));
mm[1,3]:sum(p7a13[i],i,1,length(p7a13));
mm[2,1]:sum(p7a21[i],i,1,length(p7a21));
mm[2,3]:sum(p7a23[i],i,1,length(p7a23));
mm[3,1]:sum(p7a31[i],i,1,length(p7a31));
mm[3,2]:sum(p7a32[i],i,1,length(p7a32));

gradc(u)  ?
=========

p8b:ev(p8,nouns);

transpose(gradc(u)) ?
===================

p9:tgradc=transpose(rhs(p8b));

sum (gradc(u)+transpose(gradc(u))) ?
================================

p10:rhs(p8b)+rhs(p9);

rem : symmetric matrix :  m[i,j] =m[j,i], i≠j
rem : without μ before (every element).

see above :

p10a:p3b;

deviatoric tensor:τ
================

how to use : double for loop => 

for i:1 thru 3
do(
   for j:1 thru 3
do(τ[i,j]:μ*p10[i,j])
);

p10b:ev(p10a,nouns);

σ[1,1]:σ[r,r];
σ[1,2]:σ[r,θ];
σ[1,3]:σ[r,z];
σ[2,1]:σ[θ,r];
σ[2,2]:σ[θ,θ];
σ[2,3]:σ[θ,z];
σ[3,1]:σ[z,r];
σ[3,2]:σ[z,θ];
σ[3,3]:σ[z,z];

p10c:ev(lhs(p10b),nouns) = rhs(p10b);

empty list : to fill
==============

m6:[];

how to find element : ?
====================

p10c1:lhs(p10c)[1,1];

for i:1 thru 3
do(
   for j:1 thru 3
      do(
        m6:endcons(lhs(p10c)[i,j] = rhs(p10c)[i,j],m6)
          )
        
);

list of matrix , 'p10c'
=================

m6;