Estimation of the concentration boundary layer adjacent to a flat surface using computational fluid dynamics

.


Introduction
In modern drug discovery, molecular lipophilicity is a key factor in the success of an oral drug development program (Arnott and Planey, 2012).This has influenced how active pharmaceutical ingredients and their corresponding dosage forms are evaluated during the preclinical stage where lead candidate selections are made.Molecules with high lipophilicity intrinsically have poor dissolution characteristics in water, meaning that compendia methods (standardized methods without an absorption component) are not adequate for evaluating dissolution performance and understanding the full impact of the formulation on systemic performance.As a result, academic and industrial research continues to push towards more biorelevant dissolution tests (Yukuyama et al., 2023;Keemink et al., 2019).This includes efforts to integrate absorption kinetics into the dissolution experiment by including the use of synthetic (Hedge et al., 2021;Berben et al., 2018;Hate et al., 2019;Sinko et al., 2020) or cell-based membranes (Keemink et al., 2019;Keemink et al., 2022;Li et al., 2019;Kono et al., 2019), and use non-sink dissolution conditions (Sun et al., 2016).
Permeation modeling is an important consideration in modern drug development with the inclusion of absorption kinetics and non-sink dissolution conditions.The pharmaceutical sciences are well acquainted with the film model for diffusion and this model underpins one of the most common ways to describe permeation performance, i.e. the apparent permeability (P app ).Predicting the P app requires knowing the thickness of the concentration boundary layer (CBL), because P app is the series sum of the resistances to diffusion with respect to the absorbing surface.Highly lipophilic drugs tend to be less transport-rate-limited by the hydrophobic cell membrane, thus hydrodynamic resistances begin to dominate the permeability term.While predictions of the CBL thickness surrounding particles-of importance to drug dissolution-is common practice (Chen et al., 2017;Sugano, 2008;Andersson et al., 2022;Salehi et al., 2020;Wang and Brasseur, 2019;Salehi et al., 2019), further investigation is needed for the geometry-specific hydrodynamic effects on permeation in modern dissolution-permeation (D/P) devices at the absorbing membrane surface.As it is uncommon that a detailed analysis of the hydrodynamic flow and accompanying CBL at the solid-liquid interface is performed or used in interpreting permeation data for devices that assess drug and drug product performance.
The prediction of fluid flow and mass transport conditions associated with biorelevant experiments is a natural continuation in designing and refining biorelevant methods.Two dimensionless quantities, the Schmidt number and Reynolds number, are classically used to characterize fluid and mass transport.The Schmidt number is the ratio of the momentum diffusivity over the mass diffusivity and is defined as Sc = ν/ D where ν is the kinematic viscosity and D is the diffusion coefficient of the solute.The Reynolds number is the ratio of the inertial forces and viscous forces, and is defined as Re = UL/ν where U is the fluid velocity and L is a characteristic length scale dependent upon the geometry of the problem.Literature on predicting the CBL for D/P devices is limited to classical models with simple geometries (Sugano, 2010), but prediction of CBL for more complex geometries is critical for understanding mass transfer in contemporary preclinical models.
The United States Pharmacopeia is a congress that establishes standards for executing a dissolution experiment, as well as for production of equipment to perform these experiments (Marroum, 2014).No similar scientific congress exists for permeation sciences.As a result, there is a wide variety of permeation apparatuses with complex geometries in preclinical settings.The increased accessibility to high performance computing resources means simulation of pharmaceutical problems is becoming easier to perform (Rantanen and Khinast, 2015;Vulović et al., 2018;Hopgood et al., 2018;Rahman et al., 2022).
A gap exists in the literature for standardized analysis of the CBL for the many D/P system types used in pharmaceutical laboratories.Johansson et al. (2018) used CFD to characterize a widely available, side-by-side diffusion cell used for academic and industrial pharmaceutical characterization of drug molecules and drug products (Johansson et al., 2018).The authors' model predicts the fluid velocity under a variety of factors such as the shape of the mixer and the position of in situ probes used to assay concentration.A novel aspect of their report is the use of the Péclet number, Pe = Re⋅Sc, to estimate the thickness of the CBL.This is done without simulating the concentration transport for the agitator that contains a cylindrical compact of drug in a magnetically driven rotating cylinder.The advantage of this approach is that laboratories can rapidly characterize their, often, custom-built D/P devices.This could possibly allow harmonization of D/P experimental methods or increase the translation of permeation results in between laboratories.
The aim of the work described herein was to investigate the validity of Johansson et al.'s method.To verify the Johansson et al. method we performed a benchmark simulation for flow over a flat plate using MATLAB (Fig. 1) and CFD (COMSOL Multiphysics v6.1, Fig. 2), both of which calculated the coupled flow field and concentration transport.We then applied two methods to estimate the CBL to the COMSOL solved flow (velocity) field.The laminar flow over a flat plate is a classical fluid dynamical flow case for which Blasius solution provides analytical solution for the velocity distribution.The Blasius solution was used to validate the computed flow field in COMSOL, after which the concentration transport for flow over a flat plate was assessed.The resulting CBL lengths via a Schmidt number correlation was compared to predictions made by Johansson et al.'s method.A desired outcome of this work was to present the possibility to apply CFD-based methods like the one presented herein to more complex geometries which represent dissolution-permeation apparatuses currently in-use for preclinical development of oral drug formulations.

The Blasius solution for boundary layer flow and coupled concentration transport
Evolution of the concentration field is intrinsically connected to evolution of the velocity field through the ratio of momentum diffusivity and mass diffusivity.This section describes the procedure to solve the coupled 2D incompressible Navier-Stokes equations and convection-diffusion equation in MATLAB to generate a reference solution.These coupled equations were used to calculate the velocity field (fluid flow), capture the boundary layer flow, and simultaneously model concentration transport (Batchelor, 2000).The detailed equations solved are found in Supplemental Section A.
Both the Navier-Stokes and convection-diffusion equation were nondimensionalized and transformed from partial differential equations thereby forming series of ordinary differential equations compatible with the ode45 ordinary differential equation solver in MATLAB 2019a.Boundary conditions, in non-dimensional form, for these equations took the form of: , where f, f ' , f '' , g, and g ' are the nondimensionalized functions related to the stream function (f), velocity (f ′ ), shear stress (f ″ ), concentration (g), and concentration gradient (g ′ ).ode45 was used to numerically evaluate the entire set of equations and was iteratively solved using the bi-section method.The bisection  method is a root-finding method used to iteratively solve the ode45 solver for the desired far-field values to a specified tolerance.The solution was achieved when the total error in the flow and concentration equations was less than or equal to 1 × 10 − 6 , or reached 100 iterations.If the solver was terminated due to the iteration limit, the error was manually checked and only recorded if the total error was less than 1 × 10 − 3 .An example of the numerically evaluated non-dimensionalized functions related to the stream function (f), velocity (f ′ ), shear stress (f ″ ), concentration (g), and concentration gradient (g ′ ) is shown in Fig. 1, where η is the non-dimensional distance from the plate surface.The MATLAB script further used the dimensionless velocity and concentration profiles to calculate the value of the non-dimensional distance from the plate surface at which u U∞ = 0.99 or c C0 = 0.99, which is the standard definition of the momentum or concentration boundary layer, using the function interp1.
The curvature of the non-dimensional concentration and velocity profile is large at the point defined as the top of the boundary layer.Therefore, two high order options for discretizing the concentration and velocity functions in the intersection algorithm interp1, 'pchip' and 'spline' methods, were investigated to ensure accurate calculation for the interpolation of the top edge of momentum boundary layer (η f ′ ) and the top edge of the CBL (η g ).The results of the interpolation from the pchip and spline methods were compared as a part of the interpolation algorithm to check for discrepancies.However, the discrepancy was insignificant (Δη ≅ 0.0001), so the average of the two methods was used as the final value of the CBL and momentum boundary layer.
'pchip' is a shape-preserving, piece-wise cubic interpolation method whereby the value at the query point is interpolated between neighboring grid points.'spline' interpolation uses not-a-knot end conditions and the interpolated value at a query point is based on a cubic interpolation of the values at neighboring grid points in each respective dimension.Data obtained with this method is herein referred to as 'Blasius Solution'.
Below is the step-by-step procedure used to generate data herein: (1) Define the Navier-Stokes and convection-diffusion equations as a set of first order differential equations in a MATLAB function file.(2) Define the solver settings and initial conditions for the ode set.
(3) Iteratively solve ode45 using the set of equations using the bisection method.(4) Extract the final eta vs non-dimensional parameter curves once the solutions reach convergence criteria.(5) Fit solutions generated in step 4 using the 'pchip' and 'spline' interpolation methods (parent function 'interp1′), extract the value of eta at 0.99 of the far-field or free-stream value, and then check for self-consistency between the two methods.

Computation of the flow and concentration fields
The steady-state 2D incompressible Navier-Stokes equations (Eqs.( 1) and ( 2)) were solved with the software COMSOL v6.1 for a 2D simulation domain using a first order (linear) discretization of the pressure and velocity variables.
where ρ is the density [kg/m 3 ], u is the velocity vector [m 2 /s], p is the pressure [Pa], I is the identity matrix, μ is the dynamic viscosity [Pa⋅s], and T is the transpose.In Fig. 1 the geometry and physical conditions are presented for the laminar flow over a flat plate scenario.Where (U ∞ ) is the free-stream velocity flowing over a flat plate.For the solid surface, no-slip boundary condition was applied.The flow was assumed to be fully developed and parallel to the plate upon interacting with the plate at the leading edge, where n is the normal vector.At the end of the plate, pressure was set to zero (p 0 = 0) together with a standard concentration outflow condition (n⋅D∇c = 0), where D is the molecular diffusivity and c is the concentration.A symmetry boundary was included in the computational domain as the plate also has flow occurring below the plate.However, for the purposes of this study, only solving one side of the plate was required.The flat plate was prescribed a constant concentration value of C 0 , simulating steady state diffusion of concentration out of the surface.The upper boundary of the computational domain was prescribed an open boundary condition [pI +K]n = 0, where p is the pressure, I is the identity matrix, and K is the stress tensor.The computational domain was scaled using the free-stream velocity, kinematic viscosity (ν), and critical Reynolds number to calculate the appropriate height (H) and length (L cr ) of the 2D domain.
The domain was created to be scaled automatically based on the Reynolds number, and CBL calculations were always made under laminar conditions, based on the critical Reynolds number for flat plate flow (Re < 5⋅10 5 ) (Lienhard and V., 2020).The solutions were nondimensionalized by normalizing the velocity and concentration components by their respective maximums.The maximum length of the plate was calculated using the standard assumption of the critical transitional Reynolds number for flat plate turbulent flow.The final domain size was selected to be half of this length to ensure laminar flow was always being evaluated, satisfying the Blasius assumptions for boundary layer flow.Even if the inlet velocity was changed, the domain geometry automatically scaled due to the parametric equations used to design the geometry based on the Reynolds number.Sampling of the velocity and concentration profiles perpendicular to the plate surface occurred at evaluation location, x E = 0.75⋅ , in the streamwise direction.The factor 0.75 was chosen to set the evaluation location sufficiently distant from the leading edge of the plate, but not too close to the end of the simulation domain.

Computing the concentration fields
The steady-state concentration field was calculated using the nonconservative form of the convection-diffusion equation using a firstorder discretization of Eqs. ( 3) and (4) in COMSOL Multiphysics v6.1. (3) where , and c i is the scalar concentration of species i [mol ⋅m − 3 ].It should be emphasized that the appearance of the velocity vector, u, appears in equation ( 3), demonstrating the clear interdependency of the evolution of the concentration field with respect to the velocity field.The steep gradients expected in the concentration field near the surface of the plate required a higher density of computational nodes.A second and more refined grid was created to minimize the number of elements needed to complete the coupled flow and concentration simulation.The velocity field was passed to the concentration solver as a 'set of initial values not solved' by the concentration solver.
The velocity and concentration profiles at x E were then constructed using Eq. ( 5) to non-dimensionalize the distance from the plate surface where y is the Cartesian y-coordinate above the plate, U ∞ is the inlet/ free stream velocity, and ν is the kinematic viscosity of the fluid (water at 310.15 K for this study).
The dimensionless velocity and concentration profiles were exported to MATLAB to determine the value of η f ′ →f ′ ( u U∞ ) = 0.99 and η g →g( c C0 ) = 0.99, the commonly understood limit at which the boundary layer ends (i.e., the top edge of the boundary layer or the boundary layer thickness).For both sets of data, the function interp1 was used as described in section 2.1 to determine the top edge of the boundary layer.The data obtained using this method is hereafter referred to as CBL-1.The same velocity solution from Section 2.2 was used for calculations with Johansson et al.'s method (hereafter referred as CBL-2), and the distance from the plate was treated in the same manner as described in Section 2.2.1.CBL-2 uses the Péclet number to estimate when the diffusion time scale dominates over the advection time scale.We set the Péclet number for this to 100, the same as in original published simulations.Subsequently, the velocity profile was converted into a local Péclet number along the perpendicular line probe at x E by Pe y = u⋅y D .We then calculated the intersection of Pe vs η at Pe = 100 using the intersections function in MATLAB.

Using the
Step-by-step procedure summary for Methods in Section 2.2.x: (1) Choose the flow solver settings for 2D incompressible Navier-Stokes fluid flow using a first order element discretization (CBL-1 & CBL-2).(2) Choose the solver settings for the transport of diluted species, used to solve the convection-diffusion equation using a first order element discretization (CBL-1 only).
(3) Solve the model and export the position (CBL-1 & CBL-2), velocity (CBL-1 & CBL-2), and concentration data (CBL-1 only) to a. csv file using the mesh to set the export resolution of data.(4) Import the data from step 3 into MATLAB and convert data into non-dimensional values.(5) Filter the data to ensure that no artifacts are present in the final data set.(6) Apply the interpolation function, 'interp1′ to each dataset, once using the 'pchip' option and once using 'spline' option.Compare the self-consistency of the fitted value of eta at 0.99 of the farfield or free-stream value (CBL-1 only).( 7) Apply the 'intersections' function, to the Pe transformed velocity dataset to obtain of the fitted value of eta at Pe = 100 (CBL-2 only).

Grid independence test
The grid independence study focused on resolving the flow and concentration fields for the most challenging case simulated in this study (Fig. 3).In particular the concentration transport requires a higher degree of spatial discretization to resolve the sharper gradient near the surface.Both the flow and concentration solutions were primarily assessed for grid independence by calculating the grid convergence index (GCI).Additionally, the validity of the CFD results were directly evaluated by comparing the calculated boundary layer values to the Blasius solution (see Fig. 6A).Moreover, GCI was applied to estimate the error arising from the spatial discretization of the domain.Expressed as a percentage, GCI is the difference between the result obtained from the most resolved mesh and the expected result from a mesh of infinite size (asymptotic solutions).To estimate the asymptotic solution, a series of meshes of increasing density must be simulated.The parameters of interest for the grid independence test were the area under the curve (AUC) for the non-dimensional velocity and concentration profiles.The AUC was used because both the velocity and concentration curves need to be precise enough to capture the point in dimensionless space where the local value of velocity or concentration is 99 % of the free stream/ far-field value.Both the dimensionless velocity and concentration area under the curve were calculated using Graphpad Prism to integrate ∫ 10 0 u U∞ dη, C far− field dη for a series of meshs in COMSOL with increasing spatial resolution.A domain of η from 0 to 10 for both the flow and concentration integrals was chosen as it is known a priori that both the non-dimensional velocity and non-dimensional concentration profiles reach the asymptotic value of 1 before an η value of 6.We restrict the domain for comparing of the AUC between η equal to zero and ten.For the GCI calculation, the parameter of interest was expressed as f (3, 2, 1) where 1 is the most refined mesh, 2 is the intermediate, and 3 is the most coarse.r is the refinement factor between mesh sizes, which in this study was constant:r = 4. Finally, ∊ is the relative error between the coarsest and intermediate meshes (Eq.( 6)), and p is the order of convergence (Eq.( 7)).These constitute the inputs for the GCI calculation for meshes with equal order refinement (i.e.r 1− 2 = r 2− 3 ) (Eq. ( 8)).
The main factor in concentration transport affecting the quality of the CBL is the diffusion coefficient.Thus, a range of representative Schmidt numbers was estimated, from 6.96⋅10 2 to 6.96⋅10 4 , on the basis of diffusivities of typical pharmaceutical entities (1⋅10 − 11m 2 s < D < 10⋅10 − 10m 2 s ) at physiological temperature.However, this work limited the range of the studied Schmidt number to ~1000, which approximated the representative low-molecular-weight drug ibuprofen.Entities with a diffusivity lower than 1⋅10 − 10m 2 s are most likely assemblies of monomeric species such as micelles or are macromolecules, such as peptides and proteins.The grid with refinement factor of 8 (~320,000 elements) was used to solve the fluid flow and the grid with refinement factor 16 (~1.3× 10 6 elements) was used to simulate the concentration transport.

Results
The highlighted results of the MATLAB interpolation used to determine the CBL are shown in Fig. 4A for CBL-1, and Fig. 4B for CBL-2.Johansson et al.'s original publication set the CBL as the distance away from the plate at which the Péclet number equals 100.In this setup, the y-component of the velocity is much smaller than the xcomponent.Thus, either the velocity magnitude or the x-component can be used for the Péclet number calculation with less than a ~0.1 % difference for the specific geometry studied herein.
The non-dimensional thickness of the CBL determined by CBL-1 and CBL-2 is directly compared in Fig. 5 using the same evaluation location on the plate for both methods.
Finally, the ratio of the momentum boundary layer and CBL was compared for the Blasius solution, CBL-1, and CBL-2 (Fig. 6).This ratio is critical for determining the exponent of the Schmidt number in a Sherwood number model.The Sherwood number can be used as dimensionless mass transfer model where Sh∝C⋅Re a Sc b , C is a constant unique to the geometry, Re is the Reynolds number, and a & b are constants unique to the system.For the flat plate, it is generally accepted that the average Sherwood number can be predicted via Sh = 0.644⋅Re 1/2 Sc 1/3 ; thus the ratio of the momentum and CBL thickness should follow a trend of Sc 1/3 (Sugano, 2010).

Discussion
CFD was applied to predict the CBL.The flow over a flat plate was used as a model geometry due to having an exact solution for the flow field, the Blasius solution, that could verify the CFD based approach.Two methods of determining the CBL, CBL-1 and CBL-2, were compared directly using the well-understood flat-plate geometry case.While this allowed for convenient comparison of the methodologies in this study, geometries of emerging dissolution-permeation apparatuses may not be simple enough for simple analytical solutions.These new devices would benefit from a CFD-based method with the ability to predict fluid flow and concentration flow in irregular geometric situations.
The ratio of the momentum and concentration boundary layer follows the cube root of the Schmidt number for convective mass transfer from a flat plat for Schmidt numbers between 1 and roughly 60, i.e.
where the exponent, n, in the relationship η f′ /η g = Sc n is equal to 1/3.The estimated lower limit to the range of values the Schmidt number takes for traditional small molecular weight pharmaceuticals is Sc ~ 700-1000 and is consistent with diffusion coefficients observed for single molecules (Sugano, 2010) The range of Schmidt number simulated in this study ranged from 1 to 1,000 and was chosen on the basis of the kinematic viscosity of water (0.   for CBL determination, it should be noted that the kinematic viscosity in the fasted state of the human gastrointestinal tract is slightly greater (0.7602 mm 2 /s) than that used here (Litou et al., 2016).
In general, applications of the Schmidt number in the human gastrointestinal tract for pharmaceuticals tend to have a larger value than what is studied here.However, most preclinical in vitro permeation experiments are initially performed in the fasted state or at least in a simulated, fasted-state intestinal fluid (Vinarov et al., 2023).It is commonplace to perform further experiments on formulations in fedstate media as the drug candidate selection progress progresses.While the fasted intestinal fluid has a viscosity similar to that of water, the human fed state fluid has higher viscosity and depends upon the content of the meal.The Schmidt number depends not only on the diffusivity of the molecule of interest in that specific fluid and the viscosity of the fluid, but also the density of the fluid.Increases in viscosity can increase the Schmidt number, but an increase in density lowers the Schmidt number.In the case of complex media, like the high-fat liquid meals, there are different entities within the media or that the gastrointestinal tract produces in response to that media which could modify the diffusivity of a molecule, by processes such as complexation, dimerization, or adsorption to or partitioning into micelle structures.This could dramatically increase the Schmidt number to an estimated range of 7000-100,000.From the computational perspective, as the Schmidt number increases, the gradient becomes significantly steeper across the boundary layer.This means that when simulating the mass transfer process through the boundary layer, finer spatial discretization is required.In this work we ensured that a large enough number of computational cells were located in the CBL region (approximately 3 cells per 10 µm in the final mesh) to capture the steep gradients near the solid surface.In situations where the Schmidt number is larger than what was studied here, this spatial resolution my not be sufficient and would require a study of the boundary layer gradient through a metric such as the grid convergence index.This presents a challenge for 3D geometries as the dense mesh required to resolve the boundary layer physics adds a lot of additional degrees of freedom, which results in a more computationally expensive calculation.Fig. 3 demonstrates the refinement of the computational grid and shows that the solutions to the dimensionless velocity (A) and concentration (B) profiles approach an asymptotic limit, as expected when refining the mesh.This figure demonstrates that the finest mesh adequately captures the gradient near the plate surface at the most refined meshes.The GCI calculations showed that the n = 8 grid was sufficient to resolve the velocity field, while the concentration field required a finer one (n = 16).This two-step process allowed the required interpolation of the velocity values on the new concentration grid without having to run a denser simulation for the velocity field.
When visually comparing the two approaches to determine the CBL, it is easy to see where the CBL is defined in CBL-1, at the point on the curvature of the concentration profile just before it heads towards the far-field value (Fig. 4A).Fig. 4B visualizes the CBL-2 method.It shows that as the Schmidt number increases, the spacing between the Péclet number curves decreases in a similar manner as in Fig. 4A, albeit at a much lower magnitude than the Blasius solution.It is clear that the CBL thickness estimated by CBL-2 is much smaller than what is generally accepted (Fig. 5).In fact, if the value of the intersection Péclet number is set to ~50,000-the value determined by trial and error curve fitting-the curves for CBL-2 and CBL-1 generally overlap while retaining slightly different slopes at low (Sc = 1-30) and high (Sc > 200) Schmidt numbers.This indicates a major change to the prescription of CBL-2.Johansson et al. do state in their discussion that the diffusion coefficient of drug molecules is ~10 − 3 m 2 /s.However, typical small molecular weight drug molecules have diffusivities on the order of ~10 − 10 m 2 /s (~10 − 6 cm 2 /s) and this discrepancy may contribute to their original conclusions.
The results of both CBL methods were compared against the Blasius analytical solution for a flat plate (Fig. 6).For CBL-1, computing the CBL via the ratio of the momentum boundary layer to the concentration boundary layer, agreed well with the Blasius solution, especially at low, potentially pharmaceutically irrelevant Schmidt numbers 1 < Sc < 100-200.At larger Schmidt numbers, in the range expected for most pharmaceutical applications, CBL-1 and the accepted Schmidt number correlation differed up to 6.6 % for Sc = 900.The MATLAB solver method used in the Blasius analytical solution was not optimized, as the goal of this study was to investigate the approach of using CFD methods, implemented in validated software packages like COMSOL, to form a basis for CBL predictions.However, we compared CBL-1 and the analytically determined Blasius solution and saw that they differed up to ~14 % at Sc = 900 (the highest Schmidt number simulated using CFD in this work).It is difficult to generalize given the multitude of dissolutionpermeation devices in use for preclinical evaluation of API and dosage forms, nevertheless, our analysis does show that the CBL-1 method estimates the CBL thickness to a reasonable degree of accuracy.The success of CBL-1 relies on the mesh grid in the CBL region being appropriately refined, so that the sharper gradient of concentration can be captured at high Schmidt numbers.Therefore, using CFD to solve for the CBL thickness could require two computational grids-one for velocity and one for concentration.This will avoid over-discretization of the velocity field and under-discretization of the concentration field.
In future applications of the CBL-1 methodology, analyzing the concentration boundary layers of the complex 3D geometries of real dissolution-permeation devices is the end goal.The same method to capture the edge of the boundary layer is directly translatable to the 3D space.In this work a line probe is used at one or potentially many positions on the length of the plate to capture the dimensional variables in dimensional space.For the 3D case, we simply create additional line probes on 2 axes to describe a plane instead of a line, sampling in the direction perpendicular to the absorption surface.We anticipate that we can then plot a boundary layer surface and then use an average value of the thickness as an input into traditional modeling strategies for permeation.Since it is possible that many dissolution-permeation devices operate in a regime where the fluid flow is partial or fully turbulent, and because turbulence is a highly time dependent phenomenon, we hypothesize that time-averaged analysis will be necessary to extract useful data from these more complex simulations.In early evaluations of the Johansson et al. method for 3D geometries, a work flow was developed to capture the time dependent evolution of the CBL and the work flow is being modified to be compatible with the CBL-1 method.Our hypothesis, which is currently being tested in the second phase of this work, is that steady (time-averaged) or time dependent solution can be used with this CBL-1 approach as long as the flow structures that dominate the convection process are adequately resolved.

Conclusions
The use of a CFD approach as presented here is useful for studies of mass transport performance.Estimating the CBL thickness, via the correlation Sc b = η f′ /η g , a maximum error of 14 % at high Schmidt numbers (ca.Sc = 900) was acquired, representing small pharmaceutical molecules.Although this was a simple flow case, this approach is promising and served as a proof of concept for characterization of more complicated and non-standardized device geometries used in preclinical evaluations of drug dissolution and permeation.

Fig. 1 .
Fig. 1.Example solutions of the Blasius boundary layer solutions for flow and concentration in MATLAB at arbitrarily chosen Schmidt numbers equal to 75 and 750.(A) Plot of the non-dimensional functions related to: the stream function (f), velocity (f'), and shear stress (f″) plotted against the non-dimensional distance from the plate surface (η).(B) Plot of the non-dimensional functions related to: concentration (g) and concentration gradient (g′) plotted against the non-dimensional distance from the plate surface (η).

Fig. 2 .
Fig. 2. Computational domain of the flat plate geometry and boundary conditions in COMSOL v6.1.The light blue dots on the upper and lower sides of the domain indicate points that split the domain in two.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Péclet number to calculate the concentration-boundary layer thickness according to Johansson et al Johansson et al. proposed that the CBL could be calculated strictly from the velocity profile by means of the Péclet number, Pe = L 2 /D L/u → Characteristic diffusion time Characteristic advection time .Where L is a characteristic length of the system [m], u is the velocity (fluid) [m/s] and D is the diffusivity [m 2 /s].Here, authors Johansson et al. hypothesize that when the time scale for advection is much shorter than the time scale for diffusion then advection transport dominates, specifically when the rate of transport caused by diffusion is 100 times lower than the rate of transport caused by advection.At the physical location in space where this occurs (i.e.η at Pe = 100) was defined as the outer edge of the concentration boundary layer."

Fig. 3 .
Fig. 3. Grid independence test for the COMSOL mesh using a Schmidt number of 929 used in (A) the flow field solver and (C) the concentration field solver.The grid convergence index was calculated based on the parameter area-under-the-curve (AUC) for the (B) velocity solution and (D) for the concentration solution.Note: The AUC was calculated with concentration or velocity on the ordinate, and η on the abscissa.Panels (E) and (F) represent an example solution at the leading edge and end of plate respectively for the velocity solution on calculated on the finest grid.Panels (G) and (H) represent an example solution at the leading edge and end of plate respectively for the concentration solution calculated on the finest grid.
6959 mm (Yukuyama et al., 2023)/s) at biological temperature (37 • C), the typical diffusivity of small molecular weight drugs, and what has been classically validated for this geometry.Keeping in mind future applications of CFD-based methods

Fig. 4 .Fig. 5 .
Fig. 4. CBL determined via (A) CBL-1 and (B) CBL-2.For CBL-2, η g was calculated from the intersection of the local Péclet number (Pe y = u⋅y D at x = x E ) at a y-distance from the surface of the plate where Pe = 100 (the dotted vertical line in panel B).The distance of the y-intersection was non-dimensionalized by η = y

Fig. 6 .
Fig. 6.Ratio of momentum and concentration boundary layers calculated using (A) the Blasius solution of the momentum and concentration boundary layers (green dots), and the corresponding solution from COMSOL (CBL-1, red dots).(B) Determination of the CBL via the Péclet number as proposed by Johansson et al. (CBL-2) is shown as blue diamonds.Points that deviate in the Blasius solution (green dots) were terminated by iterations (n = 100).The total error for these deviated simulations was less that 1 × 10 − 3 , but greater than the convergence criteria (1 × 10 − 6 ).The vertical dashed line indicates a Schmidt number of 60.For CBL-2, η f ′ was exactly same as for CBL-1, because CBL-2 does not calculate η f ′ .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Sinko: Writingreview & editing, Writingoriginal draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Louis Parker: Writingreview & editing, Visualization, Supervision, Methodology, Conceptualization.Lisa Prahl Wittberg: Writingreview & editing, Writingreview & editing, Visualization, Resources,