PRES22_0231.docx DOI: 10.3303/CET2294176 Paper Received: 13 May 2022; Revised: 02 July 2022; Accepted: 02 July 2022 Please cite this article as: Yegenova A., Izmbergenov N., Sultanov M., Brener A., 2022, Computer Simulation of Nonlinear Waves in Liquid Films with Surface Activity and Mass Sources at the Bottom, Chemical Engineering Transactions, 94, 1057-1062 DOI:10.3303/CET2294176 CHEMICAL ENGINEERING TRANSACTIONS VOL. 94, 2022 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Petar S. Varbanov, Yee Van Fan, Jiří J. Klemeš, Sandro Nižetić Copyright © 2022, AIDIC Servizi S.r.l. ISBN 978-88-95608-93-8; ISSN 2283-9216 Computer Simulation of Nonlinear Waves in Liquid Films with Surface Activity and Mass Sources at the Bottom Aliya Yegenovaa, Nurjan Izmbergenova, Murat Sultanovb, Arnold Brener a,* aAuezov University of South Kazakhstan, Tauke Khan 5, Shymkent, Kazakhstan bAkhmet Yassawi University, Turkistan arnold.brener@auezov.edu.kz The paper's primary goal and scientific contribution is the mathematical simulation of nonlinear surface waves in moving thin layers of perfect liquids with accounting for the combined influence of mass sources and surface activity at the free-moving boundaries. The original inverse numerical method has been developed and implemented in a computer code. The main results of the work lie in establishing certain conditions when, under the weak intensity of the bottom sources and free surface activity, the propagation of solitary waves with slowly changing characteristics and surface ripple is possible. The estimate of orders of the main control parameters in the mathematical model was also given. For selected sets of control parameters, computer simulation has been carried out, and the shapes of the surface of propagating nonlinear waves have been calculated. It was established the amplitude of the carrier wave increased 1.5 times by the joint influence of the weak mass source and surface activity at the moving boundary. In addition, the ripple intensity increased significantly, accounting for surface activity. 1. Introduction The nature of the propagation and evolution of nonlinear waves characteristics in systems with mass sources substantially depends on both the intensity of the sources and the form of the boundary conditions. The previous work showed that the structure of the waves propagation equation for thin liquid layers in the presence of mass sources could be destroyed by the perturbation on the right-hand side of the KdV-type equation (Brener et al., 2020). The fundamental structure of the equations was established, and initial asymptotic analysis was carried out to estimate the orders of control parameters. At the same time, computer simulation in the mentioned work was not carried out, and the effect of control parameters on the characteristics of propagating nonlinear waves was not studied. The literature review and the known results analysis show that the mass gain in the condensate film, for example, and the no isothermal nature of the process, can have a significant influence on the stability of the waveless flow of a thin layer of liquid (Tiwari and Davis, 2010). Despite many works devoted to the formation and propagation of nonlinear waves in liquid layers, almost all deal with constant flow rates. The experimental work of Kapitza (Dietze, 2016) made an outstanding fundamental contribution in understanding the capillary waves. This work’s unique contribution lies in the study of the role of capillary waves on the free surface. However, the flow problem with variable consumption remains largely open (Abcha et al., 2018). This is explained not by the low significance of studies of flows with variable fluid flow and variable physical properties but by a wide variety of observed effects of nonlinearity and dispersion (Karpman, 2011). At the same time, some factors make it difficult to use the well-known technique correctly applied to the flows with mass sources (Ivanova and Ivanov, 2016). These difficulties are related to the fact that the known methods are based on the existence of constant solutions of the basic unperturbed transport equations in the stationary case. In the presence of mass sources, the flow rate changes (increases or decreases), and there are no essential solutions to the constant (Onorato et al., 2010). As shown in some old works, (Newell,1985) and recently (Brener et al., 2020), even in the case of a constant flow rate but a variable bottom profile (Newell, 1985), the correct analysis of the situation turns out to be very complicated and time-consuming and leads to evolutionary equations that do not have a soliton-type solution. Recently, exact wave solutions of high-order nonlinear equations have been 1057 in active search. Such equations appear naturally in models of fluid motion in shallow water and thin layers with mass sources. Here it is necessary to mention a very recent article (Arnous et al., 2022), which considers a fifth- order equation with non-linearity of the convective type and with possible source terms. However, despite high- order derivatives in the considered equations, an analysis shows that the dispersion relations for the considered equations do not contain imaginary parts. In this respect, it can be assumed that the nature of the propagation of nonlinear waves will be similar to that of the KDV equation solutions (Yazhou Shi, 2018). However, in some accurate models that consider surface activity’s influence on wave formation (Brener et al., 2020), features make it difficult to apply the methods used in the mentioned work (Arnous et al., 2022). The next section of this paper discusses just such a situation. In this work, a more detailed analysis is carried out in order to correctly establish the type of equations describing thin-layer flows with mass sources and weak surface activity at moving boundaries. The scientific contribution and novelty of the results lie in the asymptotic analysis and computer simulation of flows in a thin layer of liquid in the presence of sources of mass and surface activity at the free boundary. The joint influence of the source of mass and surface activity on the characteristics of the wave flow is studied. 2. Theoretical details For the purpose of writing formulas more compactly, partial derivatives are denoted as subscripts for functions wherever this should not lead to misunderstandings. The equation of continuity for potential flow ϕ of perfect liquid with free surface over a plate reads (Figure 1) (Newell, 1985). 0=+ yyxx ϕϕ . (1) The boundary conditions at the bottom with allowing for the bottom mass source read ( )xhyqh yxx −==+⋅ ;ϕϕ (2) The density of mass stream q depends on the actual mechanism of physical processes in the neighbourhood of the bottom (Figure 1). The simplest form of the appropriate dependence that can be obtained from the boundary conditions reads τVkVn ⋅= . (3) Here nV is the normal component of liquid velocity nearby the bottom and τV is the tangent component. Figure 1: Scheme of liquid flow with mass source at the bottom Condition (3) can be interpreted as a velocity of washout, and it has physical meaning for a pseudo-perfect fluid when the adhesion condition is not set. It can be established (Brener et al., 2021) that ( )xyxyxx hkh ⋅+=+⋅ ϕϕϕϕ . (4) In the long-wave approach for surface perturbations, the wave length l is supposed to be much bigger than average thickness of the liquid layer 0h . Such an assumption is quite correct for thin layers, when the characteristic longitudinal scale prevails (Iele et al., 2020). In addition, suppose that the perturbation amplitude is also small that is necessary in order to stay in framework of weak nonlinearity. The set of small parameters reads as follows (Brener et al., 2020): 1220 <<= µlh , 10 <<= εha , 2 11 εµε ⋅≈⋅⋅= kkk . In framework of weak non-linearity, the bottom shape is a function of the slow variable (Newell, 1985) xX ε= . The following dimensionless variables are suitable: x y nV τV U 1058 lxx ⋅→ , 0 0 hg h a ⋅⋅⋅→ϕϕ , 0hg l tt ⋅ ⋅→ , 0hyy ⋅→ , a⋅→ηη , 0hhh ⋅→ . (5) The Taylor expansion of the velocity potential in the vicinity of the bottom gives ( ) ( ) ( ) ( ) ( ) ... 24 1 6 1 2 1 , 432 ++⋅++⋅++⋅+++= hyhyhyhytxF yyyyyyyyyy ϕϕϕϕϕ . (6) Using (6) and neglecting terms of higher than ε orders, (13), (14) take the forms 2 xxxtt FFHF ⋅−⋅⋅+−= εεη , (7) ( )     ⋅−⋅−⋅−⋅+⋅−=⋅− 321 6 1 2 HFFFFFFHFkHFHF xxxxxxtxtxxxttxxxxtt ε , (8) where hH +=1 . For more comfortable use of the methods of secular perturbations theory, it is reasonable to look for the function ( )txF , in the form (Brener et al., 2020) ( ) ( )2O,,, ε ε θ ε ε θ +      +      = XfXFtxF , (9) where θ is a self-similar variable depending on so-called slow coordinates (Newell, 1985) xX ε= , tT ε= . In order to satisfy this equation in a zero order the following dispersion relation should be fulfilled 022 =⋅− XT H θθ . (10) For eliminating secular terms in the next order, it is necessary to suppose ( ) U H HkH UH H UU H U X XXXXT XT XTX X ⋅      ⋅ −⋅+⋅− =⋅      ⋅− ⋅ +⋅⋅ ⋅ ⋅− θ θθθ θθ θθθ θθθθ 23 1 42 3 122 , (11) where θFU = . It can be concluded that for the accepted order of mass source intensity at the bottom its influence on non-linear waves propagation in thin liquid layer may be described by Eq(11). Under choosing 0Xθ this equation has a structure of the left side of Eq(11) which bears a resemblance to the Korteweg-de-Vries equation (KdV) (Tsvelodub et al, 2021). However, in any case, the complete structure of KdV equation is destroyed by the perturbation in the right-hand side, and that can be interpreted as damping influence of mass source (Tseluiko and Kalliadasis, 2011). Structure of the propagation equation for non-linear waves may be also essentially changed under the influence of phenomena taking place at the free surface (Dietze, 2016). Namely, Eq(11) after the rearrangements and transforms which pursue an aim to eliminate secular growth of perturbations can be expanded by the term describing a surface activity (Brener et al., 2020) in the following form ( ) . 2 23 1 42 3 3 122 θθ θθθθ θ θσ θ θθθ θθ θθθ U H U H HkH UH H UU H U X T X XXXXT XT XTX X ⋅ ⋅⋅ ⋅ +⋅      ⋅ −⋅+⋅− =⋅      −− ⋅ +⋅⋅ ⋅ ⋅− (12) A preliminary analysis of the references allows suggesting in this case possibility of the development of solitary waves, as well as the wave breaking at various ratios of control parameters (Trifonov, 2017). The presence of the perturbation described by second-order derivative leads to existence in the equation of derivatives of 1059 different total parity, which in turn leads to the fact that the dispersion relation contains a nonzero imaginary part. An undamped wave solution can exist only on the neutral line that is accompanied by the amplitude growth (Brener et al., 2020). By this, the possibility of obtaining exact solutions using the methods of work (Arnous et al, 2022) also becomes doubtful. 3. Computer simulation The well-conditioned numerical methods do not work well with third-order equations, as the error correct evaluation becomes difficult (Hossenny, 2012). In the course of the analysis of main equation, the following features can be mentioned. Second and third order derivatives exist only for the modified time variable; and only the first derivative for a modified space variable is used. Taking into account these features, an inverse method for the solutions behavior analyzing has been used. It is based on the Finite Differences method and the uncertainty approach, which made it possible to submit the specific algorithm, to write a code for numerical analysis and built the characteristic plots of the sought-for function. In this case, Eq(12) can be written in the form θθθθθθ UDUCUBUUAU X ⋅+⋅=⋅+⋅⋅+ . (13) Here DCBA , , , are slowly changing functions, which, in the first approximation, are taken as constants on the studied intervals of the spatial and temporal accepted arguments. Control parameters in Eq(13) can be easily identified by comparing Eq(12) and Eq(13). At the same time, when setting up a numerical experiment in this work, the control parameters were not determined based on the values of their constituent physical quantities. This decision was made because at the first stage of the study, the task was to investigate the qualitative features and the general picture of the flow of liquid layer. Therefore, the control parameters were chosen based on considerations of values that were acceptable in terms of orders. 3.1 Description of the algorithm of numerical experiment The first step of the algorithm is to generate a pseudo-random function )(θU , at the fixed point 0X . For calculating the derivatives of the first, second and third order the well-known implicit Euler scheme has been used. In order to end up with some rectangular matrix of values of the function ),( XU θ that satisfy the solution of the equation, the first set of dependencies with a number of values exceeding the size of the final solution matrix should be generated. Indeed, for calculating the derivative using the Finite Difference method requires two values of the original function. At the last point of the generated sequence, it is not able to calculate the derivative. The extra point should be discarded in subsequent calculations, and since the derivatives up to the third order should be calculated, that 3 points will be discarded at step 2. The initial sequence must be of such length that the desired size of the rectangular decision matrix remains acceptable at each step in X . Then, to get a rectangular solution matrix, it is enough to discard redundant points to bring the matrix to the required form. The second step is the calculation of the first, second and third derivative with respect to time. As it is described above, the extra points are discarded, and the lengths of the sequences decrease each time. The optimal step for the modified time variable is selected by the way of iterative procedure and by comparing the calculated values with those obtained using the approach of Doubly Recursive Multivariate Automatic Differentiation (Kalman, 2002). The third step is the choice of values for parameters A, B, C, D. By using ready- made values of the time-derivatives as terms of the equation, the selection of values of parameters A, B, C, D does not univocal, and it is carried out arbitrarily from considerations of an acceptable order of parameters magnitude. The fourth step is the calculation of the spatial derivative from the original equation. The sought derivative with respect to X can be calculated from the original equation having the set parameters A, B, C, D with the help of calculated at the previous step all three time derivatives at the point 0X . A series of values XU ∂∂ for all points of the time variable θ can be defined. The fifth step is the calculation of the U value at point XXX ∆+= 01 . Now, having a series of XU ∂∂ values for all points of the time variable θ and taking the acceptable step for the spatial variable dx, all the differences between the values of U at point X0 and point X1 can be calculated. That is, all values of sought function ),( 1XU θ for the entire time row can be calculated also. As a result, the series of parallel values ),( 1XU θ for the next point X1 will be calculated. In this way, the entire surface that describes the function can be obtained. 1060 Figure 2: Flowchart After constructing the entire surface that describes the function in this way, it is necessary to calculate the residuals (the differences between the left and right sides of Eq(13)) over the entire surface based on a comparison with the control equation. If the maximum residual over the entire surface is higher than the given value, you must return to step 3 in order to correct the parameter values. Through such iterations, an acceptable maximum of residuals is achieved. 3.2 Results of simulation Figures 3, 4 depict some characteristic results of the numerical experiment. The numerical study showed a clear increase both in the waves amplitude and in the contribution of small waves of the "ripple" type on the surface of carrier waves with an increase in the coefficient D, which is responsible for the effect of surface activity, i.e. influence of capillary forces, on the propagation of waves. I II Figure 3: Characteristic wave plots (I- U ; II- XU ∂∂ / ) for numerical solution of Eq(13). Control parameters set M1: A = 2.9E-4, B = 4.5E-8, C = 1.8E-5, D = 2.3E-6. Step Time = 0.01, Step X = 0.1, Points:100. I II Figure 4: Characteristic wave plots (I- U ; II- XU ∂∂ / ) for numerical solution of Eq(13). Control parameters set M2: A = 2.9E-4, B = 4.5E-8, C = 1.8E-5, D = 0. Step Time = 0.01, Step X = 0.1, Points:100 As shown by the numerical experiment (Figures 3 and 4), when the value of the control parameter D, which is responsible for surface activity (capillary forces), was different from zero, the amplitude of the carrier wave increased 1.5 times. The ripple intensity increased significantly too. The contribution of the surface pressure gradient, which is due to the variable curvature of the film, manifests itself as a main dispersive factor. Waves of higher orders evolve, but are not, strictly speaking, dispersive. Upon reaching a certain D, destruction of a single wave become probably possible. 1061 4. Conclusions As a result of considering the problem of the flow of a thin liquid layer over an inclined surface in the long- wavelength approximation, small parameters have been singled out by asymptotic analysis and their hierarchy has been established. An original algorithm for a numerical experiment has been developed and computer simulation of nonlinear waves for various sets of control parameters has been carried out. In the calculations, a sharp increase in the amplitude of the carrier wave by 1.5 times was observed when changing the parameter D from zero to 6103.2 −×=D . The intensity of small surface waves of high frequency, the so-called "ripple", increases essentially too. This phenomenon can be interpreted as the formation of a single wave with a simultaneous increase in the amplitudes of the surface ripples, which can lead to the destruction of the structure of the carrier wave. Nomenclature a – wave amplitude g – gravitational acceleration k – a proportionality factor l – wave amplitude x – longitudinal coordinate y – normal coordinate References ),( txη – the perturbed free surface of the liquid ϕ – velocity potential σ – the coefficient of surface activity q – the density of the mass stream across the solid bottom t – time Abcha N., Pelonovsky E., Mutabazi I., 2018, Nonlinear Waves and Pattern Dynamics, Cham, Switzerland, Springer. Arnous A.H., Mirzazadeh M., Akinyemi L, Akbulut A., 2022, New solitary waves and exact solutions for the fifth- order nonlinear wave equation using two integration techniques, Journal of Ocean Engineering and Science, doi: 10.1016/j.joes.2022.02.012. Brener A.M., Yegenova A., Botayeva S., 2020, Equations of Nonlinear Waves in Thin Film Flows with Mass Sources and Surface Activity at the Moving Boundary, WSEAS Transactions on Fluid Mechanics, 15, 149- 162. Dietze G.F., 2016, On the Kapitza Instability and the Generation of Capillary Waves, Journal of Fluid Mechanics, 789, 368-401. Hossenny C., 2012, Numerical Methods for the KdV Equation and Related Water Waves, PhD Thesis, University of Technology, Mauritius. Iele B., Palleschi B., Gallerano F., 2020, Numerical Modelling of Wave Fields and Currents in Coastal Area, WSEAS Transactions on Fluid Mechanics, 15, 41-53. Ivanova O. and Ivanov M., 2016, 1-D Mathematical Modelling of Flood Wave Propagation, Chemical Engineering Transactions (CET), 53, 139-144. Kalman D., 2002, Doubly Recursive Multivariate Automatic Differentiation, Mathematic Magazine, 75(3), 187- 202. Karpman V.I., 2011, Nonlinear Waves in Dispersive Media, Pergamon Press, Oxford, New York, Toronto, Sydney. Newell A. C., 1985, Solitons in mathematics and physics, Society for Industrial and Applied Mathematics, University of Arizona, USA. Onorato M., Proment D., Toffoli O., 2010, Freak waves in crossing seas, The European Physical Journal Special Topics, 185(1), 45-55. Tiwari N. and Davis J.M., 2010, Stabilization of thin liquid films flowing over locally heated surfaces via substrate topography, Physics of Fluids, 22, 1-12, 042106. Trifonov Yu., 2017, Nonlinear waves on a liquid film falling down an inclined corrugated surface, Physics of Fluid, 29, American Institute of Physics, 054104-(1-7). Tseluiko D. and Kalliadasis S., 2011, Nonlinear waves in counter-current gas–liquid film flow, Journal of Fluid Mechanics, 673, Cambridge University Press, 19 – 59. Tsvelodub O.Yu., Arkhipov D.G., Vozhakov I.S., 2021, Investigating waves on the surface of a thin liquid film entrained by a turbulent gas flow: modeling beyond the “quasi-laminar” approximation, Thermophysics and Aeromechanics, 28, 223–236. Yazhou Shi, Xiangpeng Li, Ben-gong Zhang, 2018, Traveling Wave Solutions of Two Nonlinear Wave Equations by -Expansion Method, Advances in Mathematical Physics, 2018, 1-8. 1062 https://link.springer.com/article/10.1134/S0869864321020050#auth-O__Yu_-Tsvelodub https://link.springer.com/article/10.1134/S0869864321020050#auth-D__G_-Arkhipov https://link.springer.com/article/10.1134/S0869864321020050#auth-I__S_-Vozhakov https://link.springer.com/journal/11510 https://link.springer.com/journal/11510 PRES22_0362.pdf Computer Simulation of Nonlinear Waves in Liquid Films with Surface Activity and Mass Sources at the Bottom