title: Learning Unknown Physics of non-Newtonian Fluids
authors: Brandon Reyes Amanda A. Howard Paris Perdikaris Alexandre M. Tartakovsky
text: 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Pacific Northwest National Laboratory, Richland, WA
  
  Pacific Northwest National Laboratory, Richland, WA
  
  University of Pennsylvania, Philadelphia, PA
  
  []amt1998@illinois.edu
  
  
  
  Pacific Northwest National Laboratory, Richland, WA; Department of Civil and Environmental Engineering, University of Illinois Urbana-Champaign, Urbana, IL
  
  
  
  
  
  
  
  
  
  
  
  We extend the physics-informed neural network (PINN) method to learn viscosity models of two non-Newtonian systems (polymer melts and suspensions of particles) using only velocity measurements. The PINN-inferred viscosity models agree with the empirical models for shear rates with large absolute values but deviate for shear rates near zero where the analytical models have an unphysical singularity. Once a viscosity model is learned, we use the PINN method to solve the momentum conservation equation for non-Newtonian fluid flow using only the boundary conditions.
  
  
  
  
  
  
  
  Learning Unknown Physics of  non-Newtonian Fluids
      Alexandre M. Tartakovsky
      December 30, 2023
  =================================================
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  In many applications, data is scarce and indirect and the governing physics is not fully known, which limits the utility of standard machine learning (ML) and physics-based methods. For example, in non-Newtonian flow experiments one can easily measure velocity, but not stress or viscosity.  This makes it impossible to use data-driven ML methods to learn stress  as a function of velocity or shear rate. Also, the momentum and mass conservation equations governing non-Newtonian flow are not fully known as one needs to assume a stress-shear-rate relationship (we refer to such relationships as unknown physics) to close the system of these equations. 
  It is important to note that standard parameter estimation methods cannot be used for learning unknown physics because the function space is infinite-dimensional.    It is this issue that the physics informed neural network (PINN) method attempts to solve. PINNs use the known underlying structure of physical laws governed by PDEs or ODEs to predict unknown functions or functionals from indirect observations. 
  By representing states of the system and hidden physics with neural networks and training using available data subject to the conservation laws, the PINN method can learn unknown physics using sparse and indirect data. 
  In the past, the PINN method was used to learn unknown physics in partially unsaturated flow in porous media <cit.> and bioreactors <cit.>. In this work, we extend the PINN method for estimating the non-Newtonian viscosity based solely on velocity data.     
  
  
  
  
  
  § PINN METHOD FOR NON-NEWTONIAN FLOW MODELS
  
  Consider a shear flow of a non-Newtonian fluid between two parallel plates satisfying the steady-state momentum conservation equation: 
  
      d/dy[ μ(u_y(y)) d u(y)/dy] = -C for y ∈Ω = (0, H),
  
  where the velocity vector is given by 𝐮 = (u(y), 0, 0)^T, u_y ≡ du/dy is the shear rate, the viscous stress has the form μ(u_y) u_y (y), μ(u_y) is the unknown shear-rate-dependent viscosity, H is the channel width, and C is a force per unit volume.  The fluid velocity u is subject to the no-slip boundary conditions (BCs): 
  
      u(0) = 0,  u(H) = 0 .
  
  We consider two cases: no measurements of μ are available and  some measurements of μ are present. In both cases we assume that there are N_u measurements of the velocity 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) ≈û(y; θ)  and  μ(u_y(y)) ≈μ̂(û_y(y;θ); γ),
  
  where θ and γ are the DNN weights. We train û and μ̂ jointly using Eqs. (<ref>) and (<ref>) as constraints. 
  This allows us to train μ̂ even without direct measurements of μ.
  
  We note that the DNNs û and μ̂ are known non-linear functions of y and θ and/or γ. Therefore,  we can analytically compute the DNN derivatives with respect to y and the weights. The former are needed to impose the physical constraints given by Eq. (<ref>), while the latter are required to update the values for the weights in the process known as backpropagation <cit.>. Here, we use automatic differentiation <cit.> to compute the derivatives.
  
  Eq. (<ref>) is enforced in the DNN training by forming an additional “auxiliary" DNN: 
  
      f̂(y; θ, γ) = d/dy[ μ̂(û_y(y;θ); γ) d û (y; θ)  /dy].
  
  We train the DNNs simultaneously by minimizing the loss function
  
  
      L(θ, γ)    = ω_1/N_u∑_i=1^N_u[  û(y_i; θ) - u^* (y_i)  ]^2  
         + ω_2 /2[û(y=0; θ)^2 + û(y=H; θ)^2  ] 
         +ω_3 /N_u∑_i=1^N_u[ f̂(y_i; θ, γ) + C ]^2 
         + ω_4/N_μ∑_i=1^N_μ[  μ̂(u_y_i; γ) - μ^* (u_y_i)  ]^2.
  
  
  In L(θ, γ), the first term forces û (y; θ) to match the velocity measurements, the second term forces û (y; θ) to match the Dirichlet BCs, and the third term forces û and μ̂ to satisfy Eq. (<ref>). The last term is present (ω_4 ≠ 0) if measurements of μ (i.e., μ^* (u_y_i) for i=1,...,N_μ) are available and forces μ̂ to match these measurements. The weights {ω_i}_i = 1, …, 4 reflect the fidelity level of the data and physics models. For example, u measurements are more accurate than viscosity measurements in general, so ω_1 ≥ω_4. We note that Eq. (<ref>) is an approximation of the momentum conservation equation because it involves assumptions about the general form of the viscous stress, therefore, ω_3 ≤ω_1. For some flows the no-slip BCs assumption might not be very accurate, which would affect the relative value of ω_2.  The relative values of ω_i can also affect the convergence rate of iterative solutions of the minimization problem  
  (θ, γ) = argmin_θ, γ L(θ, γ) 
  <cit.>. 
  
  To solve this minimization problem, we set the initial values of θ and γ using the Xavier's normal initialization scheme <cit.>. Next, we run the Adam optimizer <cit.> for a set number of steps. Finally, we run the quasi-Newton L-BFGS-B optimizer <cit.> until the desired convergence and tolerance are achieved. We find that  for the considered here problems, this combination of the optimizers increases the convergence rate and reduces the computational cost as compared to using either optimizer alone. We use DNNs with two hidden layers with sixty nodes each and a learning rate of 0.001 for the Adam optimizer. The error e_û_2 =  û(y; θ) - u^* (y) _2/   u^* (y) _2 estimates the accuracy of the DNN approximations of u relative to the u measurements and the error f _∞ = max_1 ≤ i ≤ N_u |f(y_i; θ, γ) + C | is a measure of how well the DNN approximations of u and μ satisfy Eq. (<ref>). 
  
  We refer to the PINN method that is used to evaluate the 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 momentum conservation equation without observations of u (and/or μ) if the shear rate does not exceed the maximum shear rate in the experiment used to train μ̂. 
  To train û as an approximate solution of Eq. (<ref>) we use the loss function Eq. (<ref>)  with ω_1 = ω_4 = 0 and ω_2 = ω_3 = 1. We refer to this application of PINNs as the forward PINN method. 
  
  
  
  
  
  
  
  
  
  
  § VALIDATION OF THE INVERSE AND FORWARD PINN METHODS
  
  
  We first validate the ability of the inverse PINN method to learn the unknown shear-dependent viscosity using velocity data generated with the Ostwald-de Waele power-law effective viscosity model <cit.>, 
  
  	μ_pl( u_y(y)) = K | u_y(y) | ^n-1,
  
  
  where K is the power-law consistency coefficient and n is the power-law index. This model in combination with Eqs. (<ref>) and (<ref>) allows for an analytical solution for u(y) and du(y)/dy <cit.>. 
  
  We generate two data sets by selecting N_u =501 uniformly distributed measurements of u from the analytical solution for u using both n = 0.898 (shear-thinning fluid) and n = 1.2 (shear-thickening fluid) with C = 0.75, H = 25, and K = 40.788. For both values of n we  train the û(y;θ) and μ̂(û_y(y;θ);γ) DNNs by minimizing the loss function Eq. (<ref>) with ω_1 = ω_2 =ω_3 =1 and ω_4=0. We note that the minimization problem is not convex and its solution (θ,γ) can depend on the initial values of  θ and γ. To demonstrate how different initial values for the weights affect the PINN solution, we solve the minimization problem with 100 different initializations of θ and γ and then average the resulting DNNs û(y;θ) and μ̂(û_y(y;θ);γ) to obtain the solutions for u(y) and μ(u_y), respectively. For n = 0.898 the average solutions are compared with the analytical solutions in Figs. <ref>a and  <ref>b. The average DNN û(y) solution agrees well with the analytical u(y) solution. The average DNN μ̂(û_y) solution agrees very well the constitutive model  for large shear rates. For small shear rates, the DNN solution deviates from the analytical solution and for zero shear rate has a finite value while the analytical solution has a nonphysical singularity. 
   Fig. <ref>b also shows that the standard deviation in the learned μ(u_y) is largest at u_y = 0 and is several orders of magnitude smaller than the mean value of μ at u_y = 0, indicating that the uncertainty of the PINN method due to DNN initialization is relatively small.  
  Fig. <ref>c depicts the residual f̂(y;θ,γ) + C of Eq. (<ref>) as a function of y. The small values of the residual show that the DNNs û and μ̂ approximately satisfy Eq. (<ref>).    
  The e_û_2 and f _∞ errors for both values of n are given in Table <ref>. Small errors demonstrate that the inverse PINN method is equally accurate for both shear-thinning and shear-thickening fluids. 
  
   
  
  
  
  
  Next, we validate the ability of the forward PINN method to solve Eq. (<ref>). 
  We fix the weights of the DNN μ̂ obtained from the inverse PINN with n=0.898 and train the û(y,θ) DNN by minimizing the loss function Eq. (<ref>)  with ω_1 = ω_4 = 0, ω_2 = ω_3 = 1.  Fig. <ref>a shows that the trained û(y,θ) closely agrees with the analytical solution for the power-law fluid with n=0.898. Fig. <ref>c  shows the maximum residual corresponding to this DNN is two orders of magnitude smaller than C, indicating that û approximately satisfies Eq. (<ref>).
  
  The good agreement between û and the reference solution for u and small residuals confirm the accuracy of the forward PINN method for solving non-linear differential equations with constitutive relationships given by a DNN with known weights.  
  
  
  
  
  
  § MONODISPERSE POLYMER MELTS
   
  
  
  We consider a synthetic Dissipative Particle Dynamic (DPD) fluid consisting of chains of N equal-size beads connected by springs to model polymer melts. Two-dimensional DPD simulations of such fluids between two parallel plates with chains made of N=2, 5, and 25 beads are presented in  <cit.>. In <cit.>, the DPD results were used to compute  μ_DPD(du(y)/dy) using the Irving-Kirkwood relationship <cit.>. 
  
   We use the velocity data from <cit.> and the inverse PINN method  with ω_1 = ω_2= ω_3 = 1 and ω_4 = 0 in Eq. (<ref>) to estimate μ(u_y). To match <cit.> , C = 0.75 and H=25. The relative velocity error and the maximum residual error are given in  Table <ref>. For all considered N the relative error in u is less than 0.1% and the maximum residual error is 3 orders of magnitude smaller than the driving force C, indicating that the DNN û accurately approximates data and the DNNs û and μ̂ satisfy the governing equations. 
  
  
  
  
  Figs. <ref>a and <ref>b compare the velocity profiles and viscosities estimated from the DPD simulation, μ_DPD and from the PINN method for N=2. The DNN velocity profile û(y;θ) closely matches the DPD velocity profile u_DPD(y). The agreement between μ̂(u_y;θ,γ) and μ_DPD(u_y)  is good but less accurate than the agreement for the velocities. To test whether u_DPD(y) and μ_DPD(u_y) satisfy Eq. (<ref>), we train the û(y;θ) and  μ̂(u_y;θ,γ) DNNs conditioned on both u and μ DPD measurements.
  Figs. <ref>a and <ref>b show that  conditioning of the DNNs on the DPD measurements of u and the estimates of μ produces DNNs that match well both u and μ
  data. However, conditioning on the DPD μ estimates also results in the residual errors that are two orders of magnitude larger than the residual errors in the case where no μ estimates are used to train the DNNs, as shown in Fig. <ref>c. 
  
  Next, we use the PINN method to evaluate the viscosity of the polymer melt with 25-bead chains. As for the melt with N=2, we  first train the û and μ̂ DNNs using only u_DPD(y) measurements.  Fig. <ref> shows that the û DNN agrees well with the u_DPD(y) measurements and the resulting residual point errors are nearly zero (more than four orders of magnitude smaller than C). We also see that μ̂ significantly deviates from the μ_DPD(u_y) values estimated from the DPD simulations near a shear rate of zero. Then, we train the û and μ̂ DNNs using both u_DPD(y) and μ_DPD(y) data. Fig. <ref>  demonstrates that the resulting DNNs fit the  u_DPD(y) and μ_DPD(y) data well, but the corresponding residual is very large (on the order of C). We obtain similar results for the polymer melt with N=5. 
  
  Finally, we demonstrate that once μ̂(u_y;θ,γ) is trained, the forward PINN method can be used to solve Eq. (<ref>) subject to the BC Eq. (<ref>). We use the weights γ in the μ̂(u_y;θ,γ) DNN obtained above from the inverse PINN and train the forward solution, û_f(y;θ),  DNN by minimizing the loss function Eq. (<ref>)  with ω_1 = ω_4 = 0 and ω_2 = ω_3 = 1 for C = 0.75. For N=2, Fig. <ref>a shows that the û_f(y;θ) DNN matches the experimental data corresponding to  C = 0.75 well. In addition to this, Fig. <ref>c demonstrates that the residual of the governing equation is two orders of  magnitude smaller than C confirming that û_f(y;θ) approximately solves Eq. (<ref>) subject to Eq. (<ref>).
  
  These results lead to the conclusion that the inverse PINN is capable of estimating the effective viscosity function μ(u_y), which can be used for solving the momentum conservation equation (<ref>). The predicted viscosity deviates from the viscosity obtained from the DPD simulations for small shear rates with the discrepancy increasing with the number of beads N. Our results show that velocity and viscosity data provided in  <cit.> cannot accurately be described by Eq. (<ref>).
  
   
  
  
   
  
  
  
  
  
  § DENSE SUSPENSIONS OF SPHERICAL PARTICLES
  
  
  In this section, we employ the inverse PINN method to learn the shear-rate-dependent viscosity of  densely packed spherical particles suspended in a Newtonian fluid using the velocity measurements presented in  <cit.>. 
  The considered data are obtained from the numerical simulations of suspension flows in a channel using the Force Coupling Method (FCM) <cit.>. In the considered suspensions,  the average particle volume fraction ϕ_a = 4/3π a^3 N/V ranges from  0.2 to 0.4, where a is the particle radius, N is the number of particles, and V is the volume of the domain. 
  In the FCM simulations, the particle radius was set to a = 1, the channel length to L_x = 80, the height to H = 40, and the width to L_z = 30. The channel walls were located at y=0 and 40, constant Dirichlet BCs for pressure were prescribed at the x=0 and 80 boundaries with the pressure drop over the length of the channel Δ P/L_x=0.029, and periodic conditions were used in the z direction. At the continuum level, the considered suspension behaves as a non-Newtonian fluid and can be described by Eq.  (<ref>) with C = Δ P/L_x.
  
  
  The velocity profiles for the suspension flows with ϕ_a=0.2 and 0.4 are shown in Figs. <ref>a and d, respectively, and the local volume fractions ϕ(y) are depicted in  <ref>a.
  A key feature of suspensions is irreversible shear-induced migration of particles to areas of low shear rate <cit.>. 
  
   Particles in a suspension subjected to a Poiseuille flow will migrate to the channel centerline, greatly increasing the volume fraction at the centerline until it reaches the maximum close-packing limit, as shown in Fig. <ref>a. This migration also impacts the velocity profile, resulting in a flattened parabola shape that is observed in Figs. <ref>a and d. 
  
  
  
   
  As in the analysis of polymer melts above,  we use the inverse PINN to find the viscosity μ(u_y(y)) 
   by approximating u and μ with û(y;θ) and μ̂(û_y(y;θ);γ) DNNs trained by minimizing the loss function  (<ref>) with ω_1 = ω_2 =ω_3 =1 and ω_4=0.
  We use N_u=401 measurements of the velocity profile u(y) from the FCM simulations. Because the velocity profiles from the simulations (see Fig. <ref>a and d) deviate from the flattened-parabola shape near the walls due to particle layering, a phenomenon that cannot be described by  Eq.  (<ref>), we train the PINN with velocity data in the range y ∈ [0.25y/h, 1.75y/h ], but still impose the zero Dirichlet BCs for u at y = 0, H. 
  
   Figs. <ref>a and  <ref>d compare the velocity profiles of the suspension flow observed in the numerical simulations and are approximated with the û(y;θ) DNN for ϕ_a=0.2 and 0.4, respectively.  The μ̂(û_y(y;θ);γ) DNN and the viscosity estimated from the numerical experiments are plotted in Figs. <ref>b and  <ref>e. 
   The viscosity μ(u_y) for the FCM simulations is found by computing  u_y(y) and ϕ(y) from the simulation data, assuming that μ(u_y) = η_s(ϕ(u_y)) η_f and using the Eilers formula
  η_s(ϕ) = (1 + 5ϕ/4(1-ϕ/ϕ_c))^2  <cit.>. Here, η_f is the fluid viscosity (which was set to unity in the FCM simulations) and ϕ_c is the maximum volume fraction of a suspension (ϕ_c = 0.62 in the FCM simulations.)
  We observe that the PINN method is able to accurately learn the velocity profile and captures the increase in viscosity at the channel centerline. Figs. <ref>c and <ref>f demonstrate that  the residuals are three orders smaller than C = 0.0288, indicating that the  û(y;θ) and μ̂(û_y(y;θ);γ) DNNs  satisfy Eq. (<ref>).   
   
   
  
  Finally, we employ the inverse PINN method to evaluate η_s as a function of ϕ and in Fig. <ref>b compare it with the  Eilers model and the Krieger  model 
  	η_s(ϕ) = (1-ϕ/ϕ_c)^-2.5ϕ_c <cit.>.
  In the PINN method, we compute η_s(ϕ) using the μ̂(û_y(y;θ);γ) and û(y;θ) DNN models of viscosity and velocity and ϕ(y) observed in the FCM simulations. The considered empirical models predict similar μ values for ϕ<0.35 away from the channel centerline. The empirical models assume that μ(ϕ) is independent of ϕ_a. Fig. <ref>b shows that the PINN predicted μ(ϕ) functions agree with the empirical models  for small ϕ for all considered ϕ_a. For large ϕ, the PINN estimated μ(ϕ) relationships depend on ϕ_a and deviate from all considered empirical models. There are several reasons for the disagreement between the PINN estimated and empirical viscosity models.  At high volume fractions near the close-packing limit, ϕ_c, the particle movements are highly correlated leading to non-locality of the particle forces. Therefore, Eq.  (<ref>) breaks down for dense suspensions at the centerline. Additionally, Eq. (<ref>) with a constitutive relationship of the Eilers or Krieger analytical forms predicts that the suspension will reach maximum packing with ϕ = ϕ_c at the centerline  that is independent of the average volume fraction ϕ_a <cit.>. However, FCM simulations <cit.> and experiments <cit.> show that the volume at the centerline varies with the initial average volume fraction of the system.  
  
  
   
  
  
  In conclusion, we have extended the PINN method for learning unknown physics, including the functional dependence of viscosity on the shear rate and other properties of fluids using indirect measurements such as fluid velocity and volume fraction. We have also demonstrated that once an accurate DNN approximation of the viscosity is available, the PINN method can be used to model non-Newtonian flow without any data except the boundary conditions. 
  
  
  
  
  
  
  
  
  
  
  