Multilayer Diffusion
Quick Start
It's a bit complicated to take into account different thicknesses, diffusion coefficients, concentrations and partition coefficients, but it can be done, with some practical limitations.
Credits
This is a generalisation of some multilayer apps created for other specific setups.
Multilayer Diffusion
h1 μm
K1/Next
C1
h2 μm
K2/Next
C2
h3 μm
K3/Next
C3
h4 μm
K4/Next
C4
h5 μm
K5/Next
C5
D1 cm²/s
D2 cm²/s
D3 cm²/s
D4 cm²/s
D5 cm²/s
tmax d
tstep (rel)
Remaining
"use strict"
//One universal basic required here to get things going once loaded
window.onload = function () {
//restoreDefaultValues(); //Un-comment this if you want to start with defaults
Main();
};
//Main() is hard wired as THE place to start calculating when inputs change
//It does no calculations itself, it merely sets them up, sends off variables, gets results and, if necessary, plots them.
function Main() {
//Save settings every time you calculate, so they're always ready on a reload
saveSettings();
for (let i=1; i<=5; i++) checkMinMax("K"+i)
//Send all the inputs as a structured object
//If you need to convert to, say, SI units, do it here!
const inputs = {
h1: sliders.Slideh1.value,
// K1: sliders.SlideK1.value,
K1: parseFloat(document.getElementById('K1').value),
D1: parseFloat(document.getElementById('D1').value),
C1: sliders.SlideC1.value,
h2: sliders.Slideh2.value,
K2: parseFloat(document.getElementById('K2').value),
D2: parseFloat(document.getElementById('D2').value),
C2: sliders.SlideC2.value,
h3: sliders.Slideh3.value,
K3: parseFloat(document.getElementById('K3').value),
D3: parseFloat(document.getElementById('D3').value),
C3: sliders.SlideC3.value,
h4: sliders.Slideh4.value,
K4: parseFloat(document.getElementById('K4').value),
D4: parseFloat(document.getElementById('D4').value),
C4: sliders.SlideC4.value,
h5: sliders.Slideh5.value,
K5: parseFloat(document.getElementById('K5').value),
D5: parseFloat(document.getElementById('D5').value),
C5: sliders.SlideC5.value,
tmax: sliders.Slidetmax.value,
tstepf: sliders.Slidetstep.value,
};
//Send inputs off to CalcIt where the names are instantly available
//Get all the resonses as an object, result
console.time("Calc")
const result = CalcIt(inputs);
//Set all the text box outputs
document.getElementById('remaining').value = result.remaining;
//Do all relevant plots by calling plotIt - if there's no plot, nothing happens
//plotIt is part of the app infrastructure in app.new.js
if (result.plots) {
for (let i = 0; i < result.plots.length; i++) {
plotIt(result.plots[i], result.canvas[i]);
}
}
console.timeEnd("Calc")
//You might have some other stuff to do here, but for most apps that's it for Main!
}
function checkMinMax(input){
if (parseFloat(document.getElementById(input).value) < 0.001) document.getElementById(input).value=0.001
if (parseFloat(document.getElementById(input).value) > 1000) document.getElementById(input).value=1000
}
//Here's the app calculation
function CalcIt({ h1, K1, D1, C1, h2, K2, D2, C2, h3, K3, D3, C3, h4, K4, D4, C4, h5, K5, D5, C5, tmax, tstepf }) {
let h = [], D = [], K = [], C = [], layers = [], gridC = [], gridK = [], gridD = []
if (h1 > 1) { h.push(h1 * 1e-4); D.push(D1); K.push(K1); C.push(C1); layers.push(h1 / 1) }
if (h2 > 1) { h.push(h2 * 1e-4); D.push(D2); K.push(K2); C.push(C2); layers.push(h2 / 1) }
if (h3 > 1) { h.push(h3 * 1e-4); D.push(D3); K.push(K3); C.push(C3); layers.push(h3 / 1) }
if (h4 > 1) { h.push(h4 * 1e-4); D.push(D4); K.push(K4); C.push(C4); layers.push(h4 / 1) }
if (h5 > 1) { h.push(h5 * 1e-4); D.push(D5); K.push(K5); C.push(C5); layers.push(h5 / 1) }
let i, j, n = h.length, startTotal = 0
let nLayers=0
for (i=0; i < n; i++) nLayers+=layers[i]
if (nLayers < 4) {
alert("The system is too thin")
return {
plots: [],
canvas: ['canvas'],
remaining: "Error",
}
}
let DOK = true
for (i = 0; i < n; i++) {
if (isNaN(D[i]) || D[i] > 1e-7 || D[i] < 1e-16) {
DOK = false
}
}
if (! DOK) {
alert("At least one of your D values is too small, too large or not a number")
return {
plots: [],
canvas: ['canvas'],
remaining: "Error",
}
} //Add dummy layer at the start
gridC.push(C[0]); gridK.push(K[0]); gridD.push(D[0])
for (i = 0; i < n; i++) {
for (j = 0; j < layers[i]; j++) {
gridC.push(C[i]); gridK.push(K[i]); gridD.push(D[i])
startTotal += C[i]
}
}
//Add dummy layer at the end
gridC.push(0); gridK.push(K[n - 1]); gridD.push(D[n - 1])
const nGrid = gridC.length - 1
let tmpC = []
for (i = 0; i <= nGrid; i++) { tmpC[i] = gridC[i] }
const DStep = 1e-4 //One micron steps through the grid
const DStep2 = 1 / (DStep * DStep)
const dtos = 24 * 3600
const maxD = Math.max(...D)
let tinc = tstepf/1000 * DStep * DStep / maxD
let tmin = 1; if (tmax <= 10) tmin = 0.1
let tNext = tmax / 8, rbow, Col, plotData = [], lineLabels = [], myColors = [], borderWidth = [], dottedLine = []
// let incNext= tNext
for (let t = 0; t < tmax; t += tinc / dtos) {
for (i = 1; i < nGrid; i++) {
tmpC[i] += tinc * DStep2 * (gridD[i - 1] * (gridC[i - 1] - gridC[i]) * gridK[i] / gridK[i - 1] - gridD[i] * (gridC[i] - gridC[i + 1]) * gridK[i + 1] / gridK[i])
}
for (i = 0; i <= nGrid; i++) { gridC[i] = tmpC[i] }
gridC[0] = gridC[1] //Dummy grid
// if (t > incNext) {
// tinc *=1.01
// incNext+=tmax/100
// }
if (t > tNext) {
let tPts = []
for (i = 0; i < nGrid; i++) {
tPts.push({ x: i, y: 100 * gridC[i] })
};
rbow = Rainbow(t / tmax); Col = "rgb(" + rbow.r + "," + rbow.g + "," + rbow.b + ")"
plotData.push(tPts)
lineLabels.push(t.toFixed(2) + "d")
myColors.push(Col)
borderWidth.push(2)
dottedLine.push(false)
tNext += tmax / 8
}
}
let endTotal = 0
for (i = 0; i <= nGrid; i++) { endTotal += gridC[i] }
const remaining = 100 * endTotal / startTotal
//Now set up all the graphing data detail by detail.
const prmap = {
plotData: plotData, //An array of 1 or more datasets
lineLabels: lineLabels, //An array of labels for each dataset
colors: myColors, //An array of colors for each dataset
hideLegend: false,
borderWidth: borderWidth,
xLabel: 'h&μm', //Label for the x axis, with an & to separate the units
yLabel: 'Conc&%', //Label for the y axis, with an & to separate the units
y2Label: null, //Label for the y2 axis, null if not needed
yAxisL1R2: [], //Array to say which axis each dataset goes on. Blank=Left=1
logX: false, //Is the x-axis in log form?
xTicks: undefined, //We can define a tick function if we're being fancy
logY: false, //Is the y-axis in log form?
yTicks: undefined, //We can define a tick function if we're being fancy
legendPosition: 'top', //Where we want the legend - top, bottom, left, right
xMinMax: [,], //Set min and max, e.g. [-10,100], leave one or both blank for auto
yMinMax: [100, 0], //Set min and max, e.g. [-10,100], leave one or both blank for auto
y2MinMax: [,], //Set min and max, e.g. [-10,100], leave one or both blank for auto
xSigFigs: 'F0', //These are the sig figs for the Tooltip readout. A wide choice!
ySigFigs: 'F1', //F for Fixed, P for Precision, E for exponential
};
//Now we return everything - text boxes, plot and the name of the canvas, which is 'canvas' for a single plot
return {
plots: [prmap],
canvas: ['canvas'],
remaining: isNaN(remaining)? "Error":remaining.toFixed(1) + "%",
};
}
var RB = [[0, 48, 245], [0, 52, 242], [0, 55, 238], [0, 59, 235], [3, 62, 231], [9, 66, 228], [14, 69, 225], [18, 72, 221], [20, 74, 218], [22, 77, 214], [23, 80, 211], [24, 82, 207], [25, 85, 204], [25, 87, 200], [25, 90, 197], [25, 92, 193], [25, 94, 190], [25, 96, 187], [24, 99, 183], [24, 101, 180], [24, 103, 177], [23, 105, 173], [23, 106, 170], [24, 108, 167], [24, 110, 164], [25, 112, 160], [27, 113, 157], [28, 115, 154], [30, 117, 151], [32, 118, 148], [34, 120, 145], [36, 121, 142], [39, 122, 139], [41, 124, 136], [43, 125, 133], [45, 126, 130], [47, 128, 127], [49, 129, 124], [51, 130, 121], [53, 132, 118], [54, 133, 115], [56, 134, 112], [57, 136, 109], [58, 137, 106], [59, 138, 103], [60, 139, 99], [61, 141, 96], [62, 142, 93], [62, 143, 90], [63, 145, 87], [63, 146, 83], [64, 147, 80], [64, 149, 77], [64, 150, 74], [65, 151, 70], [65, 153, 67], [65, 154, 63], [65, 155, 60], [66, 156, 56], [66, 158, 53], [67, 159, 50], [68, 160, 46], [69, 161, 43], [70, 162, 40], [71, 163, 37], [73, 164, 34], [75, 165, 31], [77, 166, 28], [79, 167, 26], [82, 168, 24], [84, 169, 22], [87, 170, 20], [90, 171, 19], [93, 172, 18], [96, 173, 17], [99, 173, 17], [102, 174, 16], [105, 175, 16], [108, 176, 16], [111, 176, 16], [114, 177, 17], [117, 178, 17], [121, 179, 17], [124, 179, 18], [127, 180, 18], [130, 181, 19], [132, 182, 19], [135, 182, 20], [138, 183, 20], [141, 184, 20], [144, 184, 21], [147, 185, 21], [150, 186, 22], [153, 186, 22], [155, 187, 23], [158, 188, 23], [161, 188, 24], [164, 189, 24], [166, 190, 25], [169, 190, 25], [172, 191, 25], [175, 192, 26], [177, 192, 26], [180, 193, 27], [183, 194, 27], [186, 194, 28], [188, 195, 28], [191, 195, 29], [194, 196, 29], [196, 197, 30], [199, 197, 30], [202, 198, 30], [204, 199, 31], [207, 199, 31], [210, 200, 32], [212, 200, 32], [215, 201, 33], [217, 201, 33], [220, 202, 34], [223, 202, 34], [225, 202, 34], [227, 203, 35], [230, 203, 35], [232, 203, 35], [234, 203, 36], [236, 203, 36], [238, 203, 36], [240, 203, 36], [241, 202, 36], [243, 202, 36], [244, 201, 36], [245, 200, 36], [246, 200, 36], [247, 199, 36], [248, 197, 36], [248, 196, 36], [249, 195, 36], [249, 194, 35], [249, 192, 35], [250, 191, 35], [250, 190, 35], [250, 188, 34], [250, 187, 34], [250, 185, 34], [250, 184, 33], [250, 182, 33], [250, 180, 33], [250, 179, 32], [249, 177, 32], [249, 176, 32], [249, 174, 31], [249, 173, 31], [249, 171, 31], [249, 169, 30], [249, 168, 30], [249, 166, 30], [248, 165, 29], [248, 163, 29], [248, 161, 29], [248, 160, 29], [248, 158, 28], [248, 157, 28], [248, 155, 28], [247, 153, 27], [247, 152, 27], [247, 150, 27], [247, 148, 26], [247, 147, 26], [246, 145, 26], [246, 143, 26], [246, 142, 25], [246, 140, 25], [246, 138, 25], [245, 137, 24], [245, 135, 24], [245, 133, 24], [245, 132, 24], [244, 130, 23], [244, 128, 23], [244, 127, 23], [244, 125, 23], [244, 123, 22], [243, 121, 22], [243, 119, 22], [243, 118, 22], [243, 116, 21], [242, 114, 21], [242, 112, 21], [242, 110, 21], [241, 109, 21], [241, 107, 21], [241, 105, 21], [241, 103, 21], [240, 101, 21], [240, 100, 22], [240, 98, 22], [240, 96, 23], [240, 95, 24], [240, 93, 26], [240, 92, 27], [240, 90, 29], [240, 89, 31], [240, 88, 33], [240, 87, 36], [240, 87, 38], [241, 86, 41], [241, 86, 44], [242, 86, 47], [242, 86, 51], [243, 86, 54], [243, 87, 58], [244, 88, 62], [245, 88, 65], [245, 89, 69], [246, 90, 73], [247, 91, 77], [247, 92, 82], [248, 94, 86], [249, 95, 90], [249, 96, 94], [250, 97, 98], [251, 99, 102], [251, 100, 106], [252, 101, 111], [252, 103, 115], [253, 104, 119], [253, 105, 123], [254, 107, 128], [254, 108, 132], [255, 109, 136], [255, 111, 140], [255, 112, 145], [255, 114, 149], [255, 115, 153], [255, 116, 157], [255, 118, 162], [255, 119, 166], [255, 120, 170], [255, 122, 175], [255, 123, 179], [255, 125, 183], [255, 126, 188], [255, 127, 192], [255, 129, 196], [255, 130, 201], [255, 132, 205], [255, 133, 210], [255, 134, 214], [255, 136, 219], [255, 137, 223], [255, 139, 227], [255, 140, 232], [255, 141, 236], [254, 143, 241], [254, 144, 245], [253, 146, 250]]
function Rainbow(v) {
const i = Math.floor((Math.min(v, 1), Math.max(v, 0)) * 255)
const r = RB[i][0]
const g = RB[i][1]
const b = RB[i][2]
return { r: r, g: g, b: b }
}
What are we trying to do?
It is known that a generalised multi-layer diffusion system can only be solved via numerics. It is also known that numerics can either "blow up" or far too slow. So everything is a compromise.
The compromise here is the assumption that the 1-5 polymers can only be a maximum of 100μm thick (so these aren't mm-thick calculations), are in steps of 1μm and that diffusion coefficients are no larger that 1e-7 cm²/s (and even here, calculation times can be large), and that you start with a small timescale (maybe just a day at first) before looking at larger time-scales when calculations can get very long. When calculations blow up, the tstep slider can be used as described below.
What are the inputs?
For each of the 5 layers we need a thickness, h, an initial concentration, C, (from 0 to 1) of the molecule of interest, a partition coefficient K which is "this layer over the next layer". Then there is the diffusion coefficient, D, in units of cm²/s, usually entered as an exponential: 5e-9 means 5.10-9. Then you choose the number of days for the simulation: start small while you sort out your setup.
Some checks are made on total thickness and D values; it's better to have an error alert than a calculation that hangs your browser.
What are the outputs?
The "Remaining" value tells you what % of the original is left. Then there is the graph.
The graph is colour coded in time - from blue, short time, to red, long time. Once you're used to it, it's a very powerful way to show what is happening in your system. If, for example, you increase K for one layer, this means that the molecule prefers to stay in that layer, so you generate a big step between this layer and the next.
It's not just some lines. Move your mouse over any line and you are told the time and concentration at that point.
tstep (rel)
A large value gives you large time steps in the calculation, giving fast calculations. But under some combinations of values, the numerics blow up and you get a strange (or blank) graph and Error in the Remaining output. So now try smaller values of tstep. Don't rush - as you get smaller, calculations take longer and once they're more than a few seconds, you're out of control. If the response is getting too slow, reduce tmax. There's no guarantee that you'll find a combination for your specific setup, but in my testing it copes with a wide range of scenarios.
