Sandor Balog |  Adolphe Merkle Institute | Université de Fribourg | Switzerland 

2022.12.24
Mathematica Version 13.0
sandor.balog@unifr.ch

This Wolfram / Mathematica notebook accompanies the article “The Influence of Fluid Menisci on Nanoparticle Dosimetry in Submerged Cells” written by Christina Glaubitz, Laetitia Haeni, Eva Sušnik, Barbara Rothen-Rutishauser, Sandor Balog, and Alke Petri-Fink

Keywords: In vitro particle dosimetry, Fluid meniscus, Mathematical modelling

The mathematical model is a continuous and deterministic model that comprises two second-order nonlinear and linear partial differential equations with corresponding boundary and initial conditions. The model addresses the impact of fluid the meniscus in cylindrical wells, and the transport equation (advection-diffusion) in three plus one dimensions: spatial coordinates and time. The advection-diffusion model is evaluated via finite element method, addressing non-interacting colloidally stable particles in a quiescent fluid.

SetDirectory[dirS]; 
dir1=CreateDirectory[fn=DateString["ISODate"]<>"_"<>DateString["Hour" ]<>"-"<>DateString["Minute"]<>"-"<>DateString["Second"]];
SetDirectory[dir1];

Print[Style["Input parameters are set and folder is created.",FontColor->Red,FontWeight->Bold,FontSize->10]];

vol= vol milli liter;
conc0=c0 ((micro g)/(milli liter));
massTot=conc0  vol;

T=273.15+TCelsius;

{a,b,c}={0.0000383897,432.9089,160.41438};
visc[T_]:=a Exp[b/(T-c)] (* Viscosity of water as a function of temperature (Pa s) *)

kB =1.38 (10^(-23));
η=visc[T];
Diff=((kB T)/(6 π η ))(1/(rh (10^(-9)) )); (* Particle translational diffusion coeff. ((m^2)/s) *)
Sv=((2gr)/(9η))(ρ-ro) (((rh (10^(-9))))^2); (* Particle settling velocity (m/s) *)

fig1=Plot[(1/Sv )√(Diff/(t 3600 )),{t,1,Tep},
PlotLabel->"Particle transport factor",
PlotRange->{{All,24+1},{0,All}},
PlotStyle->Black,

FrameLabel->{"Time (h)", "ϵ" },
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,

GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},

LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
];

If[Round[(Diff/(3600 (Sv^2)))]<Tep,
ms="Particle transport is dominated by translational diffusion until "<>ToString[Round[(Diff/(3600 (Sv^2)))]]<>" h.";
,
ms="Particle transport is dominated by translational diffusion.";
];

P={{"Particle hydrodynamic radius (nm):", 
"Particle mass density (g/(m^3)):",
"Transport: ", 
"Particle concentration (μg/mL):",
"Administered particle mass (μg): ",
"Well radius (mm):",
"Boundary model of axis-symmetric Young-Laplace equation: " ,
"Contact angle (°): ",
"Fluid volume in well (mL:",
"Fluid surface tension (N/m):",
"Fluid mass density (g/(m^3)):",
"Fluid temperature (°C):",
"Number of annuli: ",
"Particle translation diff. coeff. ((m^2)/s):", 
"Particle settling velocity (m/s):", 
"Neumann boundary value at well bottom (1/m)", 
"Time end-point (h)",
"Folder to save the results:"} ,
{1. rh ,1. ρ ,ms,1. c0,1. massTot (10^6)/g,1. R,BC,θ,1. vol/( milli liter),1. St,1. ro,1. TCelsius,1. v,1. Diff,1. Sv,1. κ,1. Tep, dir1}}//Transpose//TableForm;

Print[Style[P,FontColor->Red,FontWeight->Bold,FontSize->10]];

R=R (10^(-3));
γ=St/(ro gr);
vol=vol/(m^3);
Heq=(vol/((R^2) π)); (*  Volume-equivalent height for the 1D model (assuming flat fluid-air interface) *)
θ=θ(π/180);

Clear[Vol,zz]
If[BC<2,
Vol[R_,z0_,γ_]:=π (R^2) z0 Hypergeometric0F1Regularized[2,((R^2)/(4 γ))];
Clear[z0];
z0=NSolve[Vol[R,z0,γ]==vol,z0][[1,1,2]];
zz[r_]:=z0 BesselJ[0,((I r)/(Sqrt[γ]))];
,
Vol[R_,θ_,γ_]:=2 π R γ Cot[θ];
Clear[h];
h=((vol-Vol[R,θ,γ])/(π (R^2)));
zz[r_]:=h+((I (Sqrt[γ]) BesselJ[0,((I r)/(Sqrt[γ]))] Cot[θ])/(BesselJ[1,((I R)/(Sqrt[γ]))]));
];

Z=Interpolation[Table[{r,Chop[zz[r]]},{r,0,R,R/99}],InterpolationOrder->1];
Clear[zz];

Clear[r,z] (* Construct 2D cross-sectional mesh for solving the 3D equation *) 
Ω=ImplicitRegion[0≤z≤Z[r],{{r,0,R},{z,0,Z[R]}}];
mesh=ToElementMesh[DiscretizeRegion[Ω],MaxCellMeasure->(10^(-8))]; 
meshX=ToElementMesh[DiscretizeRegion[Ω],MaxCellMeasure->Automatic]; 

fig2=Show[
fig2a=Plot[Z[r],{r,0,R},
PlotRange->{{0,R},{0,All}},

PlotLabel->"Well: half cross-sectional view",
PlotStyle->{Blue},
Filling->Bottom,
FrameLabel->{"r (m)","z (m)"},
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->Automatic
]
,
meshX["Wireframe"["MeshElementStyle"->EdgeForm[{Thickness[0.00001],Blue}]]]
];

Clear[meshX];

X=Table[√(x/v) R,{x,1 ,v}];
fig3=Graphics[{Table[Circle[{0,0},X[[i]]],{i,v}]},ImageSize->5 cm]

Print[Style["Particle-transport variables are computed and well geometry and mesh are constructed. Computation begins.",FontColor->Red,FontWeight->Bold,FontSize->10]];

Clear[u];
fig4=Plot[{(κ u)/(1+κ u),κ u},{u,0,c0 (10^3)},
PlotRange->{{0,All},{0,All}},

PlotLabel->"Hill-Langmuir model",
PlotStyle->{Gray,{Orange,Dashed}},
FrameLabel->{"Particle concentration near cell surface (μg/mL)", "Neumann boundary value"},
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
];

tlist=Table[i,{i,0,Tep,1}];

Clear[u1,t,x];
op1=(∂_t)u1[t,x]+Sv (∂_x)u1[t,x]-Diff (∂_x)(∂_x)u1[t,x]; (* Advection-diffusion equation for the 1D model *)

tmax =tlist[[-1]];
Clear[u1,t,x,sol1Dirichlet];  (* 1D solution wiht perfectly adsorbing/absorbing bottom, Dirichlet boundary *)
sol1Dirichlet=NDSolve[{op1==0,DirichletCondition[u1[t,x]==0,x==0],u1[0,x]== c0},u1, {t,0,tmax 3600},{x,0, Heq }][[1]];

Clear[u1,t,x,sol1];  (* 1D solution wiht conc.-dependent affintify for adsorption/absorption, Neumann boundary *)
sol1=NDSolve[{op1==NeumannValue[-κ u1[t,x],x==0],u1[0,x]== c0},u1, {t,0,tmax 3600},{x,0,Heq }][[1]];

t=tlist 3600;
d1Dirichlet=(10^3)((massTot/g- (R^2)π NIntegrate[Evaluate[u1[t,x]/.sol1Dirichlet],{x,0, Heq}])/((R^2)π));
D1Dirichlet=Partition[Riffle[tlist,d1Dirichlet],2];

d1=(10^3)((massTot/g- (R^2)π NIntegrate[Evaluate[u1[t,x]/.sol1],{x,0, Heq}])/((R^2)π));
D1=Partition[Riffle[tlist,d1],2];

fig5=Show[
Plot[(10^3)(massTot/g/((R^2)π)),{t,0,tmax},
PlotRange->{{0,tmax+1},{0,(10^3)1.05 (massTot/g/((R^2)π))}},
PlotLabel->"Solution of 1D model",
PlotStyle->Gray,
FrameLabel->{"Time (h)", "Mass surface density (mg (m^(-2)))"},
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
]
,
ListLinePlot[{D1Dirichlet,D1},PlotRange->{0,All},
PlotMarkers->None,
PlotStyle->{Black,{Red,Dashed}},
LabelStyle->ls
,
PlotLegends->Placed[LineLegend[{"1D Dirichlet","1D Neumann"},
LegendLayout->{"Column",1},
LegendFunction->(Framed[#,RoundingRadius->0,FrameStyle->None,Background->White]&)],{0.3,0.75}]
]
];

Clear[u1,t,x];
op1=(∂_t)u1[t,x]+Sv (∂_x)u1[t,x]-Diff (∂_x)(∂_x)u1[t,x]; (* Advection-diffusion equation for the 1D model *)
op1A=(∂_t)u1[t,x]+(10^(-6)) Sv (∂_x)u1[t,x]-Diff (∂_x)(∂_x)u1[t,x]; (* Advection-diffusion equation for the 1D model *)
op1B=(∂_t)u1[t,x]+Sv (∂_x)u1[t,x]-(10^(-6)) Diff (∂_x)(∂_x)u1[t,x]; (* Advection-diffusion equation for the 1D model *)

tmax =tlist[[-1]];
Clear[u1,t,x,sol1Dirichlet];  (* 1D solution wiht perfectly adsorbing/absorbing bottom, Dirichlet boundary *)
sol1D=NDSolve[{op1==0,DirichletCondition[u1[t,x]==0,x==0],u1[0,x]== c0},u1, {t,0,tmax 3600},{x,0, Heq }][[1]];
Clear[u1,t,x,sol1Dirichlet];  (* 1D solution wiht perfectly adsorbing/absorbing bottom, Dirichlet boundary *)
sol1DA=NDSolve[{op1A==0,DirichletCondition[u1[t,x]==0,x==0],u1[0,x]== c0},u1, {t,0,tmax 3600},{x,0, Heq }][[1]];
Clear[u1,t,x,sol1Dirichlet];  (* 1D solution wiht perfectly adsorbing/absorbing bottom, Dirichlet boundary *)
sol1DB=NDSolve[{op1B==0,DirichletCondition[u1[t,x]==0,x==0],u1[0,x]== c0},u1, {t,0,tmax 3600},{x,0, Heq }][[1]];

t=tlist 3600;
AB=Partition[Riffle[tlist,(10^3)((massTot/g- (R^2)π NIntegrate[Evaluate[u1[t,x]/.sol1D],{x,0, Heq}])/((R^2)π))],2];
A=Partition[Riffle[tlist,(10^3)((massTot/g- (R^2)π NIntegrate[Evaluate[u1[t,x]/.sol1DA],{x,0, Heq}])/((R^2)π))],2];
B=Partition[Riffle[tlist,(10^3)((massTot/g- (R^2)π NIntegrate[Evaluate[u1[t,x]/.sol1DB],{x,0, Heq}])/((R^2)π))],2];

fig6=Show[
Plot[(10^3)(massTot/g/((R^2)π)),{t,0,tmax},
PlotRange->{{0,tmax+1},{0,(10^3)1.05 (massTot/g/((R^2)π))}},
PlotLabel->"Solution of 1D model / Dirichlet boundary",
PlotStyle->Gray,
FrameLabel->{"Time (h)", "Mass surface density (mg (m^(-2)))"},
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
]
,
ListLinePlot[{AB,A,B},PlotRange->{0,All},
PlotMarkers->None,
PlotStyle->{Black,{Red,Dashed},{Blue,Dashed}},
LabelStyle->ls
,
PlotLegends->Placed[LineLegend[{"Diffusion & Settling","Diffusion only","Settling only"},
LegendLayout->{"Column",1},
LegendFunction->(Framed[#,RoundingRadius->0,FrameStyle->None,Background->None]&)],{0.5,0.75}]
]
];

Print[Style["1D model is ready. Computing the 3D model begins.",FontColor->Red,FontWeight->Bold,FontSize->10]];

Clear[u,t,r,z];  (* Advection-diffusion equation for the 3D model *)
op = (∂_t)u[t,r,z]-Inactive[Div][Inactive[Plus][Times[Diff,Inactive[Grad][u[t, r, z], {r, z}]],Times[{0, Sv}, u[t, r, z]]], {r, z}]-1*(1/r)(∂_r) Diff u[t,r,z];
(* Computes particle concentration in well*)
sol=NDSolve[{op==NeumannValue[-κ u[t,r,z],z==0],u[0,r,z]==c0},u[t, r, z],{t,0, 3600tmax},{r,z}∈mesh,Method->{"PDEDiscretization"->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement"}}},InterpolationOrder->Automatic][[1]];

Clear[figx,figy,figz]
j=1;ProgressIndicator[Dynamic[j],{0,Length[tlist]}]
While[j≤Length[tlist],

t=tlist[[j]] 3600;

If[j>1,

figx[j]=Plot3D[Evaluate[Log[10,u[t,r,z]/.sol]],{r,z}∈mesh,
PlotLabel->"Time: "<>ToString[Round[N[t/3600 ]]]<>" h",
PlotRange->{0,All},(*c0*)
AxesLabel->{"r (m)","z (m)","c (μg (mL^(-1)))"},
LabelStyle->{FontColor->Black,FontWeight->"Light",FontFamily->"Arial Narrow",FontSize->7},
ImagePadding->Full,
(*ImageSize->2{R  (10^3)cm,Z[R] (10^3)cm},*)
ImageSize->{sx,sy},
ImageMargins->0. cm,
BoxRatios->{R,Z[R],0.4 R},
ColorFunction-> "TemperatureMap",
BoxStyle->fs
];

figy[j]=Show[
fig2a
,
DensityPlot[Evaluate[Log[10,u[t,r,z]/.sol]],{r,z}∈mesh,
(*PlotLabel->"Time: "<>ToString[Round[N[t/3600 ]]]<>" h",*)
PlotRange->All,(* c0 *)
FrameLabel->{"r (m)","z (m)","(Subscript[Log, 10]) c (μg (mL^(-1)))"},
LabelStyle->{FontColor->Black,FontWeight->"Light",FontFamily->"Arial Narrow",FontSize->12},
ImagePadding->Full,
ColorFunction-> "TemperatureMap",
PlotLegends->None
]
];

Export["fig_"<>ToString[tlist[[j]] ]<>"h.gif",GraphicsRow[{figx[j],figy[j]},ImageSize->{2sx, sy}],ImageResolution->300];

];

j++;
];
Clear[j];

(* Computes total dose, metric: mass / surface area *)
D3={}; ProgressIndicator[Dynamic[i],{1,Length[tlist]}]
i=1;While[i≤Length[tlist],
t=tlist[[i]] 3600;
Q={};
r=0;
While[r≤R,
Zr=Table[{z,2 r π Evaluate[u[t,r,z]/.sol]},{z,0,Z[r],0.005 Z[r]}];
ipZr=Interpolation[Zr,InterpolationOrder->1];
AppendTo[Q,{r,NIntegrate[ipZr[z],{z,0,Z[r]}]}];
r+=0.025R
];
ip=Interpolation[Q,InterpolationOrder->2];

AppendTo[D3,{t/3600,(massTot/g-NIntegrate[ip[r],{r,0,R}])/((R^2)π)}];
i++];

Clear[Q,ip,Zr,ipZr];

mD3=D3; mD3[[;;,2]]=(10^3)D3[[;;,2]];

fig7=
Show[
ListLinePlot[{mD3,D1},
PlotRange->{{0,tmax+1},{0,(10^3)1.05(massTot/g/((R^2)π))}},
PlotLabel->"Solution of 1D and 3D models",
PlotStyle->{RGBColor[0.2,0.7,1],{Red,Dashed}},
PlotMarkers->Automatic,
FrameLabel->{"Time (h)", "Mass surface density (mg (m^(-2)))"},
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
,
PlotLegends->Placed[LineLegend[{"3D-model","1D-model"},
LegendLayout->{"Column",1},
LegendFunction->(Framed[#,RoundingRadius->0,FrameStyle->None,Background->None]&)],{0.3,0.75}]
]
,
Plot[(10^3)(massTot/g/((R^2)π)),{t,0,tmax+1},PlotRange->All,PlotStyle->{Dashed,Black}]
];

Clear[t](* Computes total mass / g *)
ipM =Interpolation[Partition[Riffle[ D3[[;;,1]],(R^2)π D3[[;;,2]]],2],InterpolationOrder->1];

(*Plot[(10^6)ipM [t],{t,D3[[1,1]],D3[[-1,1]]},
PlotRange->{{0,All},{0,massTot/g (10^6)}},
PlotStyle->Black,
FrameLabel->{"Time (h)", "m (μg)" },
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
]*)

Δj=0.1;
L=tmax/Δj;

l=1; ΔM={};
While[l≤L,
AppendTo[ΔM,(ipM[l Δj ]-ipM[(l-1) Δj ])];
l++;
];

(*ListPlot[ΔM,
PlotRange->{{0,All},{0,All}},
PlotStyle->{{Gray,Thin}},
PlotMarkers-> {"OpenMarkers",5},
Joined->True,
InterpolationOrder->1,
FrameLabel->{"Time (h)", "dm (t-Δt,t) (g)"   },
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
]
*)
l=1;Clear[γm];
While[l≤L,
t=l Δj 3600;

Γ={};
Clear[r];
z=0;
AppendTo[Γ,Mean[Table[Evaluate[u[t,r,z]]/.sol,{r,0 ,X[[1]],X[[1]]/10}]]];
i=2;
While[i≤v,
AppendTo[Γ,Mean[Table[Evaluate[u[t,r,z]]/.sol,{r,X[[i-1]] ,X[[i]],(X[[i]]-X[[i-1]])/10}]]];
i++
];

ΔM=(ipM[l Δj ]-ipM[(l-1) Δj ]);
γm[l]=(ΔM/Total[Γ])Γ;
l++;
];

ps=Partition[Riffle[Table[Thin,{i,L}],Table[RGBColor[i,0.25,0.85],{i,0,1,1/(L-1)}]],2];

(*fig6=ListLogPlot[Table[γm[l],{l,L}],
PlotRange->{{1-0.1 v,1.1v},All},
PlotStyle->ps,
PlotMarkers-> {"OpenMarkers",3},
Joined->True,
InterpolationOrder->1,
FrameLabel->{"Area", "dm (g)"   },
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Table[i,{i,v}],None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
]*)

l=1;Clear[δμ];
While[l≤L,
s=Total[Table[γm[y],{y,l}]];
δμ[l]=s;
l++;
];

(*fig7=ListLinePlot[Table[δμ[l],{l,L}],
PlotRange->{{1-0.1 v,1.1v},{0,All}},
PlotStyle->ps,
PlotMarkers-> {"OpenMarkers",3},
Joined->True,
InterpolationOrder->1,
FrameLabel->{"Area", "m (g)"   },
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Table[i,{i,v}],None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx
]*)

fig8=ListLinePlot[Table[(10^3)(v/((R^2)π ))δμ[l],{l,1,L,10}],
PlotLabel->"Mass density averaged over the annuli",
PlotRange->{{1-0.1 v,1.1v},{0,All}},
Joined->True,
InterpolationOrder->1,
FrameLabel->{"Annulus", "M (mg (m^(-2)))"   },
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Table[i,{i,v}],None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx,
ColorFunction-> "TemperatureMap",
PlotStyle->Thick
];

fig9=Show[
ListLinePlot[Table[{Mean[(1/((R^2)π/v ))(10^3)δμ[l]],100((δμ[l][[-1]]-δμ[l][[1]])/δμ[l][[1]])},{l,L}],
PlotLabel->"Relative mass density difference on the annuli",
PlotRange->{{0,All},All},
PlotStyle->{Thickness[0.02]},
PlotMarkers->None,
Joined->True,
InterpolationOrder->1,
FrameLabel->{"M (mg (m^(-2)))", "(Subscript[Δ, M]) (%)"   },
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx,
ColorFunction-> "TemperatureMap"
]
,
ListLinePlot[Table[{Mean[(10^3)(1/((R^2)π/v ))δμ[l]],100((δμ[l][[-1]]-δμ[l][[1]])/δμ[l][[1]])},{l,1,L}],
PlotStyle->{Black,Thin,Dashed},
PlotRange->All
]
];

Print[Style["3D model is ready. Saving the results begins.",FontColor->Red,FontWeight->Bold,FontSize->10]];
SystemOpen[dir1];

χ=3.;

D1=1. (10^(-χ))Floor[(10^χ) D1];
mD3=1. (10^(-χ))Floor[(10^χ) mD3];
p1=Table[{D1[[i,1]],D1[[i,2]],mD3[[i,2]]},{i,2,Length[D1]}];
Export["Total_dose.csv",N[p1],"CSV"]
p1=PrependTo[p1,{"Time (h)","1D (mg (m^(-2)))","3D (mg (m^(-2)))"}];
Export["Total_dose.pdf",N[p1]//TableForm];

m1=Table[tlist[[i]],{i,2,Length[tlist]}];
m2=Table[(10^3)(v/((R^2)π ))δμ[l],{l,1,L,10}];
m={};
AppendTo[m,m1];
i=1;
While[i<=v,
AppendTo[m,m2[[;;,i]]];
i++
];

fig10=ListPlot[Table[Partition[Riffle[m[[1,;;]],m[[i+1,;;]]],2],{i,1,v}],
PlotLabel->"Mass density averaged over the annuli vs time",
PlotRange->{{0,All},{0,All}},
Joined->True,
InterpolationOrder->1,
FrameLabel->{"Time (h)", "M (mg (m^(-2)))"   },
Frame->True,
FrameStyle->{fs,fs,fs,fs},
FrameTicks->{{Automatic,None},{Automatic,None}},
FrameTicksStyle->th,
GridLines->None,
GridLinesStyle->{{th,Dashed},{th,Dashed}},
LabelStyle->ls,
ImagePadding->Full,
ImageSize->{sx,sy},
ImageMargins->0. cm,
AspectRatio->sy/sx,
ColorFunction-> Automatic,
PlotStyle->Thick
];

m=1. (10^(-χ))Floor[(10^χ) m];
m=m//Transpose;
m//TableForm;
Export["Dose_on_annuli.csv",N[m],"CSV"];
m=PrependTo[m,{"Time (h)",Table["Annulus "<>ToString[i]<>" (mg (m^(-2)))",{i,v}]}//Flatten]//TableForm;
Export["Dose_on_annuli.pdf",N[m]];

p1=Table[{Mean[(1/((R^2)π/v ))(10^3)δμ[l]],100((δμ[l][[-1]]-δμ[l][[1]])/δμ[l][[1]])},{l,1,L,10}]//Abs;
p1=1. (10^(-χ))Floor[(10^χ) p1];
Export["Dose_bias.csv",N[p1],"CSV"];
p1=PrependTo[p1,{"Total dose 3D (mg (m^(-2)))","Rel. diff. %"}];
Export["Dose_bias.pdf",N[p1]//TableForm];
Export["Parameters_"<>ToString[fn]<>".pdf",P]

Export["Annuli radii.csv",N[X],"CSV"];

j=1;
While[j≤10,
Export["fig"<>ToString[j]<>".gif",ToExpression["fig"<>ToString[j]],ImageResolution->300];
j++];

Print[Style["Ready!",FontColor->Red,FontWeight->Bold,FontSize->10]];
SystemOpen[dir1]
FrontEndExecute[FrontEndToken["DeleteGeneratedCells"]]