Differentiable Fluid Physics Parameter Identification Via Stirring

Humans interact with fluid every day, from mixing up coffee powder into water to squeezing toothpaste. Understanding fluid is the key to state estimation in a larger range of fluid-related manipulation tasks. In domestic scenarios, density and viscosity are two important physics properties to identify. For example, in the kitchen, adjusting density usually relates to adjusting the extent of a certain flavor (e.g., sweetness and saltness). Viscosity can indicate the solidity of the pastry after the bakery. Other physics properties such as elasticity also contribute to the fluid dynamics behaviors but do not have a significant influence on household tasks, thus they are not in our major focus. Specifically, common fluid has a bulk modulus to characterize the elasticity which ranges from 0.9GPa to 4.5GPa [EngineeringToolboxBulks]. It means at least $9$MPa should be applied to the fluid to make it compress by 1% volume, such pressure is unusual in household tasks.

Though density can be easily estimated with Archimedes’ principle, viscosity is more complex. In hydraulic mechanics, the dynamics behaviors of fluids are described by different constitutive models [harris2005fast, anderson1995computational, park2023fluid], which means they can be characterized by a few parameters given the predefined stress-strain relationship. Based on the stress-strain relationship, fluids can be categorized in Newtonian (Fig. 1-a) and non-Newtonian (Fig. 1-b) types. Newtonian fluids are characterized by a linear stress-strain relationship. While non-Newtonian fluids have varied viscosity based on the rate of deformation caused by the shear stress. These fluids are typically described by constitutive models such as the Carreau [carreau1972rheological], Cross [cross1965rheology] and Herschel-Bulkley [herschel1926konsistenzmessungen] models. In these models, viscosity is one of the major parameters.

In this work, we present a novel differentiable fitting framework named DiffStir to identify the key physics parameters through a simple action, stirring, a common operation in the kitchen. It utilizes a robot arm to grasp a rod and conducts the robotic stirring operation on a cup of certain fluid with recorded trajectories and forces. After synchronizing the action and force between the real-world configuration and a differentiable Material Point Method (diffMPM) [hu2019difftaichi]-based simulator, we can estimate the fluid parameters of the pre-set constitutive models by matching the observations from both worlds. Different non-Newtonian constitutive models have their own preferences for fluid. For example, the Carreau model is suitable for shear-thinning fluids that can approach Newtonian behaviors at high shear rates; the Cross model is good for shear-thinning and shear-thickening fluids across a wide shear rate range; and the Herschel-Bulkley model is best for fluids with a yield stress, such as muds and some toothpaste. Based on it, we design an online strategy to adaptively select the most proper constitutive models given the real-world recordings. Moreover, though MPM is a powerful computational tool often employed in fluid dynamics due to its ability to handle large deformations and topological changes in a numerically stable manner, the physical sim-to-real gap is inevitable, not to mention the real-world sensor noise, which causes inaccuracy of the parameter estimation. To address it, we propose a refining neural network to calibrate the parameters between the real world and the simulator.

To evaluate the proposed DiffStir, we conduct extensive experiments from four aspects to demonstrate the accuracy of the estimated parameters with reference values from the literature [EngineeringToolboxDynamicViscosities, EngineeringToolboxViscosity, EngineeringToolboxKinematicViscosities]: density validation, density & apparent viscosity estimation, fluid dynamics behavior classification, and fluid mixing analysis.

We summarize our contribution as:

• •

  We propose to identify the key physics parameters, density and viscosity, of fluid in household scenarios through a common operation, stirring. The identification process bridges the real world and the differentiable simulator, and can achieve fluid understanding during robotic operation.

• •

  We conduct extensive experiments on different Newtonian and non-Newtonian fluids to validate the usability and accuracy of our system. We hope our system can serve as a foundational tool for liquid-related manipulation tasks.

Fluid dynamics exhibit intricate behaviors with many unresolved issues [davidson2015turbulence, majda2002vorticity]. Different constitutive models are proposed to characterize the dynamics behavior of fluids [chen2021constitutive]. In this work, we exploit the Kelvin-Voigt model [dowling1999mechanical].

1) Kelvin-Voigt Fluid Model: The Kelvin-Voigt model, merging a spring (representing elasticity) and a dashpot (representing viscosity) in parallel (Fig. 2(a)), provides the rheological characteristics of diverse viscoelastic materials, such as polymeric solutions, foams, gels, and suspensions [meyers2008mechanical]. The model derives the total stress, $\sigma$, from the collective contributions of both the viscous and the elastic elements:

$$\sigma=E\varepsilon+\eta\dot{\varepsilon},$$

(1)

where $\varepsilon$ is the strain applied to the fluid, $E$ is the elastic modulus to characterize the stress-strain relationship of the spring, $\dot{\varepsilon}$ is the strain rate and $\eta$ is the viscosity parameter.

2) Integrating Volume-Related Energy Density: The Kelvin-Voigt model utilizes a bulk modulus $\kappa$ to quantify the elastic response to volume change. Since the fluid undergoes non-linear elastic deformation, we employ the Neo-Hookean model, which details the non-linear elastic reaction to volume changes through the energy density function $W_{v}$, related to the volume ratio $J$:

$$W_{v}(J)=\frac{\kappa}{2}\left(\frac{1}{2}(J^{2}-1)-\ln J\right)$$

(2)

Differentiating $W_{v}$, we obtain stress:

$$\tau_{v}=\frac{\partial W_{v}(J)}{\partial J}=\frac{\kappa}{2}(J-\frac{1}{J})$$

(3)

Incorporating this stress into the Kelvin-Voigt model:

$$\sigma=\tau_{v}I+\eta\dot{\epsilon}$$

(4)

Here, $\tau_{v}$ multiplies the identity tensor $I$, reflecting its isotropic effect.

As mentioned earlier, common fluid is hard to undergo significant volume change with the external force in household scenarios, in simulation, the volume-related stress is applied to keep the volume hold. Without this term, the volume will shrink in simulation.

3) Extending to Non-Newtonian Fluids: To fully grasp how fluids behave, we combine models from continuum mechanics with non-Newtonian fluid models. Using this approach, we focus on three key non-Newtonian models: Carreau, Cross, and Herschel-Bulkley [escudier2002fully]. These models help explain the changing relationship between viscosity $\eta$ and shear rate $\dot{\gamma}$ instead of treating viscosity as a fixed value. The function is named apparent viscosity and denoted as $\eta(\dot{\gamma})$. The detailed formulation of these models can be referred to in [escudier2002fully] and our supplementary materials.

As shown in Fig. 3, the parameters of these models, denoted by $\theta$ are: (a) $(\eta_{0},\eta_{\infty},\lambda,n)$ for the Carreau Model; (b) $(\eta_{0},\eta_{\infty},\lambda,n)$ for the Cross Model; (c) $(\tau_{y},n,k)$ for the Herschel-Bulkley Model, where $\eta_{0}$ is the zero shear rate viscosity, $\eta_{\infty}$ is the infinite shear rate viscosity, $\lambda$ is a multiplier constant, $n$ is the power-law index, $\tau_{y}$ is the yield stress and $k$ is a consistency index. To note, the viscosity from the Herschel-Bulkley can be derived from the estimated yield. The formulation can be referred to in our supplementary materials.

We take the Carreau model as an example to demonstrate how the extended non-Newtonian formulation shapes. The apparent viscosity is given by:

$$\eta(\dot{\gamma})=\eta_{\infty}+(\eta_{0}-\eta_{\infty})(1+(\lambda\dot{\gamma})^{2})^{(n-1)/2}$$

(5)

And the fluid’s stress, $\sigma$, is a function of the parameters $\theta=(\eta_{0},\eta_{\infty},\lambda,n)$:

$$\sigma(\theta)=\tau_{v}I+\left(\eta_{\infty}+(\eta_{0}-\eta_{\infty})(1+(\lambda\dot{\gamma})^{2})^{(n-1)/2}\right)\dot{\epsilon}$$

(6)

In this section, we describe the proposed DiffStir system. Density and viscosity are estimated individually and sequentially. The overall pipeline is illustrated in Fig. 3.

1) Density Estimation: By integrating a force sensor at a robotic arm’s end-effector, our method estimates density via force measurement on a submerged stirring rod. The density $\rho$ is calculated with Archimedes’ principle as:

$$\rho=\frac{G-R}{Vg},$$

(7)

where $G$ is the gravitational force of the stirring rod, $R$ is the equilibrium force measured by the force sensor, $V$ is the volume of the immersed part of the rod, and the $g$ is the gravitational acceleration. $G$, $R$, $V$ can all be measured from the real world, the experimental setup of measurement is described in Sec. V-C1.

2) Differentiable Simulation: To further estimate the viscosity, we build a bridge between the real and simulated observation via a differentiable simulator. We first replicate the setup of the robotic stirring in the real world to the simulator. The density estimated from the last step is used to set up the fluid in the simulated world. Then, we conduct the stirring operation in the real world, and transfer the recorded trajectories (e.g. rod positions, velocities) and the force applied to the rod to the simulated environment. The physics engine MLS-MPM [hu2018moving] in the simulator can produce the interaction force between the fluid and the rod, $F_{calculated}$. The fluid is modeled by multiple particles, and the particles are moved by the stirring operation. The particle dynamics is governed by a constitutive model. Thus, the same stirring operation can cause different particle movements due to different constitutive models. If the constitutive model contains the parameter of viscosity, the particle movement can exhibit the effect of drag and horizontal shear. In this sense, the reactive $F_{calculated}$ is parametric. In our simulation, we represent the fluid as a collection of particles. Each particle’s stress, $\sigma$, is governed by a constitutive relation. The force exerted on a single particle due to this stress can be computed as $\mathbf{F}_{\text{particle}}(\theta)=\int\sigma(\theta)\cdot\mathbf{n}\,dA$, where $\mathbf{n}$ is the outward normal of the particle’s surface and $A$ is the area over which the stress acts. The cumulative force exerted on the stirring rod, resulting from all particles in contact with it, is given by $F_{\text{calculated}}(\theta)=\sum_{i=1}^{N}\mathbf{F}_{\text{particle}_{i}}(\theta)$. Here, $N$ means the number of particles interacting with the rod.

3) Parameter Identification and Model Selection: We attach a thin-film pressure sensor to the gripper finger and measure the force from the rod, $F_{measured}$. By assuming no contact between the rod and the liquid container, the force applied to the pressure sensor is the same as the force applied to the rod from the fluid. In this way, we can compare the measured force $F_{measured}$ in the real world, and the calculated force $F_{calculated}$ in the simulator, and construct a loss term to optimize the parameters $\theta$ associated with $F_{calculated}$:

$$\text{Loss}(\theta)=\text{Mean}\left(|F_{\text{measured}}-F_{\text{calculated}}(\theta)|\right).$$

(8)

By leveraging the power of differentiable simulation, optimization of $\theta$ can be simply achieved using the gradient descent algorithm, expressed as:

$$\theta^{(t+1)}=\theta^{(t)}-lr\cdot\nabla\text{Loss}(\theta)$$

(9)

$lr$ denotes the learning rate, $t$ is the iteration, and $\nabla\text{Loss}(\theta)$ is the gradient of the loss function concerning $\theta$.

Different constitutive models prefer specific fluids; thus, our system uses the Carreau, Cross, and Herschel-Bulkley models to suit different fluids. We employ an adaptive method to select the most suitable model and parameters, basing our choice on the lowest error from Eq. 8 after three estimations. Notably, while these models predominantly describe non-Newtonian fluids, they can capture Newtonian behaviors under specific conditions, e.g., the Cross model when the power-law index $n$ nears zero.

4) Viscosity Refinement: Given the numerical differences and measurement noises between simulations and real-world setups, our system uses a multi-layer perceptron (MLP) network to correct these systematic errors. We generate 2500 estimated/reference viscosity pairs as the training samples. We leave the details of the training set construction in supplementary materials.

Our MLP rectifier has an architecture that consists of two input neurons (representing average stirring rate and preliminary viscosity), a hidden layer with 8 neurons, and a single output neuron for refined viscosity. For training, we set a learning rate of 0.01, achieving a balance between speed and stability in convergence. The model underwent 50 epochs of training, allowing for sufficient learning and weight adjustments.

In this work, we introduced DiffStir, a novel differentiable fitting framework, adept at identifying the key physics parameters through a commonly practiced action, stirring. By using a robotic arm and synchronizing real-world actions with a diffMPM-based simulator, the density and viscosity of various fluids are estimated with considerable accuracy. We’ve showcased that not only can the method adapt to different fluid dynamics behaviors but also has the potential to bridge the gap between physical simulations and real-world observations. Through comprehensive evaluations, DiffStir has exhibited its efficiency and precision in parameter estimation, marking a notable advancement in the domain of fluid dynamics understanding in everyday scenarios.