Performance prediction of magnetorheological fluid‐based liquid gating membrane by kriging machine learning method

A smart liquid gating membrane is a responsive structural material as a pressure‐driven system that consists of solid membrane and dynamic liquid, responding to the external field. An accurate prediction of rheological and mechanical properties is important for the designs of liquid gating membranes for various applications. However, high predicted accuracy by the traditional sequential method requires a large amount of experimental data, which is not practical in some situations. To conquer these problems, artificial intelligence has promoted the rapid development of material science in recent years, bringing hope to solve these challenges. Here we propose a Kriging machine learning model with an active candidate region, which can be smartly updated by an expected improvement probability method to increase the local accuracy near the most sensitive search region, to predict the mechanical and rheological performance of liquid gating system with an active minimal size of experimental data. Besides this, this new machine learning model can instruct our experiments with optimal size. The methods are then verified by liquid gating membrane with magnetorheological fluids, which would be of wide interest for the design of potential liquid gating applications in drug release, microfluidic logic, dynamic fluid control, and beyond.


| INTRODUCTION
As the demand for functionality and adaptivity of materials grows, exploiting adaptive materials responding to external stimuli has come to the fore with significant advantages. [1][2][3][4][5] Recently, liquid gating system was put forward as an emerging technology, which has given great inspiration to design adaptive smart materials with specific functions for various applications. [6,7] In such a system, the solid matrix material provides stable mechanical structure and the capillary-stabilized liquid material offers dynamic functional adaptivity. [8] Once the substrate material is set, liquid materials are key factors in the property and function of the liquid gating system. Commonly employed liquids are perfluorinated liquids, silicone oils, liquid metal, ionic liquids, and many responsive liquids. [9,10] In particular, with the advantage of remote and noncontact property, magnetorheological fluid-based liquid gating membrane (MRF-LGM), in which magnetorheological fluid (MRF) as a magnetically-responsive functional liquid in liquid gating membrane, have great potentials for the application of microfluidic control, drug release, and chemical reaction. [11] In order to accurately describe the dynamic rheology of MRF, current physical models have been developed with additional physical parameters, which greatly affect the pressure threshold of liquid gating membrane from a theoretical point of view, [12] and therefore influence its potential application. In addition, as a typical non-Newtonian fluid, a small value vibration of these properties could lead to a big performance change of MRF in liquid gating membrane, which means a high predictive accuracy of MRF is also required. [13] Traditional works have been made for accurately predicting the important properties based on the limited experimental data from the rheological tests. [14,15] Esmaeilnezhad et al. [16] used the regression equations, which consist of the linear, two-factor interaction, and quadratic model in Bingham fluid model for predicting and response surface methodology (RSM) for error optimization of MRF properties' prediction, though, only two parameters were involved. Kemerli et al. [17] applied the numerical CFD to map the relationship between the viscosity η and yield stress τ 0 , magnetic flux density B, shear rate γ, that is, established a function of η = f (τ 0 (B), γ). Cvek et al. [13] focused on evaluating the fitting capabilities of three-parameter viscoplastic models, including the Herschel-Bulkley model and Robertson-Stiff model, for the description of MRF. Lv et al. [18] also compared the accuracy of three constitutive models of MRF via simple fitting algorithms. Based on the previous study, fitting the experimental data with a sequential optimization method is widely used for obtaining essential parameters. The theoretical model with more parameters means higher dimensions and brings problems of predictive accuracy, and it will be even worse if the experimental data is insufficient, which limits the further guidance for the experiment and is not practical in some extreme situations. To address these challenges, it requires a new method that can meet these needs.
Machine learning methods, including regression method, neural networks, support vector machine (SVM), clustering algorithm (K-Means), can build the relationship between inputs and outputs based on limited data, which is promising in predicting of material science, Miao et al. [19] accelerated the discovery of new catalysts by combining active machine learning with density functional theory (DFT) calculations, and the prediction accuracy of the new neural network constructed by Jia et al. [20] is close to the calculation results of DFT. In the field of MRF material predicting, Bahiuddin et al. [21,22] developed an extreme learning machine method to predict the properties of the MRF constitutive model. They used the machine learning model and minimal norm least square optimization method sequentially, which is complex and inefficient. As a typical Gaussian regression learning method, the kriging machine learning method is widely used in fields of aircraft, [23,24] auto vehicle designs [25,26] and quantum chemical computations, [27][28][29] and so forth. However, because the size and select strategy of training samples could not be obtained before calculation, machine learning methods generally need to learn from the total database and the efficiency of the traditional kriging learning method is therefore poor. With the use of stochastic processes, Bayesian optimization (BO) is suitable for solving global optimization problems with objective functions that are unknown or implicit to designers, and computationally or experimentally expensive. [30] BO can be capable of combining the kriging machine learning to using a stochastic process to develop "figures of merit" for smartly selecting the search points, which has been reported as an efficient global optimization (EGO). [31] Efforts were made to reduce the size of the candidate region of EGO, which is another direction for improvin. Wen et al. proposed an adaptive sampling region method to avoid selecting the additional points where the probability density is low. [32] Yang et al. developed a strategy to determine the minimal or optimal size of candidate points [33] for reliability assessment.
In this paper, based on the BO and kriging machine learning method, we develop a new active kriging machine learning method (KAML) modifying with an active candidate region (ACR) technique to predict efficiently the properties of MRF-LGM with a high accuracy based on our experimental data. The efficiency of the kriging machine active learning method with the active candidate region techniques (KMAL-ACR) method is demonstrated by the famous mathematical problems. The results are compared with some current machine learning methods and different optimization techniques. Finally, we also give a strategy to obtain the most valuable tests for directing the experiments to decrease the experimental workloads for MRF-LGM.

| MAGNETORHEOLOGICAL FLUID-BASED LIQUID GATING MEMBRANE
The MRF-LGM is formed by impregnating the magnetically-responsive smart MRF into a porous solid membrane framework. The MRF suspension with microscale magnetic colloids is capillary-stabilized in the confined space. Without the magnetic field, magnetic colloids are randomly distributed in confined pores. When the immiscible fluid transports through the MRF-LGM, a certain pressure threshold is required to conquer the capillary force. Once the magnetic field is applied, which is perpendicular to the transport flow, magnetic colloids will form a chain structure along the direction of the magnetic field, so it needs higher pressure threshold when the transport fluid permeates the MRF-LGM due to the stable and robust chain structure ( Figure 1A). Our previous research has revealed that the pressure threshold is an intuitive indicator of the colloidal interaction strength in the MRF-LGM, which can provide a direct strategy for characterizing the colloids' mechanical property quantitatively. [11] For a new liquid gating membrane with unknown characteristics, acquiring some key factors including the pressure threshold of the system, as well as the rheological, mechanical property of functional gating liquid is very critical. However, some experiments cannot be completed especially in the conditions such as ultrahigh pressure or temperature, so they need a large number of basic experiments and theoretical computations with a lot of time to conduct theoretical behavior prediction. Therefore, the machine learning method with adaptive learning technique has a great capacity for solving nonlinear and implicit problems with higher accuracy and efficiency, which is shown in Figure 1B where μ is viscosity, k is the permeability of the medium, Q is the flow rate, h is the length of pressure drop and A is the cross-sectional area to flow. L is the effective shear length. The pressure threshold is therefore, The yield stress τ y depends on the magnetic strength, which can be measured separately with rheometers. Also, τ y can be estimated theoretically by the micro chain-rod model [34] : where a is the radius of magnetic particle, δ is the distance between each particle, χ is magnetic susceptibility, φ is the fraction of porosity, H is applied magnetic field. Or it can be expressed by the rheological model, and some of which are shown in Table 1, including the Bingham plastic model, [16,34] Bingham biplastic (BB) model, Herschel-Bulkley (HB) model, Mizrahi-Berk (MB) model, [13] Papanastasiou PP model, [12] and Robertson-Stiff (RS) model, [36] where η, γ, τ, and τ 0 are the viscosity, shear rate, shear stress, and yield stress of MRF, τ y1 , and τ y2 are the dynamic and apparent yield stress, η 1 and η 2 are the viscosities of the low and high shear rate regions, K is the denseness coefficient, and n p is the rheology index, respectively. However, K and n p have no definite physical meaning. Our studies mainly would focus on the improvement achieved by the PP model and the BB model.

| Kriging active machine learning method with adaptive candidate region
Kriging machine learning method as a typical Bayesian algorithm, active Kriging learning method treats uncertainty x as the normally distributed random variables. [31] The active learning function or acquisition function is used to update the new training point. For example, the expected improvement (EI) function is defined as a normal distribution.
where ϕ is the standard normal cumulative function, and  ψ ( ) is the standard normal probability density function. The point with a maximal EI deserves to be learned in the next learning step. We have proposed a new acquisition function, named EI with exponent acceleration (EIE), to avoid selecting the boundary points even when they are far from the optimum point, [37] T A B L E 1 General rheological models of MRF Rheological model

Number of parameters
Effective parameters

Features
Bingham plastic [16,34] τ τ ηγ = +0 2 2 Simple, widely used, but may not suitable for magnetorheological fluid (MRF) Two lines, cannot fit data accurately Herschel-Bulkley (HB) [13] τ τ Low efficiency but more accurate than Bingham model Papanastasiou PP [12]  Good performance but with a higher computational cost Robertson-Stiff (RS) [35]  Better performance in a steady shear region x min is the minimal point, n is the number of elements in the total candidates, which is often larger than the number of original sample points. Its convergence condition is . However, this leads to a higher probability to fall into a local optimum.
Thus, an improvement exponent distance (EID) functionthat adds a new item representing the error of selecting points too far from the current predicted optimum point: where p is a positive integer indicating the number of iterations. Convergence condition is If we want to obtain a global optimum with the error of ε, N is required above 10 2−log(ε) . For example, ε = 0.01 will need 10 4 samples. Two problems are brought: (1) a large number of candidate samples should be estimated at each iteration; (2) The accuracy will reach a limit threshold value if the real optimum is not covered by the initial large samples.
Therefore, we develop a method using an active candidate region (ACR) technique to solve these problems based on the theory of applied statistics and probability. [32,33] At one iterative additional training point, the probability when the real optimal point locates in the current ACR is is the probable sample point, n is the dimension of variables, P ACR is the required value of probability. f(x) is the joint probability density function. We can assume that the i th ACR is a region of normal distribution N μ σ ( , ) Transfer it as a standard normal distribution, we have If we treat the ACR as a high-dimensional ellipsoid, then the β i , which relates to the accuracy of the KAML model, will decrease as continuously learning from the additional training points. Therefore, the area of ACR will shrink and the accuracy of the optimal point will increase.

| Rheological model based on KAML-ACR
Here, the Papanastasiou PP model with a 5-dimensional prediction and optimization problem is demonstrated by our KAML-ACR model as an example. The traditional learning & optimization method divides into two steps, which burden the computation resources. Our KAML-ACR method itself is an optimization algorithm, which can directly and efficiently accomplish the field-independent rheological modeling and optimization of MRF ( Figure 2).
The strategy of KAML-ACR for MRF is therefore proposed and the steps are as follows.
Step ( , where i is the number of objective functions. Step (2)  . Then, we can obtain their mean value and standard derivation.
Step (3) The EI or EID value of each point and the standard deviation of the improvement distribution can be then calculated based on Equations (5) and (6) . Step , renew the candidate region by calculating β i min ( ) . Then active changed candidate points are sampled by Monte Carlo simulation (Random sampling) following these normal distributions and the sample points in the region out of [−3σ a , 3σ a ] should be deleted. Then, go to Step 2 using the new candidate region. If not, Optmodel i (P) is the final KAML model and the optimal predictive parameters are obtained.
Pseudo codes of KAML-ACR can be found in Online Supporting Information Text 1.5.

| Tests on KAML-ACR
The KAML-ACR technique, along with other prevalent machine learning algorithms, such as radial basis function neural networks (RBFNN), [38] original Kriging machine learning (KML), Kriging active learning method (KAML), is compared in performances by the test functions, [39] which are usually used to test the machine learning algorithms as benchmarks in solving problems with the multimodal, small-gradient, and extreme fea- process by the KAML-ACR method can be found in Figures S1-S5. It should be noted that the methods of RBFNN and KML essentially are only predictive modeling without an ability of optimal searching, genetic algorithm (GA) is used additionally to find the final optimal results. The results are shown in Table 2. And Table S1 compared the results of Goldstein-Price's function (GPF) by various methods. The increasing of dimensions would bring challenges to all learning methods, however, as a comparison, the KAML-ACR method still obtains a more accurate prediction in this situation. KAML-ACR method also has advantages in optimization compared with other prevalent optimization algorithms (Online Supporting Information text 1.2) based on the results shown in Table S2.

| EXPERIMENTAL METHODS
In this study, the MRF with a carbonyl iron particles weight percentage of 75% was purchased from Henan Huya Trading Co., Ltd. It was centrifuged at 11 000 r/min to obtain the carbonyl iron particles. The microscopic view of carbonyl iron particles was sputter-coated with a platinum layer for 30 s and observed by an Olympus IX73 Microscope. A thin layer of the MRF (75 wt.%) was coated on a glass slide for the observation of the chain structure with and without applying the magnetic field. The carbonyl iron particles were suspended in deionized water and then the particle size of carbonyl iron particles was analyzed by an LS-603 Laser Particle Size Analyzer. The viscosity and shear stress of MRF at various shear rates under the magnetic field is determined by the Anton Parr MCR 302 Rheometer with a self-designed electromagnetic coil. The shear rate range was from 0.1 to 1001/s, where the data point was recorded every 3 s. The magnetization versus magnetic field (B-H loop) curve of carbonyl iron particles is determined by a LakeShore 7404 Vibrating Sample Magnetometer with the magnetic field from −20 000 Gauss to 20 000 Gauss. When the external field is applied, the particles would move and finally become into several chain-like aggregates along with the field direction. That is the source of magnetic viscosity increasing ( Figure S8).

| Proposed optimization problems of MRF
Both the magnetic field and measured shear rates are used as inputs variables, while the shear stress predicted  | 163 from the KAML-ACR training is used as an output variable. Papanastasiou PP model including five parameters is used here as the constraints for having the highest accuracy and flexibility. [12] The root means square error (RMSE), which is suitable to describe the error of nonlinear model fitting, is used as the objective functions of the optimization problem. Generally, the apparent yield stress τ y1 is larger than the dynamic yield stress τ y2 , but they are not able to be larger than the maximum stress obtained by experiments τ exp . And apparent viscosity η 1 is larger than dynamic viscosity η 2 , and the range of shear rate is from 0 to 100 s −1 , which can be also set as the range of γ c in optimization.
Reference [21] has tested several nonlinear piecewise continuous functions, including sigmoid, hardlimit, and sine functions, as the hidden-node function for the machine learning model for estimating the predictive machine learning model. In this case, because Bezier spline regression could pass every data point, it is used to calculate the x γ τ ( ,) i with our KAML model as a baseline.

| Results and discussion
The plots of RMSE and shear stress v.s. the shear rate of different rheological models established by using KAML-ACR of MRF when B = 600 Gs (n = 1000) are shown in Figure 3A,B. With using ACR techniques, the RMSE decreases sharply after learning about 10 additional samples and keeps stable after learning 100 samples ( Figure 3A). As five predicted parameters are involved, the minimal RMSEs of the BB and PP model are smaller than that of other models (details can be found in Table S3). Table S3 shows the details of parameters predicted by different models in Figure 3B and the PP model has the smallest RMSE (9.9341 Pa) among all. BB and RS models also have a better performance. We notice that the predictive value of parameters K and n p are much different in HB, MB, and RS models, which indicates that K and n p have only the mathematical meaning rather than the definite physical meaning. In conclusion, the PP model performs better in accuracy than other models. Figure 3C shows that the RMSEs of all methods can be stable after about 20 training points even with a size difference from 50 to 1000. Table S4 gives the details of predictive parameters based on different n, which shows that a larger size can help to get a more accurate predictive value. Figure 3D and Table S5 show that KAML with ACR will have a good potential application in global optimizations. Besides this, the accuracy and convergence speed of results by KAML using the ACR technique can be both better than that without ACR (see Figure S7).
Then, according to experimental data, we discuss the prediction of the rheological parameters in different field strengths based on the PP model ( Figure 3E). It shows that the averaged shear stress will increase as the B increases. Examples in reference [12] show that minimal RMSE optimized by the PSO algorithm could be as high as 100 Pa if averaged shear stresses obtained from experiments are above 5 kPa. The comparison between the predictive shear stress and our experimental data is shown in Figure 3F. The discussion on the detailed rheological parameters is based on Table S6. We can find that our method could be suitable for a wide range of magnetic fields, even with limited experimental data.

| Prediction with field-dependent parameter
According to the previous study, yield stress has a directive relationship with the magnetic field [17] or a quadratic polynomial τ y = AB 2 + τ * for details, [40] where A is the unknown linear constant, τ * is the yield stress without the external field. The mathematical model can be written as, The RMSE change in history is shown in Figure 4A. This time, we set a convergence condition ε ≤ 10 −5 and optimal results can be obtained after only 44 additional learning points. We also give a 2D plot to show its searching process in Figure 4B. The ACR could be easily found to shrink and focus on the region near the global optimum, thus, additional points to train the KAML also locate and focus on this area.
To verify our model of prediction, we also obtain a group of experiment data at B = 1600 Gs, which is used as a benchmark for prediction. The predictive dynamic and apparent yield stress are 3524.3 Pa and 5587.1 Pa. Compared to the direct method, errors are small enough though only seven effective learning points, shown in Figure 4C,D. Therefore, the KAML-ACR model can also instruct the experimental design. If we treat the magnetic field induction B as an unknown parameter and find the most important value B that influences the predictive error, the experimental workloads could be decreased. The mathematical model could be revised based on this point as follows. and testing their shear stresses. After the prediction by Equation (13), the RMSE of the KAML-ACR model could decrease from 163.4842 Pa to 9.9915 Pa with only 10 sample points in total (original three samples adding with 7 additional learning samples) as shown in the inserted figure in Figure 4E, which is near the minimal RMSE (9.9341 Pa) with 40 samples in Section 5.2. That is to say, the experimental workloads are only about 1/6 of that in previous works. Also, it shows its advantages on the prediction of field-dependent variables by using the dynamic yield stress prediction as an example. For example, only two original experimental data at B = 0 and 600 Gs have been obtained, shown in Figure 4F. Conducting the machine learning, the most valuable experimental test would be approached at B = 964.21 Gs. And if the fourth training point is required to add, repeat this step and a new sample B = 1097.93 Gs could be found ( Figure 4F). The RMSE could decrease from 1.3173 Pa to 0.3729 Pa, which is near the minimal RMSE 0.2499 Pa with seven samples in Figure 4C. Therefore, the methodology obtaining the rheological properties of MRF can save experimental workloads instead of the theoretical calculation with a much lower economic cost, which would be of significance, especially for some expensive material testing.
F I G U R E 4 Prediction and strategy with field-dependent parameter by KAML-ACR method. (A) RMSE history in predicting with fielddependent parameters. (B) 2D plot of the optimal searching history, where A 1 and A 2 are the fitted parameters in Equation (13). (C) Results in predicting with field-dependent parameters for dynamic yield stress. (D) Results in predicting with field-dependent parameters for apparent yield stress. (E) The most valuable test points selected by KAML-ACR method (Shear stress vs. shear rate). (F) The most valuable test points selected by KAML-ACR method (Pressure threshold vs. magnetic field induction B). KAML-ACR, Kriging machine learning method-active candidate region; RMSE, root means square error 6 | CONCLUSIONS In order To accurately predict perform of MRF-LGM, we introduce a new kriging machine learning method that can accurately obtain the rheological properties of MRF. The method is essentially an active learning prediction of the rheological model based on the Papanastasiou PP model obtained by using KMAL-ACR. With continually adding new training points, it can increase the accuracy near the global optimum by active learning. Meanwhile, the accuracy can be further increased by decreasing the area of the search region, named the ACR method. Some conclusions are summarized as follows.
1) KAML with ACR methods only has one process to predict the rheological properties of MRF. Compared with traditional machine learning and optimization algorithms, ALK-ACR is more efficient. Besides this, the decreasing of the standard deviation of the learning function could continuously decrease the area of the candidate region, which therefore increases the accuracy of optimization. 2) In terms of global optimization, KAML-ACR performs better than some current optimization algorithms, such as particle swarm optimization algorithm (PSO), genetic algorithm (GA), and adaptive simulated annealing algorithm (ASA), especially in the performance of convergence speed. And it also has relatively higher accuracy in solving high-dimensional optimization problems. 3) By using our KAML-ACR method, we compared the fitting performance of five rheological models for describing the MRF based on our experimental data, including the Bingham biplastic (BB) model, The KAML-ACR is also very convenient for predicting accurately the field-dependent parameters and helps researchers with scientific selection strategies for the most appropriate yield stress of the smart material characterization and design of magnetorheological fluids' devices or systems. The most impressive advantage is that it can actively select the experimental sample with the help of the ACR technique, which can save lots of experimental workloads instead and would be especially suitable for some expensive material tests.
Finally, the proposed methodology, including the prediction of the yield stress, viscosity, and directions for experiments, can be used for the MRF-LGM based on limited experimental data and will accelerate the design of other potential applications in drug release, microfluidic logic, dynamic fluid control and beyond.

SUPPORTING INFORMATION
Additional supporting information may be found in the online version of the article at the publisher's website.