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
C1
K1/Next
h2 μm
C2
K2/Next
h3 μm
C3
K3/Next
h4 μm
C4
K4/Next
h5 μm
C5
K5/Next
D1 cm²/s
D2 cm²/s
D3 cm²/s
D4 cm²/s
D5 cm²/s
tmax d
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();
//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,
D1: parseFloat(document.getElementById('D1').value),
C1: sliders.SlideC1.value,
h2: sliders.Slideh2.value,
K2: sliders.SlideK2.value,
D2: parseFloat(document.getElementById('D2').value),
C2: sliders.SlideC2.value,
h3: sliders.Slideh3.value,
K3: sliders.SlideK3.value,
D3: parseFloat(document.getElementById('D3').value),
C3: sliders.SlideC3.value,
h4: sliders.Slideh4.value,
K4: sliders.SlideK4.value,
D4: parseFloat(document.getElementById('D4').value),
C4: sliders.SlideC4.value,
h5: sliders.Slideh5.value,
K5: sliders.SlideK5.value,
D5: parseFloat(document.getElementById('D5').value),
C5: sliders.SlideC5.value,
tmax: sliders.Slidetmax.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;
// document.getElementById('Emitted1').value = result.Emitted1;
// document.getElementById('Peak2').value = result.Peak2;
// document.getElementById('Emitted2').value = result.Emitted2;
// document.getElementById('Flow').value = result.Flow;
//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!
}
//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 }) {
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 / 2) }
if (h2 > 1) { h.push(h2 * 1e-4); D.push(D2); K.push(K2); C.push(C2); layers.push(h2 / 2) }
if (h3 > 1) { h.push(h3 * 1e-4); D.push(D3); K.push(K3); C.push(C3); layers.push(h3 / 2) }
if (h4 > 1) { h.push(h4 * 1e-4); D.push(D4); K.push(K4); C.push(C4); layers.push(h4 / 2) }
if (h5 > 1) { h.push(h5 * 1e-4); D.push(D5); K.push(K5); C.push(C5); layers.push(h5 / 2) }
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 //Two micron steps through the grid
const DStep2 = 1 / (DStep * DStep)
const dtos = 24 * 3600
const maxD = Math.max(...D)
let tinc = 0.5 * 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 * 2, 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: 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 2μm (this speeds things up 4x compared to 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.
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.
