Uncertainty quantification for the squeeze flow of generalized Newtonian fluids

The calibration of rheological parameters in the modeling of complex flows of non-Newtonian fluids can be a daunting task. In this paper we demonstrate how the framework of Uncertainty Quantification (UQ) can be used to improve the predictive capabilities of rheological models in such flow scenarios. For this demonstration, we consider the squeeze flow of generalized Newtonian fluids. To systematically study uncertainties, we have developed a tailored squeeze flow setup, which we have used to perform experiments with glycerol and PVP solution. To mimic these experiments, we have developed a three-region truncated power law model, which can be evaluated semi-analytically. This fast-to-evaluate model enables us to consider uncertainty propagation and Bayesian inference using (Markov chain) Monte Carlo techniques. We demonstrate that with prior information obtained from dedicated experiments - most importantly rheological measurements - the truncated power law model can adequately predict the experimental results. We observe that when the squeeze flow experiments are incorporated in the analysis in the case of Bayesian inference, this leads to an update of the prior information on the rheological parameters, giving evidence of the need for recalibration in the considered complex flow scenario. In the process of Bayesian inference we also obtain information on quantities of interest that are not directly observable in the experimental data, such as the spatial distribution of the three flow regimes. In this way, besides improving the predictive capabilities of the model, the uncertainty quantification framework enhances the insight into complex flow scenarios.


Introduction
Non-Newtonian fluids are encountered in many industrial applications such as food processing, additive manufacturing and coating.In these applications, different manufacturing techniques are used, e.g., film blowing, injection moulding and extrusion.To understand and optimize these processes, models that accurately predict the involved complex flows, i.e., flows that consist of mixed, time-dependent and spatially varying deformation modes, are of the utmost importance [1].
The traditional modeling approach for complex flows assumes that model parameters are "deterministic", i.e., they are known exactly.Typically, these parameters are obtained by (non-)linear rheological measurements in well-defined flows, such as simple shear.Based on the measurements and external factors, such as experts' knowledge, intuition and experience, a constitutive model is selected, which is calibrated using the rheological data [2].The model is subsequently used to predict complex flow behavior.
An accurate prediction of fluid flows with complex rheological behavior requires sophisticated constitutive models [1].When calibrated well, such models can replicate the measured behavior in simple flow (as opposed to complex flows defined above) with a high accuracy.However, calibrating rheological parameters using rheological measurements is a non-unique procedure [3,4].In addition, since the deformation types and flow history can be significantly different in complex flow scenarios, using the calibrated model in such scenarios does not necessarily lead to accurate predictions.Moreover, if the flow is surrounded by uncertainties, e.g., from operating conditions and flow behavior, predictions using the traditional modeling approach may be complicated further.
Uncertainty quantification (UQ) provides an alternative modeling framework for making predictions in settings where it is hard to accurately determine model parameters.In this "probabilistic" approach, it is assumed that the model parameters are stochastic, acknowledging the limited knowledge we have about them.The parametric uncertainties are determined from measurements, intuition or experts' knowledge.In UQ, the main goal is to quantify the uncertainty of the model predictions, as opposed to obtaining the most accurate predictions possible.Whether or not the uncertainty is acceptable depends on the considered problem.
There are several important aspects to UQ, which include mapping out experimental uncertainties, sensitivity analysis, uncertainty propagation and inference.Mapping out the uncertainties indicates whether model parameters should be considered stochastic or deterministic.Sensitivity analysis enables us to determine which parameters have the largest impact on the model prediction [5].In uncertainty propagation, we propagate the parametric uncertainties through the physical model to ob-tain a probabilistic model prediction.In inference, we incorporate experimental data from the complex flow to investigate the uncertainty in the model parameters with the goal of improving the probabilistic prediction, by updating the model parameters.
The concept of uncertainty quantification has been widely used in, e.g., nuclear reactor models [6], weather models [7], biological models [8] and solid mechanics [9].An extensive review for Bayesian inference in physics is given by Von Toussaint [10].In the scope of non-Newtonian fluids mechanics and rheology, uncertainty quantification is relatively unexplored.In complex flows, uncertainty propagation has been used to evaluate the uncertainty of a quantity of interest.In a study by Pereira et al. [11] uncertainty propagation has been used, where the uncertainty in blood viscosity is propagated through the momentum balance to obtain quantities of interest such as the wall shear stress.Kim et al. [12] have also used uncertainty propagation to predict the resistance and velocity of a settling sphere in a Carbopol and PVP solution, where the rheological parameters were obtained using Bayesian inference.Kumar et al. [13] have used deep neural networks (DNN) to characterize a non-Newtonian fluid and have applied uncertainty propagation to investigate the dominant parameters that affect the simulation.In a study performed by Sen et al. [14] uncertainty propagation has been applied to Newtonian and non-Newtonian fluid flows on the microscale level.
Uncertainty quantification in the form of Bayesian inference has been mainly applied to rheological studies.In a study performed by Freund and Ewoldt [15], Bayesian inference has been used for UQ and model selection in linear rheology.Ran et al. [16] have applied Bayesian inference in the rheological characterization of a kaolinite clay suspension, which can be modeled using a micro-structural viscoelastic model.Bayesian inference has also been used in the prediction of linear viscoelastic models of branched polymers in a study performed by Shanbhag [17].UQ has also been used to estimate rheological parameters in related fields such as geophysical research [18,19].
Sensitivity analysis in a complex flow setting has been performed by Freund et al. [20], who studied field sensitivities of flow predictions to rheological parameters for generalized Newtonian and thixotropic fluids.This work was later extended to viscoelastic fluids at a low Deborah number by Kim [21].
To summarize, uncertainty quantification has been applied in the field of non-Newtonian fluid mechanics in the form of sensitivity analysis, uncertainty propagation and Bayesian inference.However, inference has been mainly limited to simple flows as found in rheology, for which fast-to-evaluate models are available.To the best of our knowledge, no research has yet been done on applying Bayesian inference in a complex flow of non-Newtonian fluids.
In this work, we use the methodology of uncertainty quantification in the context of a complex flow case with a generalized Newtonian fluid.The study of this setting introduces two requirements.First, we need full control of our experiments, allowing for detailed assessment of the parametric uncertainties and model calibration.Therefore, we consider a squeeze flow for which we have developed a tailored experi-mental setup.A squeeze flow consists of extensional and shear deformations, making it a complex flow which is still computationally tractable.Second, the computation time for our physical model should be kept to a minimum, because it has to be evaluated many times to allow for uncertainty propagation and inference.To this end, we develop a semi-analytical model that allows for describing Newtonian and shear thinning flow behavior.
This manuscript is outlined as follows.To set uncertainty quantification in the scope of the current work, in Section 2 we first introduce the key concepts and terminology for UQ.In Section 3 we then introduce the squeeze flow setup and the corresponding experiments.Section 4 describes the Newtonian and shear thinning model to mimic the experiment.In Section 5 we describe the quantification of the parametric uncertainties.We then discuss the comparison between the squeeze flow model and experiments using uncertainty propagation and Bayesian inference in Section 6.Finally, we present the conclusions and recommendations in Section 7.

Uncertainty quantification
In this section we outline the fundamental ideas of uncertainty quantification, introducing the terminology and notation used throughout this work.We refer the reader to, e.g., Oden et al. [22] for a more detailed exposition.

Uncertainty modeling
A fundamental goal of science is to make predictions of events in physical reality, or, more precisely, about quantities of interest related to such events.In modern science and engineering, the systems for which predictions are required are becoming more complex.The complexity of the systems gives rise to uncertainties in experimental observations used for model validation and for the calibration of model parameters, as well as uncertainties in the models themselves.
Uncertainty quantification (UQ) aims to make predictions of events in physical reality using models and taking into account the uncertainties that are inherent to the problem.A particular challenge in many engineering problems is that calibration of models using relatively simple experiments insufficiently reduces the uncertainties associated with the predictions for the real system of interest.
In UQ, quantities that are uncertain are described by random variables, which are characterized by probability distributions.To formalize this concept, we consider a distribution g of a vector of quantities that describe the "true" event, and experimental observations y = {y 1 , ..., y n } of this event, where n represents the number of observations.Denoting the realization of reality corresponding to the i-th observation y i as g i allows us to express the observation error ε i (in an additive way) as We note that the truth, g i , is unknown in practice and that in general a model should be assumed for the observational error.Such a model is generally referred to as a noise model.
Similar to the definition of the observation error, we can also define a model error (or model bias) γ i (θ) as where d i is the model prediction corresponding to the parameter set θ.We define the parameter domain as Θ = Θ 1 ×Θ 2 ×...×Θ p , where θ ∈ Θ and p is the number of parameters.
Combining the observational error (1) and the model bias (2) enables the elimination of the unknown truth from the error definition as This expression conveys that the difference between the model predictions and the observations is a combination of two error contributions, viz. the model errors and the observational errors.
A fundamental aspect of uncertainty quantification is to postulate models of the error contributions, as these models play an essential role in the uncertainty in the problem.In practice, it is frequently decided to combine the two error contributions in equation ( 3) in a single noise model, encompassing both error sources.
In this work we consider both experiments and models, which enables us to systematically study both error contributions.In particular, we study the relation between the simple calibration experiments, i.e., the rheological measurements, and the more complex application, i.e., a squeeze flow.Before considering the problem of interest in this work, we first discuss the methods used in the framework of UQ.

Uncertainty propagation and inference
The two main approaches of UQ are uncertainty propagation and inference.In uncertainty propagation, we quantify the uncertainty in model parameters, i.e., parametric uncertainties, θ, and propagate these through a physical model, d, to obtain the uncertainty in the quantity of interest.In inference problems, the model parameters are updated based on observations.We evaluate the posterior probability density function, π(θ|d), of the parameters based on both prior information of these parameters, π 0 (θ), and on observations y, using Bayes' rule In this expression, the normalization constant, π Υ (y), is referred to as the evidence.The posterior quantifies the probability of obtaining parameter values, θ, given the experimental data y.The observations are encoded in the likelihood function, L(θ|y) ≡ π(y|θ).
By combining prior knowledge with experimental data, an improved model prediction is obtained compared to uncertainty propagation.The experimental data influences the posterior density through the likelihood term L(θ|y), which quantifies the likelihood of observing the data y given the parameters θ.In the likelihood function, the measured data y is fixed and the parameters θ are varied over the admissible domain.Because it is assumed that the data remains fixed and the parameter values can vary, the likelihood is not a probability distribution.The model errors and observational errors, as discussed in relation to equation (3), are incorporated in the likelihood function as random variables.A common way to do this is to assume the combined errors to be independent and identically distributed, following a normal distribution.We then obtain where Σ = {Σ 1 , ..., Σ n } is the variance of the noise.It is noted that alternative (e.g., multiplicative) noise models can be used [23].Furthermore, we note that, to avoid arithmetic underflow problems, in our implementation we consider the log-likelihood function l(θ|y) = ln (L(θ|y)) , instead of the likelihood function itself.The prior density in equation ( 4) is based on the acquired knowledge prior to obtaining the experimental data.This term is based on experts' knowledge, intuition, or previous experiments.In uncertainty propagation, we use the prior as the parametric uncertainty which is propagated through the model.If no relevant information is known about a parameter, one typically uses an uninformative prior, for which often a uniform density is used [24].If sufficient and convincing information is obtained for a parameter's value, the prior is informative with a relatively small uncertainty.If the uncertainty in the prior is relatively small in comparison to the experimental noise, the posterior will strongly be influenced by the prior, and vice versa.
On account of the in general high dimensionality of the parameter domain, the evidence can be computationally demanding to evaluate since it involves an integral over all parameter values.There are several ways to evaluate it, where the method depends on the number of parameters in the systems.In special cases, the term can be evaluated analytically.For lowdimensional problems, i.e., p ⪅ 4, one typically uses quadrature techniques [5].In a moderate number of dimensions, (adaptive) sparse grids are used [25].For high-dimensional problems, Monte Carlo methods, such as Markov chain Monte Carlo (MCMC) algorithms [26], are commonly used.In these samplers, only the probability ratio of two subsequent steps in the chain is needed, making the evaluation of the evidence unnecessary.

Squeeze flow setup and experiment
In this section we specify all relevant aspects of the tailored experimental setup and experiment.

Fluids
The squeeze flow experiments are performed using two types of fluids: glycerol and a Polyvinylpyrrolidone (PVP) solution.The former is in liquid form obtained from VWR Chemicals with purity ≥99.5%.It has been exposed to air for two days to be saturated with water absorbed from the air.The PVP is obtained in powdered form (average molecular weight of 360 kg/mol) from Sigma Aldrich, which is dissolved in purified water.A 23 wt.%PVP aqueous solution is obtained by stirring the mixture for three days.

Instrumentation
To have control over all aspects of the experiment, we developed a tailored setup (Figure 1) to perform the squeeze flow measurements.The setup consists of an aluminium frame and a 3D-printed parallel-beam construction.Furthermore, two Polymethylmethacrylaat (PMMA) plates are used, in between which the fluid is compressed.One plate is attached to the 3D-printed construction and the other one to the aluminium frame.The positions of the setup before and during an experiment are shown in Figure 1a and Figure 1b, respectively.The red beams act as leaf springs, ensuring the parallel motion of the PMMA plates.To capture the fluid flow behavior, a Nikon Coolpix W300 camera is used 1 .The 3D printed parts are made from Polylactic Acid (PLA) and have been developed with the Ultimaker S5, which is a 3D printer using Fused deposition modeling (FDM).Each part is designed in the software SIEMENS NX and transmitted to the software Ultimaker Cura to enable printing2 .

Procedure
In our experiments, we extract glycerol using a syringe.The PVP solution is extracted using a small spoon because the fluid has a too high viscosity to be extracted by a syringe.The bottom PMMA plate is removed from the setup to deposit the material on it.When reassembling the setup, a release mechanism is placed between the top and bottom plate to ensure a controlled height of the fluid layer.We ensure that the fluid touches both plates before the start of the experiment to properly determine the initial height.The camera is remotely focused on the sample using the SnapBridge app.When the sample is in focus, the camera can start recording and the release mechanism is triggered.The motion of the fluid is captured for a duration of five minutes.After measuring, we use water to clean the plates and dry them before positioning a new sample.Each experimental case (see Table 1) is performed ten times.
We tested several cases for both materials, where we varied the applied force and fluid volume.The cases are listed in Table 1a and Table 1b for glycerol and PVP solution, respectively.The applied force is adjusted by placing additional weights on the top plate.The amount of additional weight is denoted by Table 1: Squeeze flow experimental cases: (a) using glycerol, where V a < V b (b) using PVP solution, where V c < V d .See Section 5 for details.

Case
F add [kg] Figure 2: Image processing for case VI at t = 9.14 s. (left) Original image retrieved from the camera.(right) Binary image created by subtracting the current frame from the initial frame and using Otsu's method [27] to determine the threshold value.Using the circle Hough Transform [28], the initial fluid front (blue) and the current fluid front (red) are determined, assuming a circular shape.Subsequently, the radius of the fluid layer is evaluated.

Processing
We analyze the experimental data by extracting the radius of the fluid layer using image processing, as illustrated in Figure 2. The frames are retrieved with a frame rate of 60 Hz.Due to the delay between the recording of the camera and the start of the measurement, the initial frame is determined by comparing the frames up until the fluid is set in motion.

Squeeze flow model
In this section we present the considered squeeze flow model.A schematic of the model is provided in Figure 3.The primary quantity of interest for our model is the radius of the fluid layer over time, R(t).To model the evolution of the radius, at any given moment in time, we assume the problem to be isothermal.Furthermore, we assume negligible inertial forces and the fluid to be incompressible.As a starting point, we consider the mass and momentum balance given by where v and p are the velocity and pressure field, respectively, both defined over the volumetric fluid domain Ω(t).The extra stress τ is related to these fields through a constitutive rheological model, which will be discussed in detail below.Note that the domain Ω(t) is time-dependent, being directly related to the radius R(t).In our model, the motion of the fluid front is governed directly by the velocity field at the fluid-air interface, Γ 3 (t).We assume the solution to be axisymmetric, meaning that we assume independence of the circumferential coordinate.Assuming the fluid layer to be thin in comparison to the radius, i.e., H ≪ R, we use lubrication theory [29] to dimensionally reduce the mass and momentum balance (7) where v r is the r-component of the velocity vector v, v z the zcomponent of v and τ rz is the rz-component of the stress tensor τ.All other components of the velocity vector and extra stress tensor are neglected.
To solve these conservation equations on the evolving domain, we develop a model using the simulation domain given in Figure 3b, where the boundaries are denoted by Γ 2 , Γ 3 and Γ 4 and the axisymmetric axis denoted by Γ 1 .At Γ 2 we have a noslip condition between the fluid and the plates, i.e., v r (r, h) = 0. Furthermore, the vertical velocity v z equals half the velocity of the upper plate Ḣ, i.e., v z (r, h) = 1 2 Ḣ = ḣ.An additional condition is set at Γ 2 stating that the pressure integrated over the top plate equals the applied force, i.e., R 0 p(r, t)2πr dr = F.Note that the atmospheric pressure is omitted, which is justifiable due to the fluid incompressibilty assumption.At the axisymmetry axis, Γ 1 , the radial velocity equals zero, i.e., v r (0, z) = 0.At the horizontal symmetry plane, Γ 4 , the shear stress is equal to zero, i.e., τ rz (z, r, t) = 0.At the free boundary Γ 3 we have a traction boundary condition where the normal traction equals the Laplace pressure, i.e., p(R(t), z) = ∆p, with The Carreau model and the truncated power law model in simple shear (with shear rate γ) describe the Newtonian regimes (i.e., at the zero shear rate viscosity and at the high shear rate viscosity) and the power law regime similarly.The transition between these regimes is smooth for a Carreau model and sharp for a truncated power law model.Note that the axes are assumed logarithmic in this schematic.
where γ is the surface tension, ∇ s is the surface gradient operator defined by ∇ s = (I − n ⊗ n) • ∇, where I is the unit tensor and n is the normal to the interface, and κ the curvature.The curvature around the z-axis is assumed to be significantly smaller than the curvature between the plates (κ −1 ).The former curvature is therefore neglected.
We herein consider two constitutive models to be used in combination with the squeeze flow model: a Newtonian model and a generalized Newtonian model.The former describes fluids where the viscosity is independent of the shear rate, while the latter is able to mimic fluids with combined Newtonian and shear thinning behavior.A commonly used model for shear thinning is the Carreau model [30], which combines Newtonian behavior at high and low shear rates, and introduces a shear thinning behavior for intermediate shear rates.However, we use a truncated power law (TPL) model instead of a Carreau model because the squeeze flow model can then be solved semi-analytically, while still describing similar flow behavior, thereby saving computational time.The truncated power law model is illustrated in Figure 4 to show the main difference with the Carreau model.
In the remainder of this section the squeeze flow model will be elaborated for the two considered constitutive models.It is noted that formally the Newtonian model is a special case of the truncated power law model.For the clarity of the derivation, we however opt to first consider the Newtonian model as a separate case.

Newtonian model
To set the scene for the derivation of the truncated power law model in Section 4.2, we briefly review the fundamental steps in the derivation for the Newtonian case, including the Laplace pressure boundary condition.We note that in the absence of the Laplace pressure, the result is well established in the literature, see, e.g., [31].
The first step in the derivation of the expression for the radius (and height) over time is to use the momentum balance (8b) to acquire the through-the-thickness velocity profile as a function of the pressure gradient ∂p/∂r and the rheological parameters.For the Newtonian case, the viscosity is a constant denoted by η N , thus τ rz simplifies to where D rz is the rz-component of the rate-of-deformation tensor D. By substitution of this expression in the momentum balance (8b), we obtain the radial velocity where use has been made of the boundary conditions v r (r, h) = 0 and τ rz (r, 0) = 0.
To obtain the pressure field over time, the obtained velocity profile is substituted in the thickness-integrated mass balance (8a) as where use is made of the boundary conditions v z (r, 0) = 0 and v z (r, h) = ḣ.Evaluation of the integrated mass balance then yields where we have used ∂p ∂r (r = 0) = 0 and p(r = R) = ∆p.Note that the pressure field ( 13) is a function of the speed of the top plate.Since our model is force-driven, the pressure field must be expressed in terms of the applied force instead of the plate's speed.To this end, we use the additional condition from where we obtain the height and radius of the fluid layer over time using the forward Euler method [32] and the initial condition h(t = 0) = h 0 .
Remark 1 (Zero Laplace pressure special case).The squeeze flow model can be solved analytically in the special case where we assume the Laplace pressure to be equal to zero [31].In that case one obtains from (14) that the height and radius evolve as

Truncated power law model
The truncated power law model describes Newtonian and shear thinning flow behavior, as illustrated in Figure 5a.Note that during our experiments we do not reach the third region.However, to improve numerical stability we also consider this region.In other research the truncated power law model has been used in e.g.[33,34].However, Lee et al. [34] have modeled two regions instead of three and Lavrov [33] has modeled three regions with the pressure gradient as pre-knowledge.Therefore, to the best of our knowledge, the three-region case with arbitrary pressure gradient, which we consider here, has not been reported.
The division of the regions in simple shear (note that the axes are assumed logarithmic in this schematic), (b) The division of the regions in a squeeze flow experiment at a fixed moment in time.The heights of the interfaces between the regions are denoted by w 1 and w 2 , which are both a function of r.These interfaces are shown in dashed purple and dark red, respectively.
Due to the lubrication assumption, the squeeze flow is rheometric, i.e., it can be considered a point-wise simple shear flow with absolute shear rate |γ| = ∂v r ∂z (Appendix A).The corresponding mathematical model for the viscosity then reads where η 0 , λ cr , n and η ∞ are the viscosity at the zero shear rate plateau, the inverse of the shear rate corresponding to the transition between region 1 and 2, the power index and the viscosity at infinite shear rate.Note that for notational convenience we have rewritten the power law constitutive model from its more common representation η = K |γ| n−1 , where K is the flow consistency index, by introducing K = η 0 λ n−1 cr .During a squeeze flow experiment, the maximum shear rate decreases from the high shear rate plateau to the low shear rate plateau over time.Furthermore, at a fixed moment in time, the shear rate near the top plate is higher than the shear rate in the middle of the fluid layer, as visualized in Figure 5b.Therefore, it is essential that the three possible regions of the flow behavior are incorporated in the squeeze flow model.
As for the Newtonian case, the first step in our derivation of the radius over time is to solve the momentum balance to attain the velocity profile as a function of the pressure gradient and rheological parameters.The fundamental complication in the case of the truncated power law is that the profile is subdivided in three regions, parameterized by the heights w 1 and w 2 , where the interface positions themselves depend on the pressure drop and rheological parameters.In our derivation, we commence with solving the momentum balance in each of the three regions, after which we use the continuity of traction and velocity across the interfaces to obtain the velocity profile with assumed values for the interface heights.We then determine the analytical profile by using the continuity of the viscosity to derive expressions for the interface positions.
To start, the momentum balance is solved by implementing the relevant constitutive models per region.The extra stress is then defined as Using the boundary condition τ rz (z = 0) = 0 and traction continuity at the interfaces z = w 1 and z = w 2 , we obtain the stress from the momentum balance as To obtain the velocity profile through-the-thickness, we match the expressions provided in equation ( 17) and (18).Making use of the boundary condition v r (z = h) = 0 and the continuity condition of v r at the interfaces gives v r,2 = ∂p ∂r The velocity profile is visualized in Figure 6 for numerous pressure gradients, yielding flows with shear rates in the Newtonian plateau and power law regimes.The velocity profile flattens as the gradient of the pressure increases, which is due to the shear thinning behavior of the fluid.Next, we determine the total flux by the addition of the fluxes per region as where the regional fluxes are determined by integrating the expressions for v r over the specified height per region, given by , where ∂p ∂r λcr is the pressure gradient that yields a maximum shear rate of 1/λ cr .The velocity is normalized by the Newtonian flux for the corresponding ∂p ∂r values.The black triangles and circles denote the shear thinning and Newtonian behavior, respectively.As the pressure gradient increases, the shear thinning behavior increase throughout the velocity profile.
We incorporate the flux Q in the conservation of mass as which is a non-linear differential equation in time, where the non-linearity stems from the dependence of the interface positions, w 1 (r) and w 2 (r), on the solution.Using the viscosity continuity condition, this dependence can be expressed as To integrate (23) in time, we have developed a semi-analytical non-linear time-integrator based on fixed point iterations.The typical runtime of our model is t sim = 2 s for a simulation time of t = 350 s.We refer the reader to Appendix B for details about the solver.
Remark 2 (Limiting cases of the truncated power law).When | ∂p ∂r | ≤ η 0 hλ cr it follows that w 1 = w 2 = h and we get η = η 0 and the solutions (13) and (14).When | ∂p ∂r | → ∞ it follows that w 1 = w 2 → 0 and we get η = η ∞ and the solutions (13) and ( 14), but with η ∞ instead of η 0 .Note that in these cases, the solution can be expressed analytically, but that this is not possible in the general case.
Remark 3. Due to physical considerations, the rheological model parameters η N , η 0 , η ∞ and λ cr and squeeze flow parameters V and R 0 must be non-negative.To avoid these parameters from becoming negative in the sampling method introduced in Section 5, a log-transformation is applied.This is done by introducing the transformation θ = exp( θ), where the model parameters are denoted by θ, and assigning a distribution to the transformed parameters θ.

Characterization of parametric uncertainties
The input parameters used to determine the quantity of interest can either be probabilistic or deterministic.Parameters can be probabilistic due to a lack-of-knowledge of a deterministic value or due to physical randomness in the system.
As defined in Section 4, the input parameters of the squeeze flow model are the fluid volume V, the applied force F, the initial radius of the fluid R 0 and the rheological parameters, i.e., η for the Newtonian model and η 0 , η ∞ , n and λ cr for the truncated power law.The squeeze flow experiment introduces additional uncertainties.These include the refraction of the camera lens and PMMA plates, the pixel-to-length ratio of the squeeze flow images and the uncertainties associated with the circular fit of the fluid front.We assume that these additional parameters are deterministic, because our experience suggests that the uncertainty of the model prediction is dominated by the stochastic parameters mentioned above (i.e., the squeeze flow model parameters).
The parametric uncertainties for R 0 , V and F are determined by performing dedicated experiments for each of the parameters.Each experiment is repeated ten times, from which a parametric distribution is defined.Because it is difficult to perform independent measurements for each of the truncated power law parameters, these constitutive model parameters are determined by Bayesian inference.

Sampling method
In our Bayesian inference problems, we use a Markov chain Monte Carlo (MCMC) sampling method, which draws a number of realizations from the posterior distribution by evaluating the prior and likelihood for any particular parameter value.As the sample size (i.e., the number of realizations) is increased, the sampled distributions can be expected to move closer to the posterior distribution.In the considered Markov chain Monte Carlo method, the next step in the chain of realizations is only dependent on the latest parameter value.Multiple chains (or walkers) are used to robustly explore the parameter space, circumventing the possibility of only sampling a local region of high probability.In this work, the walkers are initialized by sampling from the prior distribution.The number of steps it takes for a walker to move from the initial state to the typical set (the region with high posterior mass probability) is known as the burn-in period.
We herein use a specific MCMC algorithm: the affine invariant MCMC ensembler sample [35], which is implemented in Python through emcee [36].The pivotal idea behind this sampling algorithm is to use multiple chains (referred to as an ensemble) to explore the posterior distribution, and to use information from another chain in the proposal step.By doing so, the sampler is effectively capable of rescaling the parameter domain (referred to as affine invariance), thereby ameliorating sampler performance bottlenecks associated with anisotropic distributions.
We employ the affine invariant MCMC algorithm with its default stretch move, which is visualized in Figure 7.The proposal stretch move is accepted when r < q, where r is sampled from a uniform distribution between zero and one and q is given by where p is the dimension of the parameter space, θ k is the current step of walker k and θk its proposal step.The proposal is defined by θk where θ j is the current step of a randomly selected walker j k and Z is a random variable sampled from the probability density We herein use the default value of a = 2.
To investigate whether the chains have sufficiently converged, we apply an autocorrelation analysis by calculating the estimated integrated autocorrelation time τ θ for each parameter θ in the concatenated chain C = {θ n } N n=1 .The term "estimated" comes from the fact that we are looking at finite chains with concatenated length N. Subsequent steps in the chain can be dependent on one another, and therefore we need to define the number of independent samples within a chain.The estimated autocorrelation time for a single model parameters is defined as for some M ≪ N, where ρθ is the estimated normalized autocorrelation function of the stochastic process that generated the chain.We use M instead of N to lower the level of noiseto-signal ratio [37].The estimated normalized autocorrelation function is defined as where with Based on the maximum autocorrelation time, τ = max θ τ θ , for each chain we remove the first N burn = 2•τ samples for burnin, and we thin the chains by selecting every N thin -th sample in each chain, with N thin = 1 2 • τ.This results in a final sample size of We note that, preferably, the integrated autocorrelation time is as small as possible to save more of the generated samples instead of discarding them.

Rheological characterization
To characterize the fluids, the simple shear response is measured by performing steady rate sweep tests on a TA Instruments ARES rotational rheometer using a 25 mm diameter cone-plate geometry.To counteract dehydration of the PVP solution, a Plexiglas ring is attached to the rheometer base around the sample.By injecting silicone oil in between the sample and the Plexiglas ring, the water inside the PVP solution cannot evaporate.We have executed ten measurements for both glycerol and PVP solution.Based on the measurement data and a constitutive model, the likelihood is defined, in accordance with equation (5).The average of the experiments is defined as where µ y denotes a vector of average values per shear rate and n denotes the number of experiments (ten in this case).The Likelihood equation ( 5) can subsequently be written as where Σ y is the covariance matrix corresponding to the measurements.

Glycerol
Based on experience, we assume the constitutive model for glycerol to be Newtonian with a viscosity η N .We define a prior distribution for the viscosity of glycerol.Based on a study performed by Segur [38], the viscosity of glycerol equals 1.4 Pa•s.Because the viscosity cannot be negative (see Remark 3), we create a prior normal distribution for the log-transform ηN with µ=1.4 Pa•s and σ=0.3 Pa•s, to make the distribution weakly informative.To determine the mean and standard deviation of the lognormal distribution we use To sample from the posterior, we define the number of walkers/chains, number of samples per walker, the burn-in period and the thinning.We have used four walkers and 10,000 samples per walker.The burn-in period is defined as 2 × τ f and the thinning as τ f /2, where the autocorrelation time τ f for ηN equals 30.The effective sample size N final then equals 3,059 samples.We evaluated N = 4 × 10, 000 = 40, 000 samples out of which 3,059 have been used to make a model prediction, indicating that a substantial number of samples has been lost due to the burn-in period and (especially) thinning.
We construct a trace plot to investigate whether sufficient samples have been used to define the posterior distribution.If the chains fluctuate around a value for ηN , we can assume that sufficient samples have been generated for the chain to be stationary.The trace plot for η is shown in Figure 8.We observe that all the walkers fluctuate around a horizontal line and therefore a certain value for ηN .
The posterior predictive distribution of the rheological model is shown in Figure 9, together with the experimental data.Note that, as discussed in Section 2, experimental observations include both parametric uncertainties and observation errors.Including the latter in the model predictions yields the posterior predictive distribution (PPD) [9].We determine the PPD by sampling from the posterior (reusing samples from the MCMC), determining the model response and adding noise according to the noise model (34).Alternatively, the uncertainty in the model results could have been visualized through the posterior distribution, which yields probabilistic model realizations that only include parametric uncertainties.
Because the model is equal to the constant value ηN , the viscosity is independent of the shear rate and equal to the posterior distribution ηN .The posterior distribution for η N is visualized in Figure 10 together with the posterior distribution ηN .The posterior distribution for η N is obtained by calculating the exponential function of ηN (η N = exp ( ηN )) for each sample and creating a distribution from these samples.We observe that both distributions have a similar shape, but have a different mean and standard deviation.Because the standard deviation of η N is relatively small and not close to zero, the right-skewness characteristic for a log-transform distribution is not clearly visible.To evaluate the level of uncertainty in η N , the coefficient of variation, CV = (σ/µ) × 100, is determined.The mean µ, standard deviation σ and coefficient of variation CV are provided in Table 2.

PVP solution
The truncated power law model parameters η 0 , η ∞ and λ cr , which can physically not be negative, are approached in a sim- ilar way as the viscosity η N in the Newtonian model.The normal distributions of the log-transforms of these parameters are denoted by η0 , η∞ and λcr .We assign a uniform distribution between zero and one to n, because the material shows shear thinning behavior and n < 0 would lead to a non-physical maximum in the shear stress.
We have used eight walkers and 20,000 samples per walker.The burn-in period and thinning are defined similarly as for the inference applied on the Newtonian model for glycerol.The autocorrelation times for the rheological model parameters are provided in Table 3.The total number of samples N, and thus the number of physical model evaluations, equals 160,000.The final sample size is equal to 12,239, which is based on the highest autocorrelation time of the model parameters.
The constitutive model prediction and experimental results are shown in Figure 11.The distribution of and correlation between the model parameters η0 , η∞ , λcr and n are shown in Fig- ure 12 [39].The prior distributions are denoted by the solid red line and the histogram denotes the posterior distribution for the model parameters.The dotted black lines represent the 95% "credibility interval", i.e., the interval containing 95% of the area under the posterior distribution.The samples and distribution of the correlation between every two parameters are provided in the two-dimensional posterior marginals.We find that n and η0 show a strong dependence on the likelihood because the prior probability is low and flat as compared to the posterior, while η∞ is mainly determined by the prior information.The influence of the likelihood and prior are both visible in the posterior distribution of λcr .The correlation plots indicate that the strongest correlation exists between n and η0 .The statistical moments are provided in Table 4.

Independent parameter experiments
Because of the different cases used for the squeeze flow experiments, several distributions are obtained for the fluid vol- Experimental data (µ exp ± 2σ exp ) Figure 11: Posterior predictive distribution of η TPL (γ) using the truncated power law constitutive model, where the subscript 'm' corresponds to the model and the subscript 'exp' to the experimental data.Because of the logarithmic scale, the error bars are asymmetric per shear rate and appear to increase as the shear rate increases.ume, applied force and initial radius.Each of the parameters is given a (log-)normal distribution where the mean and standard deviation are based on ten measurements.
• Fluid volume.The uncertainty in fluid volume is determined differently for glycerol and PVP solution.For glycerol, the fluid volume is measured separately from the squeeze flow experiment.We pipette ten samples and weigh them to find the distribution in fluid volume, using the density of glycerol as obtained from the supplier.For PVP solution, each sample is weighted just before performing the squeeze flow experiment.The density of PVP solution is experimentally determined by weighing a sample of known volume.The fluid volume distribution applicable to several squeeze flow experimental cases are given in Table 5.
• Applied force.We measure the applied force using a spring suspension attached to the moving parallel plates structure, which is denoted by the black part in Figure 1.
We perform force-displacement measurements for several combinations of force and displacement to obtain the spring constant of the construction.Due to the minimal vertical displacement, the force is assumed to be constant in this range.The mean and standard deviation are determined from the linear fit onto the measured data.The applied force distributions are provided in Table 6.
• Initial radius.The initial radius is obtained from the image processing step by averaging the values of initial radius obtained from a single set of the squeeze flow measurements.Each squeeze flow case corresponds to a different distribution of the initial radius.The initial radius distributions are provided in Table 7.
• Laplace pressure.To account for the surface tension in the squeeze flow model, we implement the Laplace pressure (equation ( 9)), which depends on the surface tension γ  8.We have quantified the surface tension for PVP solution using literature [40].
Quantifying the curvature at the fluid front is more challenging, as it is not directly observable in the squeeze flow experiment.The curvature is implemented with a dependence on the semi-height of the flow domain h as where α is a factor valued between 0 and 1.For α = 0, the interface is straight, and thus the Laplace pressure is excluded.For α = 1, we have reached the maximum curvature, i.e., κ = 1/h.

Comparison of the models and experiments
In this section we present the comparison between the experimental and numerical results for the squeeze flow of glycerol and PVP solution.The comparison is realized by uncertainty propagation, using the forward Monte Carlo sampling method, and Bayesian inference, using MCMC.
We used Bayesian inference on the constitutive model to obtain a distribution for the rheological parameters (see Section 5).Based on the simple shear measurements and constitutive model prediction, we expect that the squeeze flow behavior of glycerol can be described by the Newtonian squeeze flow model.For PVP, a prediction for the initial shear rate using the Newtonian model and zero-shear viscosity, yields maximum shear rates in the order of 10 2 1/s.Since this is in the shearthinning regime (Figure 11), we use the truncated power law model to predict the squeeze flow behavior of PVP solution.
For both of the fluids, we start by using the method of uncertainty propagation to compare the numerical model to the experimental data, keeping the parameter distributions fixed.Thereafter, we use Bayesian inference to update the parameters, and thus the model response, using the experimental data from the squeeze flow.Finally, we will compare both methods for both fluids.

Newtonian fluid
The results of the experimental cases I to IV are visualized in Figure 13.From case I to III we have an increase in applied force (see Table 1a and Table 6).In Figure 13 we observe that an increase in applied force leads to a faster growth of the radius.The applied force in case IV is similar to the applied force in case II and the volume is increased (see Table 1a and Table 5).We observe that an increase in fluid volume leads to a slight increase in initial radius.Analyzing the experimental data, we see that the uncertainty increases as time proceeds.

Uncertainty propagation
To use the method of uncertainty propagation, we assign a distribution to the model input parameters.In Section 5, we have discussed experiments to obtain the uncertainty in the input parameters F, V, R 0 and γ and have obtained the distributions for these parameters.The curvature α has been defined in two limiting cases: α = 0 and α = 1.Using Bayesian inference on the simple shear measurements, we have retrieved a distribution for the rheological parameters.Every distribution is to be propagated through the squeeze flow model to obtain the distribution of the radius over time.We have used 10,000 samples from each distribution of the model parameters to retrieve 10,000 samples for the radius per point in time.
We consider the Newtonian model prediction including the two limiting cases of the Laplace pressure (α = 0 and α = 1) in Figure 13.In the case of α = 0 (no Laplace pressure), we observe that the model prediction increasingly deviates from the experimental data as time proceeds.Moreover, the difference between the model prediction without Laplace pressure and the experimental data decreases with increased applied force.These observations can be explained by the capillary number Ca = η ḣ/γ.In the squeeze flow setting, the radius increases with time meaning that ḣ will decrease and thus also Ca.The capillary number decreases from ∼ 10 −2 at t = 0 s to ∼ 10 −7 at around t = 350 s in case I.A low capillary number indicates that the influence of the surface tension increases.This explains why the model prediction without Laplace pressure describes the short-term flow behavior well in comparison to the longterm flow behavior.It also explains why there is an increasing similarity between the experimental data and model prediction for an increasing applied force.
In Figure 13, we can see that for α = 1 (maximum Laplace pressure) for cases I to IV, denoted in yellow, the model response clearly underestimates the experimental data.Based on the observation in Figure 13, the curvature should be smaller than the maximum curvature to describe the experimental data.Because we are uncertain about the value of α throughout the  experiment, we have a relatively wide range of possible model predictions.In order to quantify the curvature we use the method of Bayesian inference.

Bayesian inference
In the framework of Bayesian inference, we use the experimental squeeze flow data to define the likelihood and use the parametric uncertainties as priors.The likelihood function is defined in a similarly way as the one used in Bayesian inference on the rheological measurements (see Section 5), where we use the experimental data to determine the noise in the likelihood function.However, in this case the noise differs per time step instead of per shear rate.The procedure of Bayesian inference is conducted using the MCMC method described in Section 5.
The prior distribution of the model input parameters is similar to the distribution used in the method of uncertainty propagation.For these parameters, we choose normal distributions.The parameters η, R 0 and V should physically be non-negative.Therefore, we use equation (35) to obtain a lognormal distribution.We use the transformation explained in Remark 3 to obtain the distribution of the true parameter values.We define α in between zero and one.Based on the observations made in Figure 13, we expect that α is closer to zero than to one.Therefore, we assign a beta distribution to α, of which the probability density function is defined as where a and b are the shape parameters valued a = 3 and b = 8, respectively and B(a, b) is a normalization constant.We evaluate the Newtonian squeeze flow model using twelve walkers/chains and 10,000 samples per walker.The burn-in period per walker is 2 × τ f and the thinning is τ f /2, similar to the formulation as described in Section 5.The autocorrelation time per model parameter for each of the cases I to IV is given in Table 9.The effective sample size per case is given in Table 10.
In Figure 14, we show the posterior predictive distribution of the model response R(t) and the experimental data for cases I to IV.The width of the model prediction is comparable to the width of the experimental data.The initial data point in case 0.02 0.04

R[m]
Case I : III shows an increased distribution width for the posterior predictive distribution as well as the experimental data, because the noise in the posterior predictive distribution is determined using the standard deviation of the experimental data, which differs per time step.The short-term flow behavior predicted using uncertainty propagation shows similar behavior to the prediction made using Bayesian inference, because the Laplace pressure's influence decreases toward the short-term flow behavior.Furthermore, the similarities in model response determined using uncertainty propagation and Bayesian inference indicates that the uncertainty in the model input parameters are well-quantified, meaning that the parametric uncertainty determined using independent parameter experiments shows a similar amount of uncertainty as the parametric uncertainty obtained using the squeeze flow experimental data.
The posterior of the model parameters and the correlation between them are visualized in Figure 15 for case I. Similar plots for the remaining cases can be found in the supplementary information.The posterior per input parameter is visualized by the histogram and the prior distribution by the green graph.Furthermore, the area in between the dashed black lines indicates the 95% credibility interval.We observe that for the parameters η N , F and γ the posterior looks similar to the prior distribution, indicating that the model parameters are not updated through the likelihood function.We also observe that the prior of the fluid volume is not similar to the posterior.This suggests that the pipetting and subsequent weighting of the sample yield predicted volumes that are inconsistent with the squeeze flow model predictions.The posterior of α is not similar to the prior distribution and is close to zero, suggesting that the influence of the Laplace pressure is minimally present.We would expect a perfect correlation between α and γ based on the expression for the Laplace pressure (see equation ( 9) and equation ( 36)).However, because we are much more certain about γ in comparison to α, the uncertainty is dominated by γ, due to which the parameters do not seem to be correlated.The posteriors in the corner plot are given as θ for the parameters to which we have applied the log-transformation.The mean, standard deviation and coefficient of variation are given in Table 11 for the actual model parameters θ.For the model parameters that are difficult to determine, such as the curvature-related parameter α, Bayesian inference can give an accurate probablistic prediction.The uncertainty in viscosity and surface tension is expected to be similar for cases I to IV.In Appendix C we show the uncertainty for the remaining cases (cases II to IV), from which it can be concluded that the uncertainty per case is comparable.

Generalized Newtonian fluid
We now consider the comparison between the model predication and experimental data for PVP.The results for the experimental cases V to VIII are visualized in Figures 16 and 17.We observe that cases V and VI show a larger initial radius in comparison to cases VII and VIII, which can be explained by the increased fluid volume in cases V and VI.Comparing cases V and VI, we see that the radius in the latter case increases faster than in the former, which is related to the larger applied force in case VI.Similar observations and conclusions are made between cases VII and VIII, where the latter shows the larger increase in radius.The uncertainty in experimental data increases over time as shown in Figure 16 and 17, similar to the observations made for the experimental data of glycerol (see Figure 13).

Uncertainty propagation
The parametric uncertainties of F, V, R 0 and γ have been determined through additional experiments and the uncertainty in the rheological parameters λ, η 0 , η ∞ and n have been determined using Bayesian inference applied on the simple shear measurements.These uncertainties are propagated through the squeeze flow model to predict the outward motion of the fluid.We have used 12,240 samples from each distribution of the model parameters to retrieve 12,240 samples for the radius per point in time.
We start by comparing the Newtonian squeeze flow model to the experimental data obtained for PVP solution, including and excluding Laplace pressure.We choose to implement the viscosity at the zero-shear rate plateau as input for the Newtonian squeeze flow model.We expect to observe a mismatch between the data and model prediction, because of the shear thinning flow behavior at high shear rates (see Figure 11).
In Figure 16, we show the Newtonian model prediction and experimental data for PVP.There is no significant distinction between the model prediction including and excluding the Laplace pressure.Furthermore, we see that initially the experimental data is underestimated by the model prediction.This can be explained by the constant viscosity at the zero shear rate plateau.Due to the high initial shear rates, the viscosity is lower than the viscosity at the zero-rate plateau, meaning that the evolution rate of the fluid front is higher.We expect that the model prediction improves as we incorporate shear thinning flow behavior in the model through the truncated power law model.The truncated power law model prediction for α = 0 and α = 1 is visualized in Figure 17.We observe that for the initial stages the model prediction has improved using the truncated power law model instead of the Newtonian model.As time proceeds, the model prediction describes the experimental data similarly as the Newtonian model prediction, which can be explained by the decrease in shear rate.Around t = 1 s for cases V to VIII, the model overestimates the experimental data, independent of the inclusion of the Laplace pressure.
We can use Bayesian inference to update the rheological parameters using the rheological data as prior knowledge.In this way, we can investigate whether the rheological parameters as determined in simple flow are representative for the squeeze flow.Furthermore, we will use this method to investigate the uncertainty in curvature of the interface at the fluid front.We will not apply Bayesian inference on the Newtonian model separately, since it is a special case (i.e., specific parameter setting) of the truncated power law model.

Bayesian inference
We apply the framework of Bayesian inference, using MCMC, on the squeeze flow model with the truncated power law model incorporated.The prior distribution of the model input parameters are similar to the ones used in uncertainty propagation.The rheological parameters, η 0 , η ∞ , n and λ cr , should be non-negative and therefore are assigned a lognormal distribution (Remark 3).The surface tension and applied force have a normal distribution and we use a lognormal distribution like the one used for the rheological parameters for the fluid volume and initial radius.The prior distribution of the curvature is a Beta-distribution, given by equation (37), where a = 3, b = 8.
We use 27 walkers/chains, 10,000 samples per walker, a burn-in period of 2 × τ f and τ f /2 for the thinning to apply Bayesian inference on the squeeze flow model.The autocorrelation time for cases V to VIII is provided in Table 12 per model parameter.The effective sample size is given in Table 13 for cases V to VIII.
In Figure 18, the experimental cases V to VIII and model predictions using the Bayesian inference framework are given.The mismatch between the experimental data and model prediction around t = 1 s is significantly decreased in the model prediction using Bayesian inference, meaning that the data available through the likelihood improved the model prediction.To get a closer look at the posterior distribution for the individual parameters, we analyze the two-dimensional posterior marginals provided in Figure 19.The posterior distributions are given by the histograms on the far right image of every row, including the prior distribution given in red.We observe that the posterior distribution corresponds to the prior distribution for the applied force F and rheological parameters η0 and η∞ and reasonably corresponds to the rheological parameters n and λcr .The posterior distribution of γ is similar to the prior distribution and the prior distribution does not correspond to the posterior distribution of α.Note that α is shifted toward zero in comparison to the prior distribution, but not as strong as in the glycerol case (Figure 15).The posterior of the initial radius corresponds to the prior distribution, where the posterior distribution is more narrow than the prior distribution and the mean remains approximately the same.Using Bayesian inference led to a decrease in uncertainty in the initial radius.A significant mismatch can be observed in the prior and posterior distribution of the fluid volume.The mismatch in prior and posterior distribution of the parameters underlies the improved prediction using Bayesian inference around t = 1 s in Figure 18  In analyzing the correlation between the parameters, we investigate the two-dimensional posterior marginals provided in Figure 19.We observe that nearly every posterior marginal does not show a significant correlation between parameters.The posterior marginal between the applied force F and fluid volume V is elliptic, which could indicate a correlation between the parameters or a difference in level of informativeness in the prior.For the parameters to which we have applied the logtransformation, the posteriors are given by θ in the corner plot.The statistical moments are given in Table 14 for the actual model parameters θ.The uncertainty in the rheological parameters and surface tension is expected to be similar for cases V to VIII.In Appendix C we show the uncertainty in the case independent parameters (η 0 , η ∞ , n, λ cr and γ) for the remaining cases (VI to VIII), from which it can be concluded that the uncertainty per case is comparable as expected.

Model-based results
In this section we analyze model-based results for the PVP solution, i.e., results acquired for non-observable quantities.We analyze two model-based results at t = 1.0 × 10 −2 s, t = 1.0 s, t = 10 s and t = 30 s: the local viscosity regimes and the velocity profile.The results are shown for case V. Using the data available in the supplementary information, the results can be obtained for cases VI to VIII as well.
First, we discuss the local viscosity regimes.In Section 4, we have discussed the three regimes in the truncated power law model: Newtonian at zero-shear rate plateau, power law at intermediate shear rates and Newtonian at infinite shear rates.Using the squeeze flow model, we can obtain the interfaces w 1 (t, r), w 2 (t, r) and h(t).In addition, we can obtain the uncertainty in the position of the interfaces.In Figure 20, we show each of the interfaces with the uncertainty.Note that the domain is half of the height of the fluid layer.In this case, we do not reach a high enough shear rate to reach the third regime, where η = η ∞ , therefore, w 2 (r, t) = h(t), due to which we cannot distinguish these two regimes in the figure.In time, we observe a decrease in height and increase in radius as expected.By comparing the local viscosity regime for t = 1.0 × 10 −2 s to t = 30 s we observe that the area covered by the Newtonian regime, where η = η 0 , increases, whereas the power law regime decreases.This shows that the shear thinning behavior decreases through time.The uncertainty in w 1 (r, t) is not constant in the radial domain.In time, the level of uncertainty in h(t) decreases.
In Figure 21 we show the velocity profile at r = R for t = 1.0 × 10 −2 s, t = 1.0 s, t = 10 s and t = 30 s.Note that we show the full height of the fluid layer.We observe a rapid decrease in velocity between t = 1.0 × 10 −2 s and t = 1.0 s, after which the velocity slightly decreases.Furthermore, the height decreases over time in the squeeze flow and we show that the height decreases the fastest in the initial stages.Near the walls, the uncertainty is equal to zero because we assume a no-slip condition.Toward the middle of the velocity profile, the level of uncertainty increases.We expect the profile to change with time, because the shear thinning behavior is only present during the initial stages.By comparing the model prediction of the velocity for the four different points in time, we do not directly observe the changes in the form of the velocity profile, which we attribute to the logarithmic scale and/or the decrease in height.The level of uncertainty decreases in time, especially from t = 1.0 × 10 −2 s to t = 1.0 s.

Conclusions and recommendations
We have used UQ to make probabilistic predictions for generalized Newtonian fluids in squeeze flow.Two methods to quantify the uncertainty in model response have been applied: uncertainty propagation and Bayesian inference.For uncertainty propagation, the parametric uncertainty is propagated through the squeeze flow model to obtain the uncertainty in model response.For Bayesian inference, the experimental squeeze flow data is used to update the model parameters, treated as probabilistic in this framework, and thereby obtain an improved model prediction.To enable the costly inference step using MCMC sampling, we have developed a semianalytical model, assuming a three-regime truncated power law for the viscosity.The UQ framework allows us to separate measurement noise and model bias, and thus make decisions as to which model is more suitable for the given experimental observations.We have seen that for a PVP solution during the initial stages of the squeeze flow, shear thinning behavior is significant, due to which the truncated power law prediction is superior to the Newtonian prediction.Furthermore, the Laplace pressure boundary condition at the fluid front should be included for improved predictions during the later stages.Moreover, we can quantify model-based predictions for quantities that are not directly observable, such as the velocity profile and the local viscosity regimes.With the obtained results, we demonstrate the capabilities of the UQ framework in the field of non-Newtonian fluid mechanics.
Because the concept of uncertainty quantification of non-Newtonian fluid mechanics is in its infancy, we recognize that many paths for future work are possible.For example, an optimization of the MCMC sampling method could lead to a significant reduction in computational time and improve the accuracy of the uncertainty in model parameters and model predictions.Furthermore, efficient parallel sampling methods would be essential in applying UQ using more advanced numerical models.With respect to the obtained squeeze flow result, advanced  rheological models and fluids can be incorporated, such as viscoelastic and viscoplastic models and fluids.Other model extensions that could be incorporated are, e.g., wall-slip and thermal effects.This raises the question whether a more complex model is justified given the experimental data, which we intend to adress in future work using Bayesian model selection.Noting that this expression essentially yields the semi-height rate as a function of the pressure-gradient field, equation (B.3) in essence has the form of a fixed point iteration, which can be solved for the pressure gradient field, and, through (B.4) for the pressure field.
To evaluate the integral in equation (B.4) we employ the midpoint rule on a discretization of the domain [0, R] in n int segments.We use backward Euler time integration, which implies that the radius R in the above derivation also varies during the fixed point iterations.
The initial time step is determined based on the Newtonian model in the high shear rate regime (η = η ∞ ) as where s is a scaling parameter (typically s = 0.001).We use an adaptive time step because in the initial stages of the squeeze flow simulation the time step is required to be much smaller than in the later stages due to high shear rates.The new time step ∆t i+1 is dependent on the Picard iterations and the current time step ∆t i as ∆t i+1 = ∆t i n target n iter , (B.8) where n target is the number of desired Picard iterations (typically n target = 20) and n iter is the number of Picard iterations used in the previous time step.

Convergence study
We here investigate the mesh convergence and time step convergence of our non-linear solver.In Table B.15 the truncated power law parameter settings used for the convergence simulations are listed.
For the convergence study, we define the error for the pressure field as e = p − p h , (B.9 where p is the analytical solution to which the approximate pressure solution p h is compared.Because no analytical solution can be obtained for the pressure, p is defined using a very fine mesh (16,384    elements.We investigate the H 1 -error of the pressure field for the number of time steps and the number of elements to analyze the convergence in the pressure p and the gradient of the pressure ∂p ∂r , given by In this section we evaluate the uncertainty obtained for the squeeze flow model parameters that are case independent using  Bayesian inference.We expect that the uncertainty in rheological parameters and surface tension for cases I to IV is similar, even though we applied Bayesian inference on the separate cases.The same hypothesis holds for cases V to VIII.We show the uncertainty in the Newtonian squeeze flow model parameters in Table C.17 and the truncated power law squeeze flow model parameters in Table C.18.We observe a similar uncertainty in each of the model parameters.

F
add .Due to the high viscosity fluids being used, determining the deposited volume is challenging.The intended volume for glycerol is indicated by markings on the pipette, represented by V a and V b , where V a < V b .For the PVP solution, the intended volume is indicated by drawing two differently sized circles on the PMMA plate denoted by V c and V d , where V c < V d .The values of these fluid volumes and their uncertainty are discussed in Section 5.

Figure 3 :
Figure 3: Overview of the squeeze flow model, where two parallel plates (darkred) compress a fluid with volume V (light blue) by applying a force F. The distance between the parallel plates H = 2h and the radius of the fluid layer R evolve over time.The model domain is denoted by the dashed black box and cylindrical coordinates are used for solving the problem.

Figure 5 :
Figure 5: Schematic of the characterization of the three possible flow regions, where we have three different constitutive models: 1) (pink) η = η 0 , 2) (yellow) η = η 0 (λ cr |γ|) n−1 and 3) (green) η = η ∞ .(a)The division of the regions in simple shear (note that the axes are assumed logarithmic in this schematic), (b) The division of the regions in a squeeze flow experiment at a fixed moment in time.The heights of the interfaces between the regions are denoted by w 1 and w 2 , which are both a function of r.These interfaces are shown in dashed purple and dark red, respectively.

Figure 7 :
Figure 7: Visualization of a stretch move in the parameter space (θ 1 , θ 2 ).The pink dots represent the current coordinates of each walker.The proposal step for walker k, θk , shown in blue, is based on the current position of walkers k and j k.

Figure 8 :
Figure 8: Trace plot of of the log-transformed Newtonian viscosity ηN using four walkers and 10,000 samples.

Figure 9 :Figure 10 :
Figure 9: Posterior predictive distribution of η N for glycerol using a Newtonian constitutive model, where the subscript 'm' corresponds to the model and the subscript 'exp' to the experimental data.Because of the logarithmic scale, the error bars are asymmetric.

Figure 12 :
Figure 12: Corner plots showing the correlation between the truncated power law parameters.The histogram and dark-red graphs denote the posterior and prior, respectively.The region in between the black dashed lines defines the 95% credibility interval.

Figure 13 :
Figure 13: Comparison between the experimental data of glycerol and the model prediction using uncertainty propagation.The subscript 'm' denotes the model prediction and the subscript 'exp' denotes the experimental data.The model prediction including the Laplace pressure is the case where we have the maximum possible curvature.The two extreme cases encompass the experimental data for all cases (I to IV).

Figure 14 :
Figure 14: Comparison between the posterior predictive distribution and experimental data of the squeeze flow using Bayesian inference.The subscript 'm' corresponds the model and the subscript 'exp' to the experimental data.

αFigure 15 :Figure 16 :
Figure 15: Corner plot showing the correlation between the posterior of the model parameters for case I, where we use glycerol.The histogram and dark-red graphs denote the posterior and prior, respectively.The region in between the black dashed lines defines the 95% credibility interval.

Figure 17 :
Figure 17: Truncated power law model prediction compared to the experimental squeeze flow data of PVP using uncertainty propagation.The subscript 'm' denotes to the model and the subscript 'exp' to the experimental data

Figure 18 :
Figure 18: Truncated power law model prediction compared to the experimental squeeze flow data of PVP using Bayesian inference.The subscript 'm' corresponds to the model and the subscript 'exp' to the experimental data in comparison to Figure 17 for cases V to VIII.

4 αFigure 19 :
Figure 19: Corner plot showing the correlation between the posterior of every model parameter for case V, where we use PVP solution.The histogram and dark-red graphs denote the posterior and prior, respectively.The region in between the black dashed lines defines the 95% credibility interval.

Figure 20 :
Figure20: Uncertainty in the local viscosity regimes of case V for various points in time.The error bars cover the µ ± 2σ region.Due to the relatively small uncertainty in r in comparison to the uncertainty in w 1 (r, t), w 2 (r, t) and h(t), we plot the mean of the distribution of r.

Figure B. 22 :
Figure B.22: Convergence study based on the pressure error ||e|| H 1 ; a) time step convergence, b) mesh convergence.

||e|| H 1
.22a the time step convergence is presented.The rate of convergence deviates for less than two time steps, which is caused by the time step being too coarse, resulting in preasymptotic behavior.In FigureB.22b the mesh convergence is visualized.We observe asymptotic convergence behavior for both the mesh size and time step size.A detailed convergence study is beyond the scope of this work.Appendix C. Uncertainty in the squeeze flow model parameters using Bayesian inference

Table 2 :
Uncertainty in the rheological parameter of glycerol.

Table 3 :
Autocorrelation times per model parameter for PVP solution using eight walkers and 20,000 samples per walker.

Table 4 :
Uncertainty in rheological parameters of PVP solution.

Table 5 :
Uncertainty in fluid volume used for the squeeze flow cases.V a is used for cases I, II and III.V b is used for case IV (see also Table1).

Table 6 :
Uncertainty in applied force for three different additions of weight.F add,1 , F add,2 and F add,3 correspond to 0.0 kg, 0.25 kg and 0.5 kg, respectively.

Table 7 :
Uncertainty in initial radius for the squeeze flow cases.

Table 8 :
Uncertainty in surface tension.
and curvature of the interface κ.The uncertainty in surface tension of glycerol has been measured by pendant drop experiments, of which the mean, standard deviation and coefficient of variation are provided in Table

Table 9 :
Autocorrelation times per model parameter and case using twelve walkers and 10,000 samples per walker.

Table 10 :
Effective sample size per case based on the highest autocorrelation time per case for a sample size N = 120, 000.

Table 11 :
Uncertainty in the Newtonian squeeze flow model parameters θ for case I.

Table 12 :
Autocorrelation times per model parameter and case using 27 walkers and 10,000 samples per walker.

Table 13 :
Effective sample size N final per case based on the highest autocorrelation time per case, for sample size N = 270, 000.

Table 14 :
Uncertainty in the truncated power law squeeze flow model parameters θ for case V.

Table B .
15:Model input parameters for a truncated power law model used in the convergence study.
(1,ments) and a high number of time steps(1,024).In Table B.16a, five simulation inputs are given with an increase in number of time steps.In Table B.16b, five simulation inputs are provided with an increase in number of Table B.16: Convergence study inputs for the number of time steps n t and number of elements m.

Table C .
17: Uncertainty in case independent Newtonian squeeze flow model parameters θ for the cases II to IV.