Making Dimensionless Models

Anton Antonov
IMM at DTU, UNI-C, NERI, Denmark
Anton.Antonov@uni-c.dk

When are solved problems in physics and chemical engineering, often it is more convenient to transform their model equations in dimensionless form. The package presented here contains some functions that help that process. We will show the package usage by examples.

First we should load the package

<< "~/MathFiles/dimensionless/DimLessPack.m"

The package can be also loaded from the section "The Mathematica functions"at the end of the notebook.

Notation definitions

<< Utilities`Notation`

Notation[x_ ⟨y___⟩ ⟺ Measure[x_, y___]]

Notation["\[DifferentialD]" x___ ⟺ DifferentialD[x___]]

The notation definitions above define that, for example,  C ⟨M⟩ will be interpreted as M to be the measure of C. We use angle brackets,  so we should type C < M > .

An example equation

\(\(Next\ let\ us\ take\ an\ equation\ from\ the\ book\ \ ' 
          Applied\ Mathematics\ and\ Modeling\ for\ \ Chemical\ Engineers'\ \ by\ Rice\  &\)\ Do, \ Chapter\ 1, \ page\ 30, \ 
    problem\ 1.4\_1\)

eqn = 
        m ⟨Kilogram⟩*(Subscript[C, p]) ⟨((4.1868`  (Meter^2))/((Second^2) )) Temperature ⟩*(\[DifferentialD](T/\[DifferentialD])t) ⟨(Temperature/Second)⟩ == 
          h ⟨((Kilogram  Temperature )/((Second^3) ))⟩*
            A ⟨(Meter^2)⟩*((Subscript[T, f]) ⟨Temperature⟩ - 
                T ⟨Temperature⟩);

It is instructional to see the  Mathematica internal presentation of this equation.

Here is the FullForm of the variable eqn defined above

eqn // FullForm

We can multiply arbitray physical quantities, but we can just add quantities that have the same measure. Therefore we can extract from the equation the relation imposed on its terms. This realtion is extracted by the function MaxRelation.

This transforms an equation to its relational form

MakeRelation[eqn]

We might also want to extract the units used in the equation. This is done by the function Units.

This gives the units used in the equation

Units[eqn]

Applying the Pi-theorem

The Buckingham's Pi theorem states that physical laws are independent of the form of the units. Therefore, acceptable laws of physics are homogeneous in all dimensions. (See Eric Weisstein's http://www.treasure-troves.com/physics/BuckinghamsPiTheorem.html.)

After the model equations are derived, we can transform them using units inherent for the problem modeled. This "inherentness" leads to the dimensionless form of the equations.

One way to look at a transformation of an equation to its dimensionless form is to think that  it is changed the basis in which its relational form is expressed. This algebraic point of view is taken in G. Birkhoff's "Hydrodynamics: A Study in Logic, Fact and Similitude". Princeton University Press, Revised edition, 1960. I have used the proof of the Buckingham's Pi theorem to make the function DimlessRel demonstrated bellow. The theorem states existence of a dimensionless form of an equation, so the Birkhoff's proof, although constructive at large extend, should be made fully constructive in order to be implemented in a computer algebra system. (So, actually, the code of DimlessRel is another proof of the theorem.)

DimlessRel[MakeRelation[eqn], {A, m, T}]

After we have found the new units we can substitute them in the relation form of the equation.

(MakeRelation[eqn] /. #1 // Simplify) & /@ %

Following the basis change point of view, we can make the function DimlessRel to produce several new bases when are passed to it non-independent system of vectors (i.e. units).

DimlessRel[MakeRelation[eqn], {(Subscript[C, p]), T, A, m}]

Finding the dimensionless form

The steps described in the previous sections are packed in the function Dimless. Dimless gives the dimensionless form of an equation using a list of units with which to express the units kept the global variable FUNDAMENDAL. So far the package works with CI units i.e. the equation we want to make dimensionless shoud be in CI units.

Global units variable

FUNDAMENTAL

The two equations are derived replacing the CI units, with a set of quantities that contain subset of {T,m}. Check the result of DimlessRel[MakeRelation[eqn],{T,m}]

Dimless[eqn, {T, m}]

The two equations are derived replacing the CI units, with a set of quantities that contain {(Subscript[C, p])}.

Dimless[eqn, {(Subscript[C, p])}]

All dimensionless forms chosen by the function ChooseMinor

Dimless[eqn, {}]

Choosing the simplest dimensionless equation

We can choose the simplest equations just looking and them as trees and choosing the shortest. The tree for the expression (a+b)*c is

(a + b)*c  // TreeForm

and the number of nodes in that tree is 5  (Times, Plus,a,b,c)

(a + b)*c  // LeafCount

This is the same as the last evaluation above

dimLessEqs = Dimless[eqn, {}];

Select[dimLessEqs, 
      Min[LeafCount /@ dimLessEqs] == LeafCount[#1] &]

The Mathematica  functions

Import["~/MathFiles/dimensionless/DimLessPack.ps","EPS"]//Show;

Dependency graph

<< Miscellaneous`Units`
<< Miscellaneous`PhysicalConstants`

\(\(ToFundamental\  = \ 
\ \ Dispatch[{Radian\  -> \ 
              1, \ 
\ \ \ \ Steradian\  -> \ 
              Radian^2, \ 
\ \ \ \ Newton\  -> \ \((Kilogram*Meter)\)/
                Second^2, \ 
\ \ \ \ Pascal\  -> \ 
              Newton/Meter^2, \ 
\ \ \ \ Joule\  -> \ 
              Meter*Newton, \ 
\ \ \ \ Watt\  -> \ 
              Joule/Second, \ 
\ \ \ \ Coulomb\  -> \ 
              Ampere*Second, \ 
\ \ \ \ Volt\  -> \ 
              Watt/Ampere, \ 
\ \ \ \ Ohm\  -> \ 
              Volt/Ampere, \ 
\ \ \ \ Siemens\  -> \ 
              Ampere/Volt, \ 
\ \ \ \ Farad\  -> \ 
              Coulomb/Volt, \ 
\ \ \ \ Weber\  -> \ 
              Second*Volt, \ 
\ \ \ \ Henry\  -> \ 
              Ohm*Second, \ 
\ \ \ \ Tesla\  -> \ 
              Weber/Meter^2, \ 
\ \ \ \ Lumen\  -> \ 
              Candela*Steradian, \ 
\ \ \ \ Lux\  -> \ 
              Lumen/Meter^2, \ 
\ \ \ \ Hertz\  -> \ 
              Second^\((\(-1\))\), \ 
\ \ \ \ Becquerel\  -> \ 
              Second^\((\(-1\))\), \ 
\ \ \ \ Gray\  -> \ 
              Joule/Kilogram}];\)\)

Clear[AtomicExponent]
AtomicExponent[expr_, 
        x_] :=  
If[#1 == {}, {0}, #1] &@
            If[#1 == {}, Cases[expr, s : x -> 1, 1], #1] &[
        Cases[expr, s : Power[x, n_] -> n, Infinity]]

Clear[AtomicVariables]
AtomicVariables[x_] :=  
Union[
        Cases[x, 
          s_ /;  (! NumericQ[s]  &&  AtomQ[s] || 
                MatchQ[s, 
                  Subscript[
                    s1_ /; (! NumericQ[s1]  &&  AtomQ[s1]), _]]), 
          Infinity]]

Clear[CleanSubscripts]
CleanSubscripts[
        x_] := 
Module[{pos}, 
Take[
                x, #1] & /@ (Position[#1, {}] & @(Flatten /@ 
                        Outer[Cases[#2, #1, Infinity] &, #1, #1]) &[
                x]) // 
          Flatten
]

FUNDAMENTAL = 
      AtomicVariables[(#1[[2]] //. 
                ToFundamental) & /@ Definition[
                ToFundamental][[1]][[1]]]

Clear[FindMinors, ChooseMinor]
FindMinors[mat_, highestRank_, 
        lowestRank_:  1] := 
Module[{d, mmat, amat, rules, 
          a}, 
{rows, cols} = Dimensions[mat]; 
        mmat = If[rows < cols, 
            Join[mat, Table[0, {i, rows + 1, cols}, {j, 1, cols}]]]; 
        mmat = If[rows > cols, 
            Transpose @ 
              Join[Transpose[mat], 
                Table[0, {i, cols + 1, rows}, {j, 1, 
                    rows}]]]; 
amat = 
          Array[a, Dimensions[mmat]]; 
rules = 
          Thread[Flatten[
                Table[Minors[amat, i, Identity], {i, highestRank, 
                    lowestRank, -1}], 2] -> 
Flatten[
                Table[Minors[mmat, i, Identity], {i, highestRank, 
                    lowestRank, -1}], 2]]; 
Select[
            rules, Det[#1[[2]]] ≠ 0 &] /. 
          a[x_, y_] -> {x, 
              y}
]

ChooseMinor[minorRules_] := 
First[
        minorRules]

Clear[MakeDimlessRules]
MakeDimlessRules[eqn_, minor_, basisIndex_, 
        units_] := 
Thread[units -> (Exp /@ 
                  LinearSolve[
                    minor[[2]], -Log[eqn[[#1]] //. 
                                ToFundamental /. 
                              Thread[
                                units -> 
                                  Table[1, {Length[units]}]]] & /@ 
                      basisIndex]) // PowerExpand // 
          Simplify]

Clear[DimlessBasis]
DimlessBasis[eqn_, ChooseMinor_, 
        DoWithBasis_] := 
Module[{B, units, minor, basis, 
          basisIndex}, 
units = 
          Intersection[FUNDAMENTAL, 
            AtomicVariables[
              eqn //. ToFundamental]]; 
B = PowerExpand[#1 //. ToFundamental] & /@ List @@ 
              eqn; 
B = 
          Flatten  /@  
            Outer[AtomicExponent[#1, #2] &, B, 
              units]; 
minor = 
          ChooseMinor[
            FindMinors[B, 
              Length[units]]]; 
basisIndex = #1[[1, 
                  1]] & /@ 
            minor[[1]]; 
basis = B[[#1]] & /@ 
            basisIndex; 
DoWithBasis[eqn, minor, 
          basisIndex, 
          Extract[units, ({#1[[2]]} & /@ 
                minor[[1, 1]])]]
]

Clear[BMatrix]
BMatrix[
        eqn_] := 
Module[{B, units, minor, basis, 
          basisIndex}, 
units = 
          Intersection[FUNDAMENTAL, 
            AtomicVariables[
              eqn //. ToFundamental]]; 
B = PowerExpand[#1 //. ToFundamental] & /@ List @@ 
              eqn; 
B = 
          Flatten  /@  
            Outer[AtomicExponent[#1, #2] &, B, 
              units]
]

Clear[Units]
Units[eqn_] := 
        Intersection[FUNDAMENTAL, 
          AtomicVariables[eqn //. ToFundamental]];

Clear[Inclusion]
Inclusion[l1_, 
        l2_] := #1 == Sort[l1] || #1 == Sort[l2] &[
        Sort[Intersection[l1, l2]]]

Clear[DimlessRel]
DimlessRel[eqn_, constants_: {}] /; 
        And @@ (MemberQ[eqn, #1, Infinity] & /@ constants) || 
          constants === {} := 
Module[{units, minors, pos, 
          B}, 
units = 
          Units[eqn]; 
pos = Position[eqn, #1][[1, 
                  1]] & /@ constants; 
B = 
          BMatrix[eqn]; 
minors = 
          Select[FindMinors[B, Length[units], Length[units]], 
            Inclusion[#1[[1, 1]] & /@ #1[[1]], 
                pos] &]; 
DimlessBasis[eqn, 
              Function[x, Evaluate[#1]], MakeDimlessRules] & /@ 
          minors
]

(*<< Utilities`Notation`
Notation[x_ ⟨y___⟩ ⟺ Measure[x_, y___]]
Notation["\[DifferentialD]" x___ ⟺  DifferentialD[x___]]*)

Clear[MakeRelation]
MakeRelation[eqn_] := 
f @@ 
        Cases[eqn, Measure[g_, dim_] -> g*dim, ∞]

Clear[Dimless]
Dimless[eqn_, 
          constants_: {}] /; (And @@ (MemberQ[eqn, #1, Infinity] & /@ 
                  constants) || 
            constants === {}) := 
Module[{relrules}, 
relrules = 
          DimlessRel[MakeRelation[eqn], 
            constants]; 
(eqn /. #1 /. {Measure[g_, 
                        dim_] -> g*dim} // Simplify) & /@ 
          relrules
]