awPowrWB.nb
Water Activity Power-Law Model with Wet Basis Input

This Mathematica notebook estimates the equilibrium water activity, awStar, of a dry multicomponent food mixture stored in a hermetically sealed container. It is assumed that the ingredients do not interact chemically, that the temperature is constant and that the amount of moisture absorbed by or released to the atmosphere trapped in the container is negligible.  The ingredients' moisture sorption isotherms are described by nEqn power-law equations each containing nTrm terms.  Input moisture values should be entered on a wet basis and will be converted to a dry basis for calculation and output.

Written by:
Mark D. Normand                 and                 Micha Peleg
EDP Programmer II                                      Professor
Tel: 413-545-2365                                       Tel: 413-545-5852
Fax: 413-545-1262                                       Fax: 413-545-1262
E-mail: (MDN) see http://www-unix.oit.umass.edu/~adva000/
             (MP)    see http://www-unix.oit.umass.edu/~aew2000/
of the
Department of Food Science
Chenoweth Laboratory
University of Massachusetts
Amherst, MA  01003-1410

Last modified: April 16, 2004

User Input.

To apply the model  to a different multicomponent food mixture, changes should be made to any variables in User Input cells that need updating and then the notebook should be reevaluated.

Assign Values.

Define the number of equations (nEqn, one for each ingredient) and the number of terms (nTrm) in each equation.

nEqn = 3;

nTrm = 2;

Each initWetWt[i] is the initial wet weight in grams (including the water) of the ith ingredient in the mixture.

Array[initWetWt, nEqn];

initWetWt[1] = 8.8;

initWetWt[2] = 5.25;

initWetWt[3] = 23;

Each initMoistWetB[i] is the initial moisture content on a wet basis (expressed as grams of water per 100 grams of total weight) of the ith ingredient in the mixture.

Array[initMoistWetB, nEqn];

initMoistWetB[1] = 9.090899999999999;

initMoistWetB[2] = 4.761899999999999;

initMoistWetB[3] = 13.04349999999999;

Define and initialize an array to hold the size of the dashing of the line showing the moisture sorption isotherm curve for each ingredient.  Note: A dashing value of 0 produces a solid line in Mathematica version 3.0.1 on the Macintosh but causes a freeze-up in Mathematica version 2.2.2.1.  Use a very small positive value such as .001 to produce a solid line instead.

Array[dashSize, nEqn];

dashSize[1] = 0.001;

dashSize[2] = 0.01;

dashSize[3] = 0.02;

Define and initialize an array to hold the RGB color (possible values are from 0 (black) to 1 (white)) of the line showing the moisture sorption isotherm curve for each ingredient.  Example colors: {0,0,0}=black, {1,1,1}=white, {1,0,0}=red, {0,1,0}=green, {0,0,1}=blue, {0,1,1}=cyan, {1,0,1}=magenta, {1,1,0}=yellow, {.5,.5,.5}=medium gray.

Array[lineRGB, nEqn, 3];

lineRGB[1, 1] = 1; lineRGB[1, 2] = 0; lineRGB[1, 3] = 0;

lineRGB[2, 1] = 0; lineRGB[2, 2] = 1; lineRGB[2, 3] = 0;

lineRGB[3, 1] = 0; lineRGB[3, 2] = 0; lineRGB[3, 3] = 1;

Define and initialize arrays to hold the title text for each ingredient's moisture sorption isotherm graph and its initial & equilibrium points graph.

Array[mSITitle, nEqn];

Array[initEqTitle, nEqn];

Do[mSITitle[i] = "Ingredient " <> 
          ToString[i] <> " Moisture Sorption Isotherm"; 
      initEqTitle[i] = "Ingredient " <> 
          ToString[i] <> " Init. & Eq. Moisture", {i, nEqn}]

Each c[i,j] is the coefficient of the jth term of the power-law expression of the moisture sorption isotherm for the ith ingredient in the mixture.

Array[c, {nEqn, nTrm}];

c[1, 1] = 32; c[1, 2] = 25;

c[2, 1] = 14; c[2, 2] = 91;

c[3, 1] = 21; c[3, 2] = 16;

Each p[i,j] is the power of the jth term of the power-law expression of the moisture sorption isotherm for the ith ingredient in the mixture.

Array[p, {nEqn, nTrm}];

p[1, 1] = 0.53; p[1, 2] = 15;

p[2, 1] = 0.39; p[2, 2] = 6;

p[3, 1] = 0.7; p[3, 2] = 5.5;

Define the Power-Law Model.

Each mSIso[i] is the power-law expression of the moisture sorption isotherm for the ith ingredient in the mixture.

Array[mSIso, nEqn];

Do[mSIso[i][x_] = (∑_(j = 1)^nTrm) c[i, j] (x^p)[i, j]; 
      Print["mSIso[", i, "]=", mSIso[i][x]], {i, 
        nEqn}]

Model Calculations.

Convert entered wet weights to dry so that each initDryWt[i] is the initial dry weight in grams (dry matter only) of the ith ingredient in the mixture.

Array[initDryWt, nEqn];

Do[initDryWt[i] = (1/100) initWetWt[i] (100 - initMoistWetB[i]); 
      Print["The initial dry weight for ingredient ", i, " is ", 
        initDryWt[i], " grams."], {i, nEqn}]

Convert entered wet basis initial moisture contents to a dry basis so that each initMoistDryB[i] is the initial moisture content (expressed as grams of water per 100 grams of dry matter) of the ith ingredient in the mixture.

Array[initMoistDryB, nEqn];

Do[initMoistDryB[
          i] = ((100 initMoistWetB[i])/(100 - initMoistWetB[i])); 
      Print["The initial moisture content on a dry basis for ingredient \
", i, " is ", 
        initMoistDryB[
          i], " grams of water per 100 grams of dry matter."], {i, 
        nEqn}]

Calculate initMoistMixDryB, the total initial moisture content of the multicomponent food mixture on a dry basis.

initMoistMixDryB = (((∑_(i = 1)^nEqn) initDryWt[i] initMoistDryB[i])/((∑_(i = 1)^nEqn) initDryWt[i]));

Print["The total initial moisture content of the mixture on a dry \
basis is ", 
      N[initMoistMixDryB, 
        3], " grams water per 100 grams dry matter."]

Calculate awStar, the equilibrium water activity (dimensionless) of the multicomponent food mixture, by finding the root of its model equation using an initial guess (guessAw) and then extracting the final value of guessAw from the rule that FindRoot returns.

awStar = 
        guessAw /. FindRoot[
            initMoistMixDryB == (((∑_(i = 1)^nEqn) initDryWt[i] mSIso[i][guessAw])/((∑_(i = 1)^nEqn) initDryWt[i])), {guessAw, 
              0.25}, MaxIterations -> 100];

Print["The equilibrium water activity of the mixture is ", 
      N[awStar, 3]]

Calculate each aw[i], the initial water activity value for the ith ingredient, by solving its isotherm equation using its initial moisture content.

Array[aw, nEqn];

Do[aw[i] = 
        x /. FindRoot[mSIso[i][x] == 
              initMoistDryB[i], {x, 0.5}]; 
      Print["The initial water activity value for ingredient ", 
        i, " is ", N[aw[i], 4]], {i, nEqn}]

Calculate each eqMoist[i], the final (equilibrium) moisture value for the ith ingredient, by substituting the multicomponent equilibrium water activity value in each ingredient's isotherm equation.

Array[eqMoist, nEqn];

Do[eqMoist[i] = mSIso[i][awStar]; 
      Print["The equilibrium moisture value for ingredient ", 
        i, " is ", 
        N[eqMoist[i], 4], " grams water per 100 grams dry matter."], {i, 
        nEqn}]

Calculate each changeMoist[i], the ith ingredient's change in moisture.  The result is positive if moisture was gained and negative if it was lost.

Array[changeMoist, nEqn];

Do[changeMoist[i] = eqMoist[i] - initMoistDryB[i]; 
      Print["The change in moisture for ingredient ", i, " is ", 
        N[changeMoist[i], 
          4], " grams water per 100 grams dry matter."], {i, 
        nEqn}]

Plot the Results.

Define arrays to hold, for each ingredient, the moisture sorption isotherm plot, the initial & equilibrium points plot and the plot that combines those two plots with the vertical line showing the value of the equilibrium water activity of the mixture.

Array[mSI, nEqn];

Array[initEqPts, nEqn];

Array[ingredientPlot, nEqn];

Choose a realistic upper limit on the moisture values for each ingredient then find the maximum among them and use that value to set the maximum y-axis value used in subsequent graphs.

awLimit = 0.8;

Array[moistLimit, nEqn];

Do[moistLimit[i] = mSIso[i][awLimit]; 
      Print[moistLimit[i]], {i, nEqn}]

maxMoist = moistLimit[1];

Do[maxMoist = Max[moistLimit[i], maxMoist], {i, nEqn}]

maxMoist = Ceiling[maxMoist]

maxMoist = 100

Load the standard Mathematica package for plotting several data sets on a single graph using different symbols for each set.

<< "Graphics`MultipleListPlot`"

Plot the value of the root (the equilibrium water activity of the mixture) as a vertical line on a graph with the appropriate axis limits.  Suppress output of the plot until it is actually needed.

awsLine = 
      ListPlot[{{awStar, 0}, {awStar, maxMoist}}, 
        PlotRange -> {{0, 1}, {0, maxMoist}}, Frame -> True, 
        AxesLabel -> None, PlotJoined -> True, 
        RotateLabel -> False, 
        FrameLabel -> {"Aw", "Moisture"}, 
        PlotLabel -> "Moisture Sorption Isotherms of all \
ingredients", DisplayFunction -> Identity]

For each ingredient, plot its moisture sorption isotherm curve.  Next, plot the initial (open triangle) & equilibrium (solid square) moisture vs water activity points of that ingredient. Finally, combine the two plots into a third plot that also contains a vertical line marking the equilibrium water activity of the mixture.  Actually output only the combined plot for each ingredient.

Do[mSI[i] = 
        Plot[mSIso[i][x], {x, 0, 1}, 
          PlotRange -> {{0, 1}, {0, maxMoist}}, Frame -> True, 
          AxesLabel -> None, RotateLabel -> False, 
          PlotStyle -> {{Dashing[{dashSize[i]}], 
                RGBColor[lineRGB[i, 1], lineRGB[i, 2], lineRGB[i, 3]]}}, 
          FrameLabel -> {"Aw", "Moisture"}, 
          PlotLabel -> mSITitle[i], DisplayFunction -> Identity]; 
      initEqPts[i] = 
        MultipleListPlot[{{aw[i], initMoistDryB[i]}}, {{awStar, eqMoist[i]}}, 
          SymbolShape -> {PlotSymbol[Triangle, Filled -> False], 
              PlotSymbol[Box, Filled -> True]}, Frame -> True, 
          AxesLabel -> None, RotateLabel -> False, 
          PlotRange -> {{0, 1}, {0, maxMoist}}, 
          FrameLabel -> {"Aw", "Moisture"}, 
          PlotLabel -> initEqTitle[i], 
          DisplayFunction -> Identity]; 
      ingredientPlot[i] = 
        Show[{initEqPts[i], awsLine, mSI[i]}, 
          DisplayFunction -> $DisplayFunction], {i, nEqn}]

Build a list of plots showing the vertical equilbrium  water activity line and the moisture sorption isotherm curves of all ingredients with each ingredient's initial & equilibrium moisture points marked.

allPlots = {awsLine};

Do[allPlots = Append[allPlots, initEqPts[i]]; 
      allPlots = Append[allPlots, mSI[i]], {i, nEqn}]

The final plot combines all plots in the list and shows which ingredients absorb or loose moisture.  The vertical line marks the equilibrium water activity of the mixture.  For each ingredient, two points are indicated on its moisture sorption isotherm curve: the triangle identifies its initial and the square identifies its final, equilibrium moisture vs water activity point.

mixturePlot = 
      Show[allPlots, DisplayFunction -> $DisplayFunction]

End of notebook awPowrWB