Real-time control of multiphase processes with learned operators

R E A L - T I M E C O N T R O L O F M U L T I P H A S E P R O C E S S E S W I T H L E A R N E D O P E R A T O R S Paolo Guida Clean Energy Research Center (CERP) Physical Science and Engineering (PSE) Di vision King Abdullah Univ ersity of Science and T echnology Thuwal, Saudi Arabia 23955 paolo.guida@kaust.edu.sa Didier Barradas-Bautista Kaust V isualization Lab Core Lab Division King Abdullah Univ ersity of Science and T echnology Thuwal, Saudi Arabia 23955 didier.bautista@kaust.edu.sa March 27, 2026 A B S T R AC T Multiphase flows frequently occur naturally and in manuf actured devices. Controlling such phenom- ena is extremely challenging due to the strongly non-linear dynamics, rapid phase transitions, and the limited spatial and temporal resolution of av ailable sensors, which can lead to significant inaccu- racies in predicting and managing these flows. In most cases, numerical models are the only way to access high spatial and temporal resolution data to an extent that allows for fine control. While em- bedding numerical models in control algorithms could enable fine control of multiphase processes, the significant computational burden currently limits their practical application. This work proposes a surrogate-assisted model predicti ve control (MPC) frame work for regulating multiphase processes using learned operators. A Fourier Neural Operator (FNO) is trained to forecast the spatiotemporal ev olution of a phase-indicator field (the volume fraction) o ver a finite horizon from a short history of recent states and a candidate actuation signal. The neural operator surrogate is then iterati vely called during the optimisation process to identify the optimal control variable. T o illustrate the approach, we solv e an optimal control problem (OCP) on a tw o-phase Eulerian b ubble column. Here, the con- troller tracks piecewise-constant liquid le vel setpoints by adjusting the gas flo w rate introduced into the system. The algorithm maps the predicted field to a relev ant observable, the height of the liquid column in this case, and minimizes a receding-horizon objective. Because the resulting cost can be noncon ve x and nonsmooth due to the level-e xtraction operation, which in volv es ev aluating the vol- ume fraction against a pre-established threshold, we emplo y Bayesian optimisation to select the inlet velocity that minimises the finite-horizon tracking error with few surrogate rollouts. The results we obtained indicate that field-lev el forecasting with FNOs are well suited for closed-loop optimization since they have relatively low ev aluation cost. The latter provide a practical route tow ard MPC for fast multiphase unit operations and a foundation for future extensions to partial observability and physics-informed operator learning. 1 1 Introduction Multiphase flo ws are ubiquitous in engineering and characterise a lar ge fraction of unit operations, biological systems, and emerging manufacturing technologies [ 1 , 2 , 3 , 4 , 5 ]. It is common to for processes to inv olve immiscible liquid- liquid dispersions in extraction equipments [ 6 ] where they govern droplet breakup/coalescence and interfacial mass transfer . In bio-processing and environmental engineering, multiphase reactors, such as bubble columns, are widely used due to their fa vourable hydrodynamics and enhanced mass transfer [ 7 ]. Other gas-liquid systems include, for example, syngas production reactors via fermentation [ 8 ] and spray cooling [ 9 ]. Beyond “canonical” process equip- ment, multiphase microfluidics is at play in inkjet-based additi ve manufacturing via jetting, breakup, wetting, and ev aporation physics [ 4 ], while metal additive manufacturing in volves highly transient vapour -liquid-solid interactions, 1 The code is av ailable at https://gitlab.com/paolo.guida/neuralmodelpredictivecontrol.git A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 such as keyhole instability and pore formation [ 10 ]. Finally , multiphase phenomena also arise in electrochemical en- ergy systems, where gas generation in lithium-ion batteries has significant safety implications [ 11 ], and even in the human body , where liquid-liquid separation occurs in cells, effecti vely creating a multiphase system [ 12 ]. Many of the multiphase processes discussed above ev olve ov er very short characteristic time scales, leaving little room for direct measurements, calculations and actuation, therefore limiting fine control in practice. Most control strategies for mul- tiphase processes, therefore, do not inv olve feedback and are generally passiv e approaches. The latter , although often effecti ve, may struggle when, for e xample, feed v ariability is significant. It is also important to note that measurements are often not representativ e of the entire time ev olution or the whole domain, as sensors are not always fast enough and cannot be placed arbitrarily in systems that experience high temperatures and pressures, for instance [ 13 , 14 , 15 ]. Another challenge with multiphase flows is that some adopted sensors cannot be calibrated for both v apour and liquid, requiring ad hoc calibration and the ability to distinguish between the two [ 13 , 16 ]. Model predictive control (MPC) is an attractiv e approach for handling fast transients and managing constraints [ 17 , 18 , 19 ]. It is also worth mentioning the ability of such models to couple planning with constraints, as actuator limits, safety boundaries, and operational requirements can be effecti vely enforced while optimising performance over a finite horizon [ 20 ]. Ho wev er, multi- phase systems are strongly nonlinear, often exhibiting regime transitions and hysteresis, and generally inv olve tightly coupled interfacial physics across multiple scales. As a result, constructing reduced-order models or accurate surro- gates that capture all relev ant underlying features of the problem is an exceptionally difficult task [ 21 ]. At the same time, high-fidelity CFD is typically far too computationally expensi ve for real-time prediction and optimization [ 22 ]. Modelling such systems, in fact, often requires small internal time steps, repeated nonlinear/pressure-correction iter- ations, and interface reconstruction methods [ 23 ]. Embedding such solvers within the repeated rollouts required by MPC can therefore be challenging, particularly because updates must occur at the same time scale as sensing and ac- tuation. While the observ able of interest is often a low-dimensional v ariable, such as the lev el of liquid and/or the size distribution of an emulsion, that observable depends on the global ev olution of coupled PDEs and boundary forcing. As a result, computational ef ficiency is of paramount importance to make MPC a practically implementable approach. A recent work by [ 24 ] has introduced the concept of neural operators. Neural operators are a class of machine learning models that learn mappings between function spaces instead of finite-dimensional entities [ 25 ]. More in detail, while a general surrogate model like CNN or RNN learns a mapping f : R n → R m , a neural operator learns a mapping between infinite spaces G : A → U where A and U are spaces of functions. The latter have se veral advantages, most notably that they are theoretically discretisation-independent because input and output are functions [ 26 , 27 ], a feat that, while often mentioned and particularly appealing, presents some limitations in practice and is not what we e xploit in this work. The characteristic we are interested in is, in fact, the ability to infer the evolution of complex systems extremely quickly , making neural operators a v alid alternativ e to reduced-order models for real-time optimisation and control. Among NOs, Fourier Neural Operators (FNOs) are a class of neural operators that are particularly attractiv e for PDE systems, including the Navier -Stokes equations, since they represent global interaction well in the spectral domain (spectral methods are commonly used to solve the NS equations anyway) [ 28 , 26 ]. W e will therefore adopt FNOs despite the existence of other architectures that might be similarly effecti ve, such as DeepONet or U-Net [ 29 , 30 ]. As mentioned earlier , for control, the key advantage is that FNO surrogates provide fast, differentiable multi-step forecasts, enabling repeated horizon ev aluations at a cost orders of magnitude lower than CFD [ 31 , 32 ]. While we recognise differentiability as a particularly useful property of FNOs, this w ork does not exploit it, as explained later . The central idea of this effort is to couple an FNO surrogate with MPC to enable closed-loop regulation of multiphase dynamics. W e train an FNO to forecast the future ev olution of a phase-indicator field (here, the v olume fraction α ) over a finite horizon using a short history of recently sav ed states and a candidate control signal (here, inlet velocity). W e then map the predicted fields to the controlled scalar observ able (liquid lev el). The controller then selects an admissible actuation that tracks a time-varying setpoint schedule while enforcing hard bounds and penalising large control mo ves. The algorithm is tested on a bubble column reactor case in a problem that inv olves controlling the reactor lev el. It is worth noting that the current configuration has some limitations, such as the fact that, as mentioned by Bieker et al. [ 33 ], most industrial configurations do not provide access to a complete domain description, as measurements are usually localised and discontinuous. A possible ev olution of this approach might be the adoption of Physics-Informed Neural Operators, as initially introduced by [ 34 ], or other approaches that allow w orking with partial data. 2 Numerical Model Problem formulation The dynamics we intend to control are generally modelled with systems of partial differential equations. In the application described belo w , we solve a compressible multiphase flo w problem in which the two phases are solv ed independently , and the interface is tracked by transporting the volume fraction α i ( x , t ) . Rusche and other authors introduced the Eulerian-Eulerian multi-fluid formulation [ 35 , 36 ] later incorporated it into OpenFOAM 2 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 and named: twoPhaseEulerF oam . W e used the OpenFOAM solution algorithm to generate the training dataset and to simulate the controlled system. As anticipated, the phase volume fraction α i is transported with the following: ∂ α i ∂ t + ∇ ·  U α i  + ∇ · ( U r α i (1 − α i )) = 0 , i = 1 , 2 . (2.1) In which the time ev olution of α i depends on an advectiv e term and an interface-related transport contribution that is non-zero only around the interface while vanishing in the pure phase regions. The mixture velocity , weighted on the volume fraction, and the relati ve v elocity are defined as: U = α 1 U 1 + α 2 U 2 , U r = U 1 − U 2 (2.2) respectiv ely . The phase-wise Eulerian momentum equations are: ∂ ∂ t ( α i ρ i U i ) + ∇ · ( α i ρ i U i ⊗ U i ) + ∇ ·  α i ρ i R eff i  = − α i ∇ p + α i ρ i g + M i , i = 1 , 2 . (2.3) where ρ i and U i are the density and velocity of phase i , ∂ ( α i ρ i U i ) /∂ t is the accumulation of phase momentum, ∇· ( α i ρ i U i ⊗ U i ) is conv ectiv e transport of momentum with ⊗ denoting the dyadic (outer) product, ∇ · ( α i ρ i R eff i ) rep- resents the div ergence of the effecti ve stress (viscous plus turbulent/Reynolds stress), − α i ∇ p is the pressure-gradient force using a shared pressure field p , α i ρ i g is the gra vitational body force, and M i is the interphase momentum trans- fer (e.g., drag and other coupling forces) providing two-way coupling between phases. The effecti ve stress closure is R eff i = − ν eff i  ∇ U i + ( ∇ U i ) T − 2 3 I ( ∇ · U i )  + 2 3 I k i , ν eff i = ν i + ν t i . (2.4) with ν eff i the effecti ve kinematic viscosity (molecular ν i plus turbulent/eddy viscosity ν t i ), ∇ U i the velocity-gradient tensor , ( ∇ U i ) T its transpose, I the identity tensor , and ∇ · U i the v elocity di ver gence; the bracketed term corresponds to the traceless strain-rate form for a Ne wtonian stress model. In the remaining term k i is the turbulent kinetic energy of the i -th phase whose isotropic contribution is included in the normal stresses. Solution methodology As anticipated, the solver uses an Euler-Euler two-fluid formulation, treating both phases as continua. The volume fraction continuity equation is solved using the Multidimensional Universal Limiter with Explicit Solution (MULES) algorithm, which helps maintain the volume fraction within the bounds 0 and 1. The momentum equation is coupled with a Poisson-like pressure equation in the so-called ”PIMPLE” [ 37 ] algorithm, which is a combination of Pressure-Implicit Splitting of Operators (PISO) [ 38 ] and Semi-Implicit Method for Pressure- Linked Equations (SIMPLE) [ 39 ]. Neural Operators Neural Operators (NOs) are a class of operator-learning models that approximate mappings be- tween infinite-dimensional function spaces [ 28 ]: G : A (Ω; R m ) → W (Ω; R n ) , (2.5) where A (Ω; R m ) and W (Ω; R n ) are both Banach spaces, and G is a nonlinear operator mapping input functions a ∈ A to output functions w ∈ W . In the context of fluid dynamics, Ω ⊂ R d is the domain with d ∈ N ; the input and output vector spaces R m and R n may have different dimensions. In this setting, a ( x ) may represent an input field such as a coefficient, boundary condition, or source term, while w ( x ) is the corresponding output, typically the solution of a system of PDEs. Among NOs, the Fourier Neural Operator (FNO) is particularly suitable for modelling fluid-dynamics problems, since fluid fields usually exhibit spatial correlations that are well represented in the Fourier basis. FNOs learn an approximation b G θ ≈ G by composing layers that act globally across the entire domain Ω , making it efficient at capturing long-range spatial dependencies. In the definition of b G θ , θ denotes the trainable parameters. Fourier Neural Operators The methodology proposed by Li et al. [ 28 ] begins by lifting the input function a ( x ) into a higher-dimensional latent representation: v 0 ( x ) = P ( a ( x )) , (2.6) where a ( x ) : Ω → R m is the input function and P : R m → R d v is a learned local map that projects the input into the latent space R d v , enabling the network to capture more complex structures. The output v 0 ( x ) ∈ R d v serves as input to the subsequent operator layers, defined as: v l +1 ( x ) := σ ( W l v l ( x ) + ( K l v l ) ( x )) , ∀ x ∈ Ω , (2.7) 3 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 where v l ( x ) is the feature at location x in the l -th layer , l = 0 , . . . , L − 1 , W l : R d v → R d v is a learned local linear transformation, and σ is a nonlinear acti vation function applied componentwise to v ectors in R d v . The nonlocal operator K l is defined by: ( K l v l )( x ) = F − 1 ( R l ( k ) ˆ v l ( k )) ( x ) , ˆ v l = F ( v l ) , (2.8) where k ∈ Z d is a discrete multi-index of Fourier modes. In practice, F and F − 1 are implemented by the FFT on the computational grid, and only a finite set of low modes is retained: M l :=  k = ( k 1 , . . . , k d ) ∈ Z d : | k j | ≤ M l,j , j = 1 , . . . , d  , (2.9) with other modes k / ∈ M l set to zero. For a d v -channel latent field, the spectral multiplication is ( R l ˆ v l ) h ( k ) = d v X j =1 R l,hj ( k ) ˆ v l,j ( k ) , h = 1 , . . . , d v , k ∈ M l , (2.10) where R l ( k ) ∈ C d v × d v . The l -th FNO layer is therefore written as: v l +1 ( x ) = σ  W l v l ( x ) + F − 1 ( R l ( k ) ˆ v l ( k )) ( x ) + b l ( x )  , l = 0 , . . . , L − 1 . (2.11) where b l : Ω → R d v is an implementation-specific bias field. The final latent representation is projected to the output by w ( x ) = Q ( v L ( x )) , (2.12) where Q : R d v → R n is a learned local transformation; in our implementation, Q is linear . The model parameters θ are trained by minimising the empirical loss: L MSE ( θ ) = 1 N N X i =1 1 n N y N x    b Y i − Y i    2 F , (2.13) In the NMPC application described below , the generic input a is identified with an m -channel spatial field whose channels encode the state history , the coordinate channels, and the planned control sequence ¯ u k . On the computational grid, this field is represented as a channel-stacked tensor . The generic output w is identified with the predicted v olume- fraction trajectory ov er the horizon. In this setting, temporal context is incorporated through the channel dimension rather than by enlarging the spatial domain: Ω remains purely spatial, while both the K past snapshots and the H planned control moves are encoded through channels. This yields H d u control channels in total. The generic output w comprises the H predicted volume fraction snapshots { α k +1 , . . . , α k + H } , so that b G θ ( α k − K +1: k , ( x x , x y ) , ¯ u k ) ≈ α k +1: k + H , (2.14) with m = K + 2 + H d u input channels and n = H output channels. Neural Model Predictive Control In the following, we propose a neural model predictiv e control (NMPC) algo- rithm, in which the trained FNO surrogate b G θ is used to predict the future ev olution of the system and Bayesian optimization (BO) is used to select the control action minimizing a finite-horizon objectiv e. Let T s > 0 denote the sampling interval and H ∈ Z > 0 the prediction horizon. The plant is the PDE system introduced abov e, which gener- ates the volume fraction field α ( x, t ) on the spatial domain Ω ⊂ R d , with spatial coordinate x ∈ Ω , under the action of a scalar input u ( t ) ∈ U := [ u min , u max ] . The control is applied with zero-order hold, u ( t ) = u k , t ∈ [ k T s , ( k + 1) T s ) , (2.15) where u k ∈ U is the input applied at step k . Let the measured field at time step k be α k ( x ) := α ( x, k T s ) , (2.16) and define the surrogate conditioning v ariable as the stack of the K most recent snapshots, z k := col( α k − K +1 , . . . , α k ) ∈ [ L 2 (Ω)] K . (2.17) The controlled output, is extracted from the v olume fraction field through a functional: y k := ℓ ( α k ) ∈ R . (2.18) Let r k ∈ R denote the reference value at time step k . At each control step, the FNO surrogate receiv es as input the state history z k , the spatial coordinate channels ( x x , x y ) , and a planned input sequence ov er the prediction horizon, ¯ u k := col  u k | k , u k +1 | k , . . . , u k + H − 1 | k  ∈ U H , (2.19) 4 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 and returns the predicted future fields b G θ : ( z k , ( x x , x y ) , ¯ u k ) 7→ b α k +1: k + H := col  b α 1 | k , . . . , b α H | k  , b α i | k ≈ α k + i , (2.20) where the relativ e index i | k denotes the i -step-ahead prediction made at step k . In our implementation, the control sequence is embedded as boundary-supported input channels, so that the control acts only on the inlet boundary ∂ Ω in . The corresponding predicted output trajectory is b y i | k := ℓ  b α i | k  , b ¯ y k := col  b y 1 | k , . . . , b y H | k  . (2.21) Because only a single scalar inlet input is optimized, we parameterize the planned sequence as constant ov er the horizon, u k + i | k = u k , i = 0 , . . . , H − 1 , (2.22) so that ¯ u k = u k 1 H ∈ U H . (2.23) The optimal control problem therefore reduces to a one-dimensional search over u k ∈ U . At each time step k , we solve min u k ∈U J k ( u k ) := 1 H H X i =1  b y i | k − r k  2 + λ ( u k − u k − 1 ) 2 (2.24) s.t. b α k +1: k + H = b G θ ( z k , ( x x , x y ) , u k 1 H ) , (2.25) b y i | k = ℓ  b α i | k  , i = 1 , . . . , H, (2.26) u min ≤ u k ≤ u max , (2.27) where λ ≥ 0 penalizes aggressive input variations. W e solve ( 2.24 ) using BO because the level functional ℓ is discontinuous in the application considered below: specifically , the liquid level is obtained from α via a thresholding operation, which renders the resulting objectiv e not suitable for gradient-based optimisation. At each BO iteration, a Gaussian-process surrogate is fitted to the set of pre viously ev aluated pairs D k = n u ( i ) , J k ( u ( i ) ) o n i =1 , (2.28) yielding a posterior mean m k ( u ) and variance s 2 k ( u ) . An acquisition function, here chosen as expected improvement (EI), is then maximized to balance exploitation of lo w-cost regions and e xploration of uncertain ones: u ( n +1) = arg max u ∈U EI k ( u ) . (2.29) After a fixed number of BO iterations, we apply the minimizer u ⋆ k to the plant, measure the new field snapshot α k +1 and update the history z k +1 before repeating the procedure. Finally , when ℓ is smooth, gradient-based optimization is more appropriate. Indeed, b G θ is fully differentiable, as discussed in [ 40 ], and the gradient dJ k /du k could in principle be computed by automatic differentiation. In the present setting, howe ver , BO remains particularly attractive because the decision variable is scalar , thereby av oiding the curse of dimensionality that would otherwise limit its practicality in higher-dimensional MPC parameterisations. 3 A pplication: Bubble Column Reactor Control The framework detailed in the previous section is applied to the problem of controlling a bubble column reactor . The above consists of a system in which a column of liquid is put in motion by the injection of a gas introduced from the bottom ( 1 ). These reactors are designed to fav or heat and mass transfer between the gas and liquid phases by generating an extremely large interfacial area through the formation of small bubbles, usually through a sparger . All this complexity can only be accurately handled by sophisticated CFD models, since most diagnostics in real-life scenarios lack spatial or temporal resolution [ 41 ]. Designing a control system for such devices has therefore always been challenging, and the most common approach has been to either adopt passive control mechanisms or build control on sparse data or , e ventually , on reduced-order models. In the following paragraphs, we detail the computational domain and simulation parameters, the generated dataset information, and the training results. Finally , the control performance is highlighted. Computational domain and mesh The computational domain consists of a single structured hexahedral block. The domain spans L x = 0 . 15 m , L y = 1 . 0 m , and a finite out-of-plane thickness L z = 0 . 10 m (extruded mesh), with vertices defined at (0 , 0 , 0) , (0 . 15 , 0 , 0) , (0 . 15 , 1 , 0) , (0 , 1 , 0) and their counterparts at z = 0 . 10 m . The block is discretized into ( N x , N y , N z ) = (25 , 75 , 1) cells with uniform grading. The corresponding uniform cell spacings are ∆ x = L x / N x = 6 . 0 × 10 − 3 m , ∆ y = L y / N y ≈ 1 . 33 × 10 − 2 m , and ∆ z = L z / N z = 0 . 10 m . 5 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 Figure 1: Representativ e snapshots of the bubble-column phase field (volume fraction) under closed-loop operation, shown for different commanded liquid le vels. The interface location and the internal recirculation/b ubble-plume struc- tures evolv e nonlinearly with the inlet actuation, highlighting the strongly state-dependent h ydrodynamics captured by the CFD dataset. Initial and boundary conditions (two-phase fields) The simulations were initialized with a stratified phase distri- bution using the volume fraction of the gas phase, α air . A total of 750 cells were initialized with α air = 0 and 1125 cells with α air = 1 , corresponding to an interface located at y ≈ (30 / 75) L y ≈ 0 . 40 m (water below , air abov e). At the inlet we set α air to 1 (pure air), and a step function for the velocity . No-slip conditions ( U = 0 ) were enforced on the side walls for both phases. Dataset generation The dataset was created by randomly generating cases with different input velocities changing ov er discrete time intervals. The ev olution of the system was captured e very 0 . 1 s for a total of 100 steps over 10 seconds. For each run, we stored only the volume fraction time series alpha with shape ( T , N y , N x ) and the corre- sponding inlet-velocity signal with shape ( T , C ) (or ( T , ) for scalar control, reshaped to ( T , 1) ), where C is the control dimension. W e then se gmented these time series into fixed-length temporal sliding windo ws to form the learning pairs used for training. The input portion consists of the K most recent volume fraction fields, α t − K +1: t , together with the planned inlet sequence o ver the ne xt H control intervals, u t : t + H − 1 . The target is the future volume fraction trajectory α t +1: t + H . T o match the operator-learning architecture, we encode the network input as a channel-stacked tensor X t ∈ R C in × N y × N x with C in = K + 2 + H C , where the additional 2 accounts for the spatial coordinate channels ( x x , x y ) appended point-wise to the input. The first K channels are the previous fields. The control is encoded by reshaping the future sequence into H C channels and injecting it only at the inlet through a fixed spatial mask M ∈ { 0 , 1 } N y × N x (in our implementation, M = 1 on the bottom row corresponding to the inlet boundary ∂ Ω in at y = 0 , and 0 elsewhere). The abov e produced control channels of the form U ( h,i,j ) t = ˜ u ( h ) t · M i,j , h = 1 , . . . , H C , ( i, j ) ∈ { 1 , . . . , N y } × { 1 , . . . , N x } , (3.1) where ˜ u t ∈ R H C is the planned control sequence u t : t + H − 1 flattened into H C scalars, and M ∈ { 0 , 1 } N y × N x is the inlet mask defined abov e. These control channels are concatenated with the history and coordinate channels to form the full input tensor X t =  α t − K +1: t   ( x x , x y )   U t  ∈ R ( K +2+ H C ) × N y × N x . (3.2) The training target w as: Y t = α t +1: t + H ∈ R H × N y × N x . (3.3) T raining W e trained a Fourier Neural Operator (FNO) on the resulting windowed dataset to learn an operator that maps a short history of volume fraction fields and the planned inlet-velocity profile to future volume fraction snapshots. 6 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 0 25 50 75 100 125 150 175 200 Epoch 0.05 0.10 0.15 0.20 0.25 Val MSE lr 1e-3 lr 1e-4 lr 1e-5 lr 2e-3 lr 2e-4 lr 3e-3 lr 3e-4 lr 3e-5 lr 5e-4 lr 7e-4 Figure 2: V alidation mean squared error (MSE) as a function of training epoch for different learning rates. All con- figurations sho w a rapid initial drop in error followed by gradual con vergence tow ard a similar low-MSE plateau. Larger learning rates accelerate early con vergence, while smaller learning rates reduce the error more gradually; o ver - all, learning rates in the intermediate range provide the best balance between con ver gence speed and final validation performance. MSE MAE SSIM 0.0 0.2 0.4 0.6 Value Overall Metrics Val Test Figure 3: Comparison of overall prediction accuracy on the validation and test sets using mean squared error (MSE), mean absolute error (MAE), and structural similarity index (SSIM). The similar values across both splits suggest stable performance and good generalization of the forecasting model. T raining was executed on a single NVIDIA Quadro GV100 (32 GB) with PyT orch 2.10.0, CUDA 12.2, and driver 535.288.01. The final optimised configuration used a batch size of 512 with GPU-resident preloading and CuPy preprocessing [ 42 , 43 ]. W e performed tests at v arious learning rates achieving satisfactory results in a few epochs for a learning rate of 2e-4 2 , 3 . As illustrated in Fig. 4 , the learned operator reproduces the global interface dynamics reliably; ho wev er , it tends to under -resolve fine-scale structures within the liquid phase, leading to small discrepancies in the detailed volume fraction distrib ution. T raining acceleration benchmarks W e benchmarked training on the same V100 system across multiple batch sizes and optimization settings. The best configuration achieved a 3.09x speedup ov er the single-GPU baseline at a batch size of 512. Pre-loading the dataset to GPU memory incurred a one-time cost of 4.3 ms per sample (about 50 s for 11,900 samples) and reduced data-transfer overhead during training. Inference throughput peaks at batch size 256, with 0.041–0.051 ms per sample (19,711–24,508 samples/s). T racking quality and transient beha vior Figure 5 compares the measured liquid le vel, in blue, against the reference setpoint, in red, piecewise constant. Overall, the controller achiev es consistent tracking across repeated setpoint changes while exhibiting only short transients at switching times and is equilibrated by sudden changes in gas flowrate as shown in Fig. 6 . The response is well-damped: after each step change, the level con ver ges rapidly toward the new 7 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 Ground T ruth Prediction t+1 t+2 t+3 t+4 t+5 Figure 4: V olume fraction field forecasting over a 5-step horizon with the trained Fourier Neural Operator (FNO). T op row: ground-truth volume fraction fields at t + 1 to t + 5 from the CFD model. Bottom row: corresponding FNO predictions, demonstrating the surrogate’ s ability to propagate interface motion and large-scale structures. 0 25 50 75 100 125 150 175 200 Time 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Level Measurement Setpoint Figure 5: Closed-loop setpoint tracking performance for the bubble-column liquid lev el. The measured level (blue) is regulated to a piecewise-constant reference (red dashed) across repeated setpoint changes, with short transients at switching times and small steady-state offsets near the operating e xtremes. 8 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 Configuration Batch Time/Epoch (s) Speedup Baseline (single-GPU) 64 12.13 1.00x Optimized (preload+CuPy) 64 5.73 2.12x Baseline (single-GPU) 256 13.34 1.00x Optimized (preload+CuPy) 256 5.08 2.62x Baseline (single-GPU) 512 16.96 1.00x Optimized (preload+CuPy) 512 5.48 3.09x T able 1: Training time per epoch on a single NVIDIA Quadro GV100. The optimized configuration uses GPU-resident preloading and CuPy preprocessing. 0 25 50 75 100 125 150 175 200 Time 0.0 0.1 0.2 0.3 0.4 0.5 0.6 C o n t r o l u Figure 6: Control action applied by the surrogate-assisted MPC. The inlet velocity is implemented as a constant- hold input updated at each control interval; the optimizer selects bounded actuation levels that balance tracking error reduction with conservati ve input usage over the full e xperiment horizon. target without sustained oscillations, indicating that the closed-loop dynamics are stable over the entire test horizon. The most prominent errors occur immediately after a shift in the setpoint, when the system’ s characteristic response time prev ents faster adaptation of the le vel as better highlighted in 7a . Steady-state offsets and constraint effects Once the short transition phase is excluded, the remaining steady track- ing error is small for intermediate setpoints (e.g., 0 . 6 – 0 . 8 ), while a systematic bias becomes visible at the e xtremes. In particular , at the highest commanded le vel ( setpoint = 0 . 9 ) the response tends to settle slightly below the tar get (typ- ical mean lev el ≈ 0 . 88 ), whereas at the lowest commanded le vel ( setpoint = 0 . 5 ) the response can remain mar ginally abov e the reference in some plateaus. This behavior is consistent with the input being bounded. Importantly , the offset remains bounded and repeatable, suggesting the controller is robust but operates close to its feasible actuation en velope at the e xtremes. 4 Conclusions W e introduced a framew ork for controlling comple x multiphase phenomena that occur over short characteristic times. W e applied a model predictive control strategy with neural operators for fast forecasting of complex dynamics ex- ploiting Fourier Neural Operators’ ability to accurately and promptly predict the spatio-temporal evolution of systems described by PDEs. The method’ s low latency is essential in many engineering applications, including a bubble column reactor , which is the case we chose to validate the method. W e demonstrated the ability to adapt the gas inlet velocity to match a giv en lev el of gas-entraining liquid in the system, a feat that requires capturing the nonlinear relationship between the amount of gas that escapes the free surface and the amount introduced. Rather than relying on reduced- order models tied to a fixed discretisation or operating point, the controller ev aluates candidate actuations through field-lev el rollouts and maps the predicted states to a task-relev ant observable, such as the liquid level in this case, using a receding-horizon objectiv e with input bounds and regularisation. Across repeated piecewise-constant setpoint changes, we achiev ed consistent tracking with short transients at switching times and stable behaviour throughout the ev aluation horizon, indicating that the surrogate-assisted MPC can regulate nonlinear multiphase dynamics without embedding high-fidelity CFD in the optimisation loop. Notably , the characteristic time scale of the computation is extremely small, ranging from 1e-3 to 1e-2 at each time step throughout the optimization process. On the Quadro GV100 system, FNO inference throughput reaches 19,711–24,508 samples/s (0.041–0.051 ms per sample) depending 9 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 (a) (b) Figure 7: MPC performance for liquid-lev el control under a time-varying setpoint. (a) T op row: level tracking (left) and control input u (right). The measured lev el follows the setpoint closely , with only brief overshoot and undershoot during transitions, while the inlet-velocity control action remains bounded. (b) Bottom ro w: tracking error (left) and control increments ∆ u (right). The tracking error remains close to zero for most of the simulation, with larger deviations occurring mainly at setpoint changes, while the control increments indicate generally smooth actuator behavior with occasional sharp corrections. on batch size (T able 1 ) . Overall, the results support the use of learned operators as an enabler of real-time control while a direct wall-clock comparison against an embedded CFD-MPC implementation is left for future work. By forecasting the full field, the framew ork naturally generalises beyond scalar regulation. It provides a foundation for future extensions to partial observability [ 44 , 45 ], uncertainty-aware decision making [ 46 ], and physics-informed op- erator learning [ 47 , 48 ], which are expected to further improve robustness and applicability to industrial multiphase unit operations. W e also want to note that, beyond speed, forecasting a field rather than a scalar output might en- able more refined applications, such as controlling the distribution of the gas phase across particular structures (e.g., more or less activ e regions of a catalyst bed) in the system. While MPC ultimately ev aluates performance through a task-relev ant functional, an accurate prediction of that functional still requires capturing transport, accumulation, and interfacial motion over the whole domain. A field-le vel surrogate retains the spatial structure necessary to represent these mechanisms, while the output functional provides a consistent and robust bridge from predicted fields to control 10 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 decisions. Moreover , the differentiability of the surrogate enables (in principle) gradient-based optimization. Howe ver , this possibility has not been explored in this work. Acknowledgements The authors sincerely thank the open-source communities that made this work possible. In particular , we ackno wledge the developers and contributors of the Python ecosystem, and OpenFO AM community . This work used Python 3.11.14, PyT orch 2.10.0, CuPy 13.6.0, RAPIDS 25.12, CUDA 12.2, and NVIDIA driv er 535.288.01 on a single Quadro GV100 (32 GB). Special recognition is giv en to the developers of twoPhaseEulerFoam and related two- phase solvers. This research was supported by King Abdullah University of Science and T echnology (KA UST) . References [1] S Balachandar and John K Eaton. T urbulent dispersed multiphase flow . Annual r evie w of fluid mechanics , 42(1):111–133, 2010. [2] Martin W ¨ orner . Numerical modeling of multiphase flows in microfluidics and micro process engineering: a revie w of methods and applications. Micr ofluidics and nanofluidics , 12(6):841–886, 2012. [3] Jiafeng Y ao and Masahiro T akei. Application of process tomography to multiphase flow measurement in indus- trial and biomedical fields: A revie w . IEEE Sensors J ournal , 17(24):8196–8205, 2017. [4] Detlef Lohse. Fundamental fluid dynamics challenges in inkjet printing. Annual r eview of fluid mechanics , 54(1):349–382, 2022. [5] Sriram Manoharan, Steven Summerville, Lucas Freiberg, Matthe w Coblyn, Jad G T ouma, Goran Jov anovic, and Brian K Paul. Manufacturing process design of a micro-scale liquid-liquid extractor and multi-phase separator . Journal of Manufacturing Pr ocesses , 56:1381–1391, 2020. [6] F atemeh Goodarzi and Sohrab Zendehboudi. A comprehensive revie w on emulsions and emulsion stability in chemical and energy industries. The Canadian Journal of Chemical Engineering , 97(1):281–309, 2019. [7] Dale D. McClure, Nora Aboudha, John M. Kavanagh, David F . Fletcher, and Geof frey W . Barton. Mixing in bubble column reactors: Experimental study and CFD modeling. Chemical Engineering Journal , 264:291–301, 2015. [8] Xu Zhao, Hang Ren, and Long Luo. Gas bubbles in electrochemical gas ev olution reactions. Langmuir , 35(16):5392–5408, 2019. [9] Y ong Hu, Hernan Olguin, and Ev a Gutheil. A spray flamelet/progress variable approach combined with a trans- ported joint pdf model for turbulent spray flames. Combustion Theory and Modelling , 21(3):575–602, 2017. [10] Jihun Oh and Carl V Thompson. The role of electric field in pore formation during aluminum anodization. Electr ochimica Acta , 56(11):4044–4051, 2011. [11] Y ang Liu and N Dinh. T reatment of nucleation and bubble dynamics in high heat flux boiling. Experimental Thermal and Fluid Science , 26(6-7):793–810, 2002. [12] Simon Alberti and Dorothee Dormann. Liquid–liquid phase separation in disease. Annual re view of genetics , 53(1):171–194, 2019. [13] Jamal Chaouki, Fa ¨ ıcal Larachi, and MP Dudukovic. Non-in vasive monitoring of multiphase flows . Elsevier , 1997. [14] W uqiang Y ang. Sensors and instrumentation for monitoring and control of multi-phase separation. Measur ement and contr ol , 39(6):178–184, 2006. [15] Aluisio do Nascimento Wrasse, Eduardo Nunes dos Santos, Marco Jose da Silva, Hao W u, and Chao T an. Capacitiv e sensors for multiphase flo w measurement: A revie w . IEEE Sensors Journal , 22(22):21391–21409, 2022. [16] Xue wei Shi, Chao T an, Feng Dong, Eduardo N dos Santos, and Marco Jose da Silva. Conductance sensors for multiphase flow measurement: A revie w . IEEE Sensors journal , 21(11):12913–12925, 2020. [17] Eleni Aggelogiannaki and Haralambos Sarimveis. Nonlinear model predictive control for distributed parameter systems using data driv en artificial neural network models. Computers & Chemical Engineering , 32(6):1225– 1237, 2008. 11 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 [18] Lars Gr ¨ une. Analysis and design of unconstrained nonlinear mpc schemes for finite and infinite dimensional systems. SIAM J ournal on Contr ol and Optimization , 48(2):1206–1228, 2009. [19] Hassan Arbabi, Milan K orda, and Igor Mezi ´ c. A data-driven koopman model predicti ve control framew ork for nonlinear partial differential equations. In 2018 IEEE Conference on Decision and Contr ol (CDC) , pages 6409–6414. IEEE, 2018. [20] Eduardo F Camacho and Carlos Bordons. Constrained model predictiv e control. In Model predictive contr ol , pages 177–216. Springer , 2007. [21] Himakar Ganti and Prashant Khare. Data-driven surrogate modeling of multiphase flows using machine learning techniques. Computers & Fluids , 211:104626, 2020. [22] Zhe W u, Andres Aguirre, Anh T ran, Helen Durand, Dong Ni, and Panagiotis D Christofides. Model predictiv e control of a steam methane reforming reactor described by a computational fluid dynamics model. Industrial & Engineering Chemistry Resear ch , 56(20):6002–6011, 2017. [23] Johan Roenby , Bjarke Eltard Larsen, Henrik Bredmose, and Hrvoje Jasak. A new volume-of-fluid method in openfoam. In VII International Conference on Computational Methods in Marine Engineering. Nantes: Inter- national Center for Numerical Methods in Engineering , 2017. [24] Zongyi Li, Nikola Ko vachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andre w Stuart, and Anima Anandkumar . Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03485 , 2020. [25] Nikola B K ov achki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew M Stu- art, and Anima Anandkumar . Neural operator: Graph kernel network for pdes. arXiv pr eprint arXiv:2103.12806 , 2021. [26] Nikola Ko vachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar . Neural operator: Learning maps between function spaces with applications to pdes. JMLR , 24(1), 2023. [27] Zongyi Li, Nikola K ovachki, Chris Choy , Boyi Li, Jean Kossaifi, Shourya Otta, Mohammad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Azizzadenesheli, et al. Geometry-informed neural operator for large-scale 3d pdes. Advances in Neural Information Pr ocessing Systems , 36:35836–35854, 2023. [28] Zongyi Li, Nikola K ov achki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew M Stuart, and Anima Anandkumar . Fourier neural operator for parametric partial differential equations. arXiv pr eprint arXiv:2010.08895 , 2020. [29] Lu Lu, Pengzhan Jin, and Geor ge E Karniadakis. Learning nonlinear operators via deeponet based on the univ ersal approximation theorem. Nat. Mach. Intell. , 3:218–229, 2021. [30] Olaf Ronneberger , Philipp Fischer , and Thomas Brox. U-Net: Con volutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted Intervention (MICCAI) , pages 234–241. Springer , 2015. [31] P aulo Alexandre Costa Rocha, Samuel Joseph Johnston, V ictor Oliv eira Santos, Amir A Aliabadi, Jesse V an Griensven Th ´ e, and Bahram Gharabaghi. Deep neural network modeling for cfd simulations: Benchmarking the fourier neural operator on the lid-driv en cavity case. Applied Sciences , 13(5):3165, 2023. [32] Ali Sayghe, Mohammed Mousa, Salem Batiyah, and Abdulrahman Husawi. Fourier neural operators for fast multi-physics sensor response prediction: Applications in thermal, acoustic, and flow measurement systems. Sensors , 26(4):1165, 2026. [33] Katharina Bieker , Sebastian Peitz, Stev en L Brunton, J Nathan Kutz, and Michael Dellnitz. Deep model predic- tiv e flow control with limited sensor data and online learning. Theoretical and computational fluid dynamics , 34(4):577–591, 2020. [34] Zongyi Li, Hongkai Zheng, Nikola Ko vachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar . Physics-informed neural operator for learning partial differential equations. ACM/IMS Journal of Data Science , 1(3):1–27, 2024. [35] Henrik Rusche. Computational fluid dynamics of dispersed two-phase flow at high phase fractions. Ph. D. thesis, University of London , 2002. [36] P aulo J Oliveira and Raad I Issa. Numerical aspects of an algorithm for the eulerian simulation of two-phase flows. International Journal for Numerical Methods in Fluids , 43(10-11):1177–1198, 2003. [37] Henry G W eller , Gavin T abor , Hrv oje Jasak, and Christer Fureby . A tensorial approach to computational contin- uum mechanics using object-oriented techniques. Computers in physics , 12(6):620–631, 1998. 12 A P R E P R I N T - M A R C H 2 7 , 2 0 2 6 [38] Raad I Issa. Solution of the implicitly discretised fluid flow equations by operator -splitting. Journal of computa- tional physics , 62(1):40–65, 1986. [39] Suhas V Patankar . A calculation procedure for two-dimensional elliptic situations. Numerical heat transfer , 4(4):409–425, 1981. [40] Ze Cheng, Zhuoyu Li, Xiaoqiang W ang, Jianing Huang, Zhizhou Zhang, Zhongkai Hao, and Hang Su. Accel- erating pde-constrained optimization by the derivati ve of neural operators. arXiv preprint , 2025. [41] Y ihuan Y an, Krishna Mohanarangam, William Y ang, and Jiyuan T u. Experimental measuring techniques for industrial-scale multiphase flow problems. Experimental and Computational Multiphase Flow , 6(1):1–13, 2024. [42] Ryosuke Okuta, Y uya Unno, Daisuke Nishino, Shohei Hido, and Crissman Loomis. Cupy: A numpy-compatible library for n vidia gpu calculations. In Pr oceedings of W orkshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Confer ence on Neural Information Pr ocessing Systems (NIPS) , 2017. [43] Adam Paszk e and et al. Pytorch: An imperative style, high-performance deep learning library . In Advances in Neural Information Pr ocessing Systems 32 , pages 8024–8035. Curran Associates, Inc., 2019. [44] Frank L Lewis and Kyriakos G V amvoudakis. Reinforcement learning for partially observable dynamic pro- cesses: Adaptiv e dynamic programming using measured output data. IEEE T ransactions on Systems, Man, and Cybernetics, P art B (Cybernetics) , 41(1):14–25, 2010. [45] Ste ven L Brunton, Joshua L Proctor , and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Pr oceedings of the National Academy of Sciences , 113(15):3932– 3937, 2016. [46] Christopher B ¨ ulte, Philipp Scholl, and Gitta Kutyniok. Probabilistic neural operators for functional uncertainty quantification. arXiv pr eprint arXiv:2502.12902 , 2025. [47] Sumanth Kumar Boya and Deepak Subramani. A physics-informed transformer neural operator for learning generalized solutions of initial boundary value problems. arXiv pr eprint arXiv:2412.09009 , 2024. [48] Y usuke Y amazaki, Ali Harandi, Mayu Muramatsu, Ale xandre V iardin, Markus Apel, Tim Brepols, Stefanie Reese, and Shahed Rezaei. A finite element-based physics-informed operator learning frame work for spatiotem- poral partial differential equations on arbitrary domains. Engineering with Computers , 41(1):1–29, 2025. 13