Learning Unknown Physics of non-Newtonian Fluids
Brandon Reyes, Amanda A. Howard, Paris Perdikaris, Alexandre M. Tartakovsky
LLearning Unknown Physics of non-Newtonian Fluids
Brandon Reyes and Amanda A. Howard
Pacific Northwest National Laboratory, Richland, WA
Paris Perdikaris
University of Pennsylvania, Philadelphia, PA
Alexandre M. Tartakovsky ∗ Pacific Northwest National Laboratory, Richland,WA; Department of Civil and Environmental Engineering,University of Illinois Urbana-Champaign, Urbana, IL (Dated: September 4, 2020)We extend the physics-informed neural network (PINN) method to learn viscosity models of twonon-Newtonian systems (polymer melts and suspensions of particles) using only velocity measure-ments. The PINN-inferred viscosity models agree with the empirical models for shear rates withlarge absolute values but deviate for shear rates near zero where the analytical models have anunphysical singularity. Once a viscosity model is learned, we use the PINN method to solve themomentum conservation equation for non-Newtonian fluid flow using only the boundary conditions.
INTRODUCTION
In many applications, data is scarce and indirect andthe governing physics is not fully known, which limits theutility of standard machine learning (ML) and physics-based methods. For example, in non-Newtonian flow ex-periments one can easily measure velocity, but not stressor viscosity. This makes it impossible to use data-drivenML methods to learn stress as a function of velocityor shear rate. Also, the momentum and mass conser-vation equations governing non-Newtonian flow are notfully known as one needs to assume a stress-shear-raterelationship (we refer to such relationships as unknownphysics) to close the system of these equations. It isimportant to note that standard parameter estimationmethods cannot be used for learning unknown physics be-cause the function space is infinite-dimensional. It is thisissue that the physics informed neural network (PINN)method attempts to solve. PINNs use the known un-derlying structure of physical laws governed by PDEsor ODEs to predict unknown functions or functionalsfrom indirect observations. By representing states ofthe system and hidden physics with neural networks andtraining using available data subject to the conservationlaws, the PINN method can learn unknown physics usingsparse and indirect data. In the past, the PINN methodwas used to learn unknown physics in partially unsat-urated flow in porous media [1] and bioreactors [2]. Inthis work, we extend the PINN method for estimating thenon-Newtonian viscosity based solely on velocity data.
PINN METHOD FOR NON-NEWTONIAN FLOWMODELS
Consider a shear flow of a non-Newtonian fluid betweentwo parallel plates satisfying the steady-state momentumconservation equation: ddy (cid:20) µ ( u y ( y )) du ( y ) dy (cid:21) = − C for y ∈ Ω = (0 , H ) , (1)where the velocity vector is given by u = ( u ( y ) , , T , u y ≡ du/dy is the shear rate, the viscous stress hasthe form µ ( u y ) u y ( y ), µ ( u y ) is the unknown shear-rate-dependent viscosity, H is the channel width, and C is aforce per unit volume. The fluid velocity u is subject tothe no-slip boundary conditions (BCs): u (0) = 0 , u ( H ) = 0 . (2)We consider two cases: no measurements of µ are avail-able and some measurements of µ are present. In bothcases we assume that there are N u measurements of thevelocity profile u ( y ) for y ∈ Ω: u ∗ ( y i ) for i = 1 , . . . , N u .We approximate the viscosity µ ( u y ) and the velocity u ( y )with fully connected feed-forward deep neural networks(DNNs), u ( y ) ≈ ˆ u ( y ; θ ) and µ ( u y ( y )) ≈ ˆ µ (ˆ u y ( y ; θ ); γ ),where θ and γ are the DNN weights. We train ˆ u and ˆ µ jointly using Eqs. (1) and (2) as constraints. This allowsus to train ˆ µ even without direct measurements of µ .We note that the DNNs ˆ u and ˆ µ are known non-linearfunctions of y and θ and/or γ . Therefore, we can ana-lytically compute the DNN derivatives with respect to y and the weights. The former are needed to impose thephysical constraints given by Eq. (1), while the latterare required to update the values for the weights in theprocess known as backpropagation [ ? ]. Here, we use au-tomatic differentiation [ ? ] to compute the derivatives. a r X i v : . [ phy s i c s . c o m p - ph ] A ug Eq. (1) is enforced in the DNN training by forming anadditional “auxiliary” DNN:ˆ f ( y ; θ, γ ) = ddy (cid:20) ˆ µ (ˆ u y ( y ; θ ); γ ) d ˆ u ( y ; θ ) dy (cid:21) . (3)We train the DNNs simultaneously by minimizing theloss function L ( θ, γ ) = ω N u N u (cid:88) i =1 [ˆ u ( y i ; θ ) − u ∗ ( y i )] (4)+ ω (cid:2) ˆ u ( y = 0; θ ) + ˆ u ( y = H ; θ ) (cid:3) + ω N u N u (cid:88) i =1 (cid:104) ˆ f ( y i ; θ, γ ) + C (cid:105) + ω N µ N µ (cid:88) i =1 (cid:2) ˆ µ ( u yi ; γ ) − µ ∗ ( u yi ) (cid:3) . In L ( θ, γ ), the first term forces ˆ u ( y ; θ ) to match the ve-locity measurements, the second term forces ˆ u ( y ; θ ) tomatch the Dirichlet BCs, and the third term forces ˆ u andˆ µ to satisfy Eq. (1). The last term is present ( ω (cid:54) = 0)if measurements of µ (i.e., µ ∗ ( u yi ) for i = 1 , ..., N µ ) areavailable and forces ˆ µ to match these measurements. Theweights { ω i } i =1 ,..., reflect the fidelity level of the dataand physics models. For example, u measurements aremore accurate than viscosity measurements in general,so ω ≥ ω . We note that Eq. (1) is an approximation ofthe momentum conservation equation because it involvesassumptions about the general form of the viscous stress,therefore, ω ≤ ω . For some flows the no-slip BCs as-sumption might not be very accurate, which would affectthe relative value of ω . The relative values of ω i canalso affect the convergence rate of iterative solutions ofthe minimization problem ( θ, γ ) = arg min θ,γ L ( θ, γ ) [3 ? ]. To solve this minimization problem, we set the initialvalues of θ and γ using the Xavier’s normal initializa-tion scheme [4]. Next, we run the Adam optimizer [5]for a set number of steps. Finally, we run the quasi-Newton L-BFGS-B optimizer [6] until the desired con-vergence and tolerance are achieved. We find that forthe considered here problems, this combination of the op-timizers increases the convergence rate and reduces thecomputational cost as compared to using either optimizeralone. We use DNNs with two hidden layers with sixtynodes each and a learning rate of 0.001 for the Adam op-timizer. The error (cid:107) e ˆ u (cid:107) = (cid:107) ˆ u ( y ; θ ) − u ∗ ( y ) (cid:107) / (cid:107) u ∗ ( y ) (cid:107) estimates the accuracy of the DNN approximations of u relative to the u measurements and the error (cid:107) f (cid:107) ∞ =max ≤ i ≤ N u | f ( y i ; θ, γ ) + C | is a measure of how well theDNN approximations of u and µ satisfy Eq. (1).We refer to the PINN method that is used to evaluatethe unknown viscosity function given the measurements u (or u and µ ) as the inverse PINN. Once ˆ µ is trained, the PINN method can also be used to solve the momentumconservation equation without observations of u (and/or µ ) if the shear rate does not exceed the maximum shearrate in the experiment used to train ˆ µ . To train ˆ u as anapproximate solution of Eq. (1) we use the loss functionEq. (4) with ω = ω = 0 and ω = ω = 1. We refer tothis application of PINNs as the forward PINN method. VALIDATION OF THE INVERSE ANDFORWARD PINN METHODS
We first validate the ability of the inverse PINNmethod to learn the unknown shear-dependent viscosityusing velocity data generated with the Ostwald-de Waelepower-law effective viscosity model [7], µ pl ( u y ( y )) = K | u y ( y ) | n − , where K is the power-law consistency co-efficient and n is the power-law index. This model incombination with Eqs. (1) and (2) allows for an analyti-cal solution for u ( y ) and du ( y ) /dy [8].We generate two data sets by selecting N u = 501 uni-formly distributed measurements of u from the analyti-cal solution for u using both n = 0 .
898 (shear-thinningfluid) and n = 1 . C = 0 . H = 25, and K = 40 . n we trainthe ˆ u ( y ; θ ) and ˆ µ (ˆ u y ( y ; θ ); γ ) DNNs by minimizing theloss function Eq. (4) with ω = ω = ω = 1 and ω = 0.We note that the minimization problem is not convexand its solution ( θ, γ ) can depend on the initial valuesof θ and γ . To demonstrate how different initial valuesfor the weights affect the PINN solution, we solve theminimization problem with 100 different initializationsof θ and γ and then average the resulting DNNs ˆ u ( y ; θ )and ˆ µ (ˆ u y ( y ; θ ); γ ) to obtain the solutions for u ( y ) and µ ( u y ), respectively. For n = 0 .
898 the average solutionsare compared with the analytical solutions in Figs. 1aand 1b. The average DNN ˆ u ( y ) solution agrees well withthe analytical u ( y ) solution. The average DNN ˆ µ (ˆ u y ) so-lution agrees very well the constitutive model for largeshear rates. For small shear rates, the DNN solutiondeviates from the analytical solution and for zero shearrate has a finite value while the analytical solution has anonphysical singularity. Fig. 1b also shows that the stan-dard deviation in the learned µ ( u y ) is largest at u y = 0and is several orders of magnitude smaller than the meanvalue of µ at u y = 0, indicating that the uncertainty ofthe PINN method due to DNN initialization is relativelysmall. Fig. 1c depicts the residual ˆ f ( y ; θ, γ ) + C of Eq.(1) as a function of y . The small values of the residualshow that the DNNs ˆ u and ˆ µ approximately satisfy Eq.(1). The (cid:107) e ˆ u (cid:107) and (cid:107) f (cid:107) ∞ errors for both values of n are given in Table I. Small errors demonstrate that theinverse PINN method is equally accurate for both shear-thinning and shear-thickening fluids.Next, we validate the ability of the forward PINNmethod to solve Eq. (1). We fix the weights of the DNN
12 131.1451.150 a) c) FIG. 1. Inverse and forward PINN solutions for the syntheticdata generated from the analytical solution for a power-lawfluid with n = 0 . dudy of the 100 runs versus the average µ . The variance is given by the grey area. (c) Average errorin satisfying the ODE. n (cid:107) e ˆ u (cid:107) (cid:107) f (cid:107) ∞ × − × − × − × − TABLE I. Mean errors computed from 100 inverse PINN so-lutions for the synthetic data generated from the analyticalsolution for a power-law fluid with n = 0 .
898 and 1.2. ˆ µ obtained from the inverse PINN with n = 0 .
898 andtrain the ˆ u ( y, θ ) DNN by minimizing the loss functionEq. (4) with ω = ω = 0, ω = ω = 1. Fig. 1a showsthat the trained ˆ u ( y, θ ) closely agrees with the analyti-cal solution for the power-law fluid with n = 0 . C , indi-cating that ˆ u approximately satisfies Eq. (1). The goodagreement between ˆ u and the reference solution for u andsmall residuals confirm the accuracy of the forward PINNmethod for solving non-linear differential equations withconstitutive relationships given by a DNN with knownweights. MONODISPERSE POLYMER MELTS
We consider a synthetic Dissipative Particle Dynamic(DPD) fluid consisting of chains of N equal-size beadsconnected by springs to model polymer melts. Two-dimensional DPD simulations of such fluids between twoparallel plates with chains made of N = 2, 5, and 25beads are presented in [9]. In [9], the DPD resultswere used to compute µ DP D ( du ( y ) /dy ) using the Irving-Kirkwood relationship [10].We use the velocity data from [9] and the inverse PINNmethod with ω = ω = ω = 1 and ω = 0 in Eq.(4) to estimate µ ( u y ). To match [9] , C = 0 .
75 and H = 25. The relative velocity error and the maximumresidual error are given in Table II. For all considered N the relative error in u is less than 0.1% and the maximumresidual error is 3 orders of magnitude smaller than the N (cid:107) e ˆ u (cid:107) (cid:107) f (cid:107) ∞ × − × − × − × −
25 4.396 × − × − TABLE II. Errors for N = 2, 5, and 25 beads. driving force C , indicating that the DNN ˆ u accuratelyapproximates data and the DNNs ˆ u and ˆ µ satisfy thegoverning equations.Figs. 2a and 2b compare the velocity profiles andviscosities estimated from the DPD simulation, µ DP D and from the PINN method for N = 2. The DNNvelocity profile ˆ u ( y ; θ ) closely matches the DPD veloc-ity profile u DP D ( y ). The agreement between ˆ µ ( u y ; θ, γ )and µ DP D ( u y ) is good but less accurate than the agree-ment for the velocities. To test whether u DP D ( y ) and µ DP D ( u y ) satisfy Eq. (1), we train the ˆ u ( y ; θ ) andˆ µ ( u y ; θ, γ ) DNNs conditioned on both u and µ DPD mea-surements. Figs. 2a and 2b show that conditioning of theDNNs on the DPD measurements of u and the estimatesof µ produces DNNs that match well both u and µ data.However, conditioning on the DPD µ estimates also re-sults in the residual errors that are two orders of magni-tude larger than the residual errors in the case where no µ estimates are used to train the DNNs, as shown in Fig.2c.Next, we use the PINN method to evaluate the viscos-ity of the polymer melt with 25-bead chains. As for themelt with N = 2, we first train the ˆ u and ˆ µ DNNs usingonly u DP D ( y ) measurements. Fig. 3 shows that the ˆ u DNN agrees well with the u DP D ( y ) measurements andthe resulting residual point errors are nearly zero (morethan four orders of magnitude smaller than C ). We alsosee that ˆ µ significantly deviates from the µ DP D ( u y ) val-ues estimated from the DPD simulations near a shearrate of zero. Then, we train the ˆ u and ˆ µ DNNs usingboth u DP D ( y ) and µ DP D ( y ) data. Fig. 3 demonstratesthat the resulting DNNs fit the u DP D ( y ) and µ DP D ( y )data well, but the corresponding residual is very large(on the order of C ). We obtain similar results for thepolymer melt with N = 5.Finally, we demonstrate that once ˆ µ ( u y ; θ, γ ) istrained, the forward PINN method can be used to solveEq. (1) subject to the BC Eq. (2). We use the weights γ in the ˆ µ ( u y ; θ, γ ) DNN obtained above from the inversePINN and train the forward solution, ˆ u f ( y ; θ ), DNN byminimizing the loss function Eq. (4) with ω = ω = 0and ω = ω = 1 for C = 0 .
75. For N = 2, Fig. 2ashows that the ˆ u f ( y ; θ ) DNN matches the experimentaldata corresponding to C = 0 .
75 well. In addition tothis, Fig. 2c demonstrates that the residual of the gov-erning equation is two orders of magnitude smaller than C confirming that ˆ u f ( y ; θ ) approximately solves Eq. (1)
12 131.1401.1451.150 a) b) c)
FIG. 2. Inverse and forward model results for N = 2. (a)Resulting velocity profiles. (b) Resulting viscosity profile. (c)Error in satisfying the ODE. In the power-law solution wetake n = 0 .
90 and K = 40 .
11 140.7500.7750.800 a) b) c) FIG. 3. PINN results for N = 25. (a) Resulting velocityprofiles. (b) Resulting viscosity profile. (c) Error in satisfyingthe ODE. In the power-law solution we take n = 0 .
79 and K = 43 . subject to Eq. (2).These results lead to the conclusion that the inversePINN is capable of estimating the effective viscosity func-tion µ ( u y ), which can be used for solving the momentumconservation equation (1). The predicted viscosity devi-ates from the viscosity obtained from the DPD simula-tions for small shear rates with the discrepancy increasingwith the number of beads N . Our results show that ve-locity and viscosity data provided in [9] cannot accuratelybe described by Eq. (1). DENSE SUSPENSIONS OF SPHERICALPARTICLES
In this section, we employ the inverse PINN methodto learn the shear-rate-dependent viscosity of denselypacked spherical particles suspended in a Newtonian fluidusing the velocity measurements presented in [11]. Theconsidered data are obtained from the numerical simula-tions of suspension flows in a channel using the ForceCoupling Method (FCM) [11, 16, 18]. In the consid-ered suspensions, the average particle volume fraction φ a = πa NV ranges from 0.2 to 0.4, where a is the parti-cle radius, N is the number of particles, and V is the vol-ume of the domain. In the FCM simulations, the particleradius was set to a = 1, the channel length to L x = 80,the height to H = 40, and the width to L z = 30. Thechannel walls were located at y = 0 and 40, constantDirichlet BCs for pressure were prescribed at the x = 0and 80 boundaries with the pressure drop over the lengthof the channel ∆ P/L x = 0 . z direction. At the continuum level, theconsidered suspension behaves as a non-Newtonian fluidand can be described by Eq. (1) with C = ∆ P/L x .The velocity profiles for the suspension flows with φ a = 0 . φ ( y ) are depictedin 5a. A key feature of suspensions is irreversible shear-induced migration of particles to areas of low shear rate[12]. Particles in a suspension subjected to a Poiseuilleflow will migrate to the channel centerline, greatly in-creasing the volume fraction at the centerline until itreaches the maximum close-packing limit, as shown inFig. 5a. This migration also impacts the velocity profile,resulting in a flattened parabola shape that is observedin Figs. 4a and d.As in the analysis of polymer melts above, we use theinverse PINN to find the viscosity µ ( u y ( y )) by approx-imating u and µ with ˆ u ( y ; θ ) and ˆ µ (ˆ u y ( y ; θ ); γ ) DNNstrained by minimizing the loss function (4) with ω = ω = ω = 1 and ω = 0. We use N u = 401 measure-ments of the velocity profile u ( y ) from the FCM simu-lations. Because the velocity profiles from the simula-tions (see Fig. 4a and d) deviate from the flattened-parabola shape near the walls due to particle layer-ing, a phenomenon that cannot be described by Eq.(1), we train the PINN with velocity data in the range y ∈ [0 . y/h, . y/h ], but still impose the zero DirichletBCs for u at y = 0 , H .Figs. 4a and 4d compare the velocity profiles of the sus-pension flow observed in the numerical simulations andare approximated with the ˆ u ( y ; θ ) DNN for φ a = 0 . µ (ˆ u y ( y ; θ ); γ ) DNN and the viscos-ity estimated from the numerical experiments are plottedin Figs. 4b and 4e. The viscosity µ ( u y ) for the FCMsimulations is found by computing u y ( y ) and φ ( y ) fromthe simulation data, assuming that µ ( u y ) = η s ( φ ( u y )) η f and using the Eilers formula η s ( φ ) = (cid:18) φ ( − φφc ) (cid:19) [19, 20]. Here, η f is the fluid viscosity (which was set tounity in the FCM simulations) and φ c is the maximumvolume fraction of a suspension ( φ c = 0 .
62 in the FCMsimulations.) We observe that the PINN method is ableto accurately learn the velocity profile and captures theincrease in viscosity at the channel centerline. Figs. 4cand 4f demonstrate that the residuals are three orderssmaller than C = 0 . u ( y ; θ ) andˆ µ (ˆ u y ( y ; θ ); γ ) DNNs satisfy Eq. (1).Finally, we employ the inverse PINN method to eval-uate η s as a function of φ and in Fig. 5b compare itwith the Eilers model and the Krieger model η s ( φ ) = (cid:16) − φφ c (cid:17) − . φ c [21]. In the PINN method, we compute η s ( φ ) using the ˆ µ (ˆ u y ( y ; θ ); γ ) and ˆ u ( y ; θ ) DNN models ofviscosity and velocity and φ ( y ) observed in the FCM sim-ulations. The considered empirical models predict similar µ values for φ < .
35 away from the channel centerline. a) b) y / h f ( y ) + C
1e 5 c) d) e) y / h f ( y ) + C
1e 5 f) FIG. 4. Inverse PINN results for suspensions with averagevolume fraction φ a = 0 . . y / h ( y ) a =0.2 a =0.3 a =0.4 a) EilersKrieger a =0.2 a =0.3 a =0.4 b) FIG. 5. a) Final suspension local volume fraction profiles φ ( y )in steady-state. b) Suspension viscosities learned from thePINN model as a function of the local volume fraction. Re-sults are compared with the Eilers fit [19, 20] and the Kriegerfit [21]. Filled symbols represent points that occur in therange 0 h ≤ y ≤ . h , and empty symbols are in the range0 . h ≤ y ≤ h , to denote the deviations that occur from thetheoretical values in the densely packed core of the channel. The empirical models assume that µ ( φ ) is independentof φ a . Fig. 5b shows that the PINN predicted µ ( φ ) func-tions agree with the empirical models for small φ for allconsidered φ a . For large φ , the PINN estimated µ ( φ )relationships depend on φ a and deviate from all consid-ered empirical models. There are several reasons for thedisagreement between the PINN estimated and empiri-cal viscosity models. At high volume fractions near theclose-packing limit, φ c , the particle movements are highlycorrelated leading to non-locality of the particle forces.Therefore, Eq. (1) breaks down for dense suspensions atthe centerline. Additionally, Eq. (1) with a constitutiverelationship of the Eilers or Krieger analytical forms pre-dicts that the suspension will reach maximum packingwith φ = φ c at the centerline that is independent of theaverage volume fraction φ a [22]. However, FCM simu-lations [16] and experiments [13] show that the volumeat the centerline varies with the initial average volumefraction of the system. In conclusion, we have extended the PINN method forlearning unknown physics, including the functional de-pendence of viscosity on the shear rate and other prop-erties of fluids using indirect measurements such as fluidvelocity and volume fraction. We have also demonstratedthat once an accurate DNN approximation of the viscos-ity is available, the PINN method can be used to modelnon-Newtonian flow without any data except the bound-ary conditions. ∗ [email protected][1] A. M. Tartakovsky, C. O. Marrero, P. Perdikaris, G. D.Tartakovsky, and D. Barajas-Solano, Water ResourcesResearch , e2019WR026731 (2020).[2] R. Tipireddy, P. Perdikaris, P. Stinis, and A. Tar-takovsky, arXiv preprint arXiv:1904.04058 (2019).[3] S. Wang, Y. Teng, and P. Perdikaris, arXiv preprintarXiv:2001.04536 (2020).[4] X. Glorot and Y. Bengio, in Proceedings of the Thir-teenth International Conference on Artificial Intelligenceand Statistics , edited by Y. W. Teh and M. Tittering-ton (PMLR, Chia Laguna Resort, Sardinia, Italy, 2010),vol. 9 of
Proceedings of Machine Learning Research , pp.249–256.[5] D. P. Kingma and J. Ba, arXiv preprint arXiv:1412.6980(2014).[6] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, SIAM J. Sci.Comput. , 1190 (1995).[7] R. Bird, W. Stewart, and E. Lightfoot, Transport Phe-nomena , Wiley International edition (Wiley, 2006), ISBN9780470115398.[8] E. Hinch,
Lecture 3: Simple flows , URL .[9] D. A. Fedosov, G. E. Karniadakis, and B. Caswell, TheJournal of Chemical Physics , 144103 (2010).[10] J. Irving and J. G. Kirkwood, The Journal of chemicalphysics , 817 (1950).[11] A. Howard, Ph.D. thesis, Brown University, Providence,RI (2018).[12] D. Leighton and A. Acrivos, J. Fluid Mech. , 415(1987).[13] M. K. Lyon and L. G. Leal, J. Fluid Mech. , 25(1998).[14] J. E. Butler, P. D. Majors, and R. T. Bonnecaze, Phys.Fluids , 2865 (1999).[15] B. Snook, J. E. Butler, and ´E. Guazzelli, J. Fluid Mech. , 128 (2015).[16] K. Yeo and M. R. Maxey, J. Fluid Mech. , 491 (2011).[17] F. Cui, A. Howard, M. Maxey, and A. Tripathi, Phys.Rev. Fluids (2017).[18] K. Yeo and M. R. Maxey, J. Comput. Phys. , 2401(2010).[19] F. Ferrini, D. Ercolani, B. de Cindio, L. Nicodemo,L. Nicolais, and S. Ranaudo, Rheol. Acta , 289 (1979).[20] J. J. Stickel and R. L. Powell, Annu. Rev. Fluid Mech , 129 (2005).[21] I. M. Krieger and T. J. Dougherty, Trans. Soc. Rheol. ,137 (1959).[22] E. Guazzelli and O. Pouliquen, Journal of Fluid Mechan- ics (2018). [23] G. Batchelor and J.-T. Green, Journal of Fluid Mechanics56