kill(all)$

maxima:fvm:appendix:combine(a[E]/D[e],a[W]/D[w]),powerlaw, part 3,version 1.0
peter.vlasschaert@gmail.com,16/12/2022.

load packages
-------------------

load(lrats)$

see : website
-----------------
ref 1:
linearized exponential scheme (=hybrid scheme,a[E]/D[e]) part1,version 1.0
ref 2:
linearized exponential scheme (=hybrid scheme ,a[W]/D[w]) part 2,version 1.0
see : last part of part 2.

info :
point of intersection : y =1 ,both curve ( mirror curves,around y=1)
------------------------------------------------------------------------
----------------------------------------------------------------------------------------
blue curve  : part 1 → a[E]/D[e]    = max(-Pe,1-Pe/2+1,0) 
blue curve  : part 2 → a[W]/D[w] = max(Pe,1+Pe/2,0 )
-----------------------------------------------------------------------------------------

see : ref:1 & ref:2

p1a:a[p]*Φ(p) = a[W]*Φ(W) +a[E]*Φ(E);
p1b:a[p]=F[w]/(%e^Pe(Δx(w))-1)+F[e]·%e^Pe(Δx(e))/(%e^Pe(Δx(e))-1);
p1c:a[W]=F[w]·%e^Pe(Δx(w))/(%e^Pe(Δx(w))-1);
p1d:a[E]=F[e]/(%e^Pe(Δx(e))-1);

Definition : Pe(.), peclet cell number
------------------------------------------------

p1e:Pe(m1) = F[m]/D[m];

follows :

→ interface : CS =[e], e =east.

p1ee:reverse(subst([m1=Δx(e),m=e],p1e));

p1ee1 :reverse(solve(p1ee,F[e]))[1];

p1ee1a:reverse(p1ee1);

use : ref 1:p2a4

p2e:a[E]/D[e]=Pe(Δx(e))/(%e^Pe(Δx(e))-1);

→ mutiplication "D[e]"

p2e1:ratsimp(p2e*D[e]);

p2e2:lratsubst([reverse(p1ee1)],p2e1);

→ interface : CS =[w], w =west

p1ew:reverse(subst([m1=Δx(w),m=w],p1e));

→ multiplication "D[w]"

p1ew1:reverse(solve(p1ew,F[w]))[1];

use :ref 2:p2a4

p2w:a[W]/D[w]=(%e^Pe(Δx(w))·Pe(Δx(w)))/(%e^Pe(Δx(w))-1);

p2w1:ratsimp(p2w*D[w]);

p2w2:lratsubst([reverse(p1ew1)],p2w1);

with notation above:

interface : CS =[e], face of CV
----------------------------------------
blue curve  : part 1 → a[E]/D[e]    = max(-Pe,1-Pe/2+1,0) 

p3e:a[E]/D[e]    = max(-Pe,1-Pe/2,0);

p3e1:subst([Pe=Pe(Δx(e))],p3e);

p3e2:lhs(p3e1) = args(rhs(p3e1) );

p3e3:expand(p3e2*denom(lhs(p3e2)));

use :p1ee1

p3e3a:expand(lratsubst([reverse(p1ee1)],p3e3));

p3e3b:lhs(p3e3a) =apply(max,rhs(p3e3a));

interface : CS =[w], face of CV
----------------------------------------
blue curve  : part 2 → a[W]/D[w] = max(Pe,1+Pe/2,0 )

p3w:a[W]/D[w] = max(Pe,1+Pe/2,0 );

p3w1:subst([Pe=Pe(Δx(w))],p3w);

p3w2:lhs(p3w1) = args(rhs(p3w1) );

p3w3:expand(p3w2*denom(lhs(p3w2)));

use :p1ew1

p3w3a:expand(lratsubst([reverse(p1ew1)],p3w3));

p3w3b:lhs(p3w3a) =apply(max,rhs(p3w3a));

rem :a[p] = a[W]+a[E]+(F[e]-F[w])
conservation law : (F[e]-F[w]) =0

p4:a[p] =a[E]+a[W];

p4a:p3w3b;

p4b:p3e3b;

rem : powerlaw : more linearized exponential scheme
----------------------
https://www.wikiwand.com/en/Power_law_scheme

ref : Patankar (1980)

a[E] = -F[e]                                    → Pe≤ -10
a[E] = D[e]*((1+0.1*Pe)^5-Pe)       → -10≤Pe≤0
a[E] = D[e]*((1-0.1*Pe)^5)             →  0≤Pe≤10
a[E] =0                                           → Pe ≥ 10

a[E] = D[e]*(max(0,(1-0.1*Pe)^5)+max(0,-F[e]))

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