Next Article in Journal
Characterizations of the Frame Bundle Admitting Metallic Structures on Almost Quadratic ϕ-Manifolds
Previous Article in Journal
A Long-Term Prediction Method of Computer Parameter Degradation Based on Curriculum Learning and Transfer Learning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Role of Elasticity on Chaotic Dynamics: Insights from Mechanics, Immunology, Ecology, and Rheology

by
Ángela Jiménez-Casas
1,
Mario Castro
1,2,3 and
Manuel Villanueva-Pesqueira
1,*
1
Grupo de Dinámica No Lineal, Universidad Pontificia Comillas, c. Alberto Aguilera 25, E28015 Madrid, Spain
2
Instituto de Investigación Tecnológica (IIT), Universidad Pontificia Comillas, 28015 Madrid, Spain
3
Grupo Interdisciplinar de Sistemas Complejos (GISC), 28015 Madrid, Spain
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(14), 3099; https://doi.org/10.3390/math11143099
Submission received: 15 June 2023 / Revised: 10 July 2023 / Accepted: 11 July 2023 / Published: 13 July 2023
(This article belongs to the Section Difference and Differential Equations)

Abstract

:
Elasticity is commonly associated with regular oscillations, which are prevalent in various systems at different scales. However, chaotic oscillations are rarely connected to elasticity. While overdamped chaotic systems have received significant attention, there has been limited exploration of elasticity-driven systems. In this study, we investigate the influence of elasticity on the dynamics of chaotic systems by examining diverse models derived from mechanics, immunology, ecology, and rheology. Through numerical MATLAB simulations obtained by using an ode15s solver, we observe that elasticity profoundly alters the chaotic dynamics of these systems. As a result, we term the underlying equations as the elastic-Lorenz equations. Specifically, we extensively analyze a viscoelastic fluid confined within a closed-loop thermosyphon, considering general heat flux, to demonstrate the impact of the viscoelastic parameter on the model’s chaotic behavior. Our findings build upon prior research on the asymptotic behavior of this model by incorporating the presence of a viscoelastic fluid. The results highlight the non-trivial and non-monotonic role of elasticity in understanding the control, or lack thereof, of chaotic behavior across different scales.
MSC:
35B40; 35Kxx; 35Qxx; 65L05

1. Introduction

The study of chaos theory has been a significant field in both physics and mathematics for over half a century. It has provided valuable insights into a range of natural phenomena, including weather patterns and celestial motion. One of the most notable examples of chaotic systems is the groundbreaking Lorenz model, developed by Edward Lorenz in the early 1960s [1]. This model, which describes the behavior of a simple system composed of three interconnected nonlinear ordinary differential equations, has been used to explain a wide variety of phenomena, such as fluid convection and airplane wing oscillation [2].
Despite its success, the Lorenz model is just one example of a chaotic system. Many other systems may also exhibit chaotic behavior [3]. Chaotic systems have complex dynamics and are sensitive to initial conditions, meaning that even slight changes in these conditions can lead to vastly different outcomes over time. This sensitivity presents challenges for predicting the long-term behavior of chaotic systems, which has implications for fields like meteorology and climate modeling [4] and characterizing problem complexity [5].
Recently, researchers have begun exploring the role of elasticity in chaotic systems [6]. Elasticity refers to the ability of materials to deform and return to their original shape when subjected to external forces. Many systems inherently exhibit elastic behavior, such as underdamped oscillators and viscoelastic fluids, which introduce inertia effects into their dynamics. However, the influence of elasticity on chaotic systems remains largely unexplored, with previous investigations mainly focusing on specific physically motivated problems [7].
In this paper, we investigate the role of elasticity in chaotic systems by studying different models borrowed from mechanics [2], ecology [8,9], immunology [10,11], and rheology [12]. Our primary focus is on the impact of elasticity on the chaotic behavior of these systems. We perform a detailed numerical study of a viscoelastic fluid confined in a closed-loop thermosyphon with general heat flux and show the impact of the viscoelastic parameter on the model’s chaotic behavior. We also extend previous work on asymptotic behavior for this model, considering a viscoelastic fluid [6]. There are alternative approaches to model reduction. In particular, in [7], the so-called sparse identification of nonlinear dynamics (SINDy) method is applied to study the chaotic dynamics of a Newtonian fluid in a closed-loop thermosyphon. SINDy aims to discover the underlying dynamics of a system informed by data. While these methods allow extending the range of validity of the reduced-dimensionality model, we aim to provide a top-down derivation based on constitutive equations from different fields.
In particular, we consider a non-Newtonian heat transfer law as an alternative to the widely used Newton’s linear cooling model used in [13], where we also consider a viscoelastic fluid. The heat transfer law used in this paper is the main difference of this work with respect to previous work [13]. However, it is used in the thermosyphon model for a Newtonian fluid (without viscoelasticity), as [14,15]. The novelty of our work is twofold: firstly, we show—using different examples—how four wildly dissimilar fields present strikingly similar chaotic behavior through a detailed derivation of effective equations displaying chaos. Secondly, all of them contain elastic terms (although in the case of immunology, in a metaphorical way) that, despite the traditional roles attributed to them, do not necessarily provide periodic oscillations but rather create new and unpredictable transitions from chaotic to non-chaotic behavior.
The paper is organized as follows. In Section 2, we introduce the four models in their original context and show that they share similar underlying equations, which we denote as the elastic-Lorenz equations, generically. In Section 3, we numerically integrate the equations from Section 2 to illustrate the nontrivial effect of elasticity. We compare and contrast the results obtained from the three models derived from the thermosyphon equations and discuss the implications of our findings. Finally, in Section 4, we summarize our main results and conclusions. We also present appendices for the mathematical derivations.

2. Notations and Mathematical Formulation of the Problems

2.1. Mechanics: A Linear-Jerk-Controlled Chaotic Waterwheel

Inspired by the classic chaotic waterwheel built by Willem Malkus and Lou Howard at MIT in the 1970s [2], we introduce a simple generalization based on the concept of smooth control over the jerk [16], used in mechatronics in rigid robotic manipulators [17]. The jerk is the time derivative of the acceleration. Here, we take the basic chaotic waterwheel depicted in Figure 1 and add a linear-jerk controller: an actuator that works against the increase of the jerk. This sort of jerk control is used in elevators to improve the passengers’ comfort during take-up and stopping.
The waterwheel’s mechanism is simple: water is introduced from the top of the system at a consistent rate. If the flow rate is too slow, the top cups will not receive enough water to overcome friction, and the wheel will remain stationary. However, when the inflow rate increases, the top cup becomes heavy enough to set the wheel in motion. The wheel eventually reaches a stable rotation in either direction because both directions are equally likely, depending on the initial conditions. Interestingly, a further increase in flow rate causes the stable rotation to become chaotic. During this chaotic motion, the wheel rotates in one direction for a few turns before some cups become too full, slowing down or reversing the wheel’s direction. The wheel changes direction erratically, and observers have been known to place small bets on its subsequent turn after a minute.
The equations of motion for the original waterwheel are discussed in Ref. [2] and are derived from two basic principles:
I ω = ν ω damping torque + g r 0 2 π m ( θ , t ) sin θ d θ gravitational torque ( Rigid - body motion of the wheel )
m t = Q water inflow K m water loss ω m θ water transport ( Conservation of water ) ,
where
  • θ = angle in the lab frame.
  • ω ( t ) = θ , the angular velocity of the wheel.
  • α ( t ) = ω , the angular acceleration of the wheel.
  • m ( θ , t ) = mass distribution of water at (angular) position θ .
  • Q ( θ ) = water inflow rate.
  • r = the wheel’s radius.
  • K = the amount of mass lost.
  • ν = the viscous damping coefficient.
  • I = moment of inertia of the wheel.
As mentioned above, many control systems introduced a jerk-reduction feedback control. In our case, we introduce an additional mechanism into (1),
Controlling torque : λ ω λ α ,
where ω is the angular jerk of the waterwheel. Due to the circular symmetry of the wheel, we can expand
m ( θ , t ) = n = 0 a n ( t ) sin n θ + b n ( t ) cos n θ .
and
Q ( θ ) = n = 0 q n cos n θ .
Following the same derivation using amplitude equations, as in Appendix A (similar to the derivation in Ref. [2]), we arrive at:
α = ν ω + π g r a 1 I α / λ ω = α a 1 = ω b 1 K a 1 b 1 = ω a 1 K b 1 + q 1 .

2.2. Immunology: IL7 Receptor Signaling in T Cells

T-cells play a crucial role in the adaptive immune response of jawed vertebrates. The regulation of T-cell homeostasis is achieved through the secretion and degradation of cytokines, including interleukin 7 (IL7) [18]. Recent studies have suggested that T-cells exhibit elastic-type behavior, which could explain their population-level regulation [19]. This elasticity is a macroscopic property that emerges from the coordinated behavior of individuals rather than from individual agents themselves.
We aim to model the population of T-cells, denoted as x, and the IL7 cytokine, denoted as h ( t ) , by connecting the dynamics at the receptor level [18] with those at the population level [18,19,20]. Specifically, we propose the following equations:
x = k x c x + λ y ,
y = φ x μ x z δ y ,
z = ν x y β z ,
Here, k and c represent the T-cell population’s damping coefficient and elastic constant. The parameter λ measures the force exerted on the population per unit of interleukin y. The interleukin production rate is φ per cell, while the interleukin consumption rate is μ per cell per IL7 receptor, ρ . The clearance rate of IL7 due to absorption by other cells (rather than T-cells or dilution in the lymphatic system) is represented by δ . The ligand-induced expression of IL7 receptors per cell is denoted as ν , and the downregulation rate of IL7 receptors at the cell level is represented by β [18].

2.3. Ecology: The Evolution of Coexisting Populations

Consider a population of four species in an ecological niche. Mathematically, their dynamics can be modeled using the generalized Lotka–Volterra (gLV) model [8] (or generalizations [9]). According to gLV, the time course of each population follows the first-order system of non-linear equations
x i = r i x i + j = 1 4 β i j x i x j , i = 1 4 .
In this case, it is often useful to aggregate species based on their role in the ecosystem (predator, prey, …). Also, the relevance of each species is calibrated in proportion, b i j , to their diversity. Thus, one can define a meta-population as follows:
x i = i = 1 4 a i j y j , and its inverse ,
y j = i = 1 4 b j i x i ,
where the matrix [ b i j ] is the inverse of [ a i j ] , namely,
j = 1 4 a i j b j k = δ i k ,
where δ i k is the Kronecker delta function.
Plugging these definitions into (8), we find:
y j = k = 1 4 i = 1 4 b j k r k a k i y i + k = 1 4 i = 1 4 γ j k i y k y i , j = 1 4 .
where γ j k i depends on the matrices a, b, and β . Without defining a specific problem, it is trivial to see that proper choices of b, r i , and β can be mapped into the system of Equation (20) (Lorenz), specifically, choosing γ 1 k i = 0 . This might be why chaotic oscillations have been observed in ecological systems [21].

2.4. Rheology

Our final example is the most elaborated one and deserves deeper exposition. Here, we narrow our focus on examining thermosyphon models that involve viscoelastic fluids. This presents some interesting peculiarities that could impact the behavior of a Newtonian fluid, such as water. On one hand, the dynamics exhibit memory, as they depend on past events. On the other hand, small disturbances cause the fluid to act like an elastic solid with a resonant frequency that may eventually become significant.
The simplest method of dealing with viscoelasticity is through the use of the Maxwell constitutive equation [12,22]. Despite being a simplified model, it has been demonstrated to be valid even for complex fluids like blood, in which red blood cells alter their behavior based on concentration and vessel geometry [23]. In this model, both Newton’s law of viscosity and Hooke’s law of elasticity are generalized and supplemented by an evolution equation for the stress tensor, Σ . The stress tensor is involved in the conservation of momentum equation.
ρ v t + v · v = p + · Σ + ρ g ,
Here, p represents the hydrostatic pressure, ρ stands for the material density, and g denotes the acceleration caused by gravity ( g = 9.81 m/s 2 ). Additionally, the continuity equation is introduced as a complement to the aforementioned equation, assuming the material to be incompressible (which holds sufficiently true for liquids):
· v = 0 .
In the case of a thin section thermosyphon with a Maxwell fluid, the stress tensor Σ simplifies to a single independent component, denoted as σ ˜ . This component follows the evolution governed by the equation:
μ E σ ˜ t + σ ˜ = μ γ ˙ ,
Here, μ represents the viscosity of the fluid, E corresponds to Young’s modulus, and γ ˙ is the only non-zero element of the shear strain rate tensor Γ ˙ (which describes the rate of fluid deformation), defined as
Γ ˙ = v + v T .
Under steady flow conditions, Equation (14) reduces to Newton’s law, leading to the well-known Navier–Stokes equation expressed in Equation (12). Additionally, for short time scales when we expect an impulsive response from a resting state, Equation (14) simplifies to Hooke’s law of elasticity.
In order to elucidate the coupling between fluid and thermodynamic mechanisms discussed in Section 1, we shall establish a mathematical framework that encompasses the variables of interest: the velocity v, the temperature distribution T within the fluid, and the concentration S of the solute as it circulates through the loop. The behavior of the solute concentration is governed by transport equations.
ρ T t + v · T = · J T + K T
ρ S t + v · S = · J S + K S
In the given context, J T , S denote the flows or currents of heat and solute concentration, while K T , S represent external sources of temperature or solute concentration. Generally, these terms are determined by constitutive equations that involve linear combinations of gradients of other variables. For example, according to Fourier’s law of heat conduction, J T = ν T ; thus, Equation (15) yields the familiar Laplacian term ν 2 T .
To proceed, we follow the same procedure as in Ref. [24], also summarized in Appendix A, and we arrive at the following equation, which generalizes Equation (2) for a viscoelastic fluid.
ρ μ E d 2 v d t 2 s + ρ d v d t s = ρ g cos ϕ + F s p x + μ E 2 p x t ,
where the term F s represents the component of the wall law force that acts along the pipe (more details on this function are provided below).
Additional physical assumptions, constitutive equations, and nondimensionalization lead to our main system of equations (see Appendix A):
ε d 2 v d t 2 + d v d t + G ( v ) v = ( T ( t , x ) S ( t , x ) ) f ( x ) d x , v ( 0 ) = v 0 , d v d t ( 0 ) = w 0 T t + v T x = h ( t , x ) + ν 2 T x 2 , T ( 0 , x ) = T 0 ( x ) S t + v S x = c 2 S x 2 b 2 T x 2 , S ( 0 , x ) = S 0 ( x )
The expression h ( t , x ) + ν 2 T x 2 denotes the heat transfer law along the loop, where h ( t , x ) is a known function (refer to [25,26]). As stated in the introduction, this study offers an alternative approach to the widely used Newton’s linear cooling model discussed in [13]. Our specific focus lies in investigating a non-Newtonian heat transfer law that takes into account the presence of a viscoelastic fluid.
While previous works, such as [14,15], have employed this heat transfer law in modeling thermosyphons with Newtonian fluids, they did not account for the viscoelastic effects. Hence, our current investigation merges both extensions, addressing the classical thermosyphon problem in a comprehensive manner.
In Appendix A, we provide a detailed analytical study of the system, including the well-posedness of the ODE/PDE model, the existence of a global attractor, and the construction of an inertial manifold. The inertial manifold enables us to describe the system’s asymptotic behavior via a unified and explicitly reduced system of ODEs.

3. Numerical Experiments

This section presents numerical experiments using MATLAB to integrate the differential equations numerically. In particular, we use the solver ode15s for stiff systems.
As shown in Appendix A, we derive the system of equations
ε w ˙ = 2 a 1 2 d 1 G ( v ) v w v ˙ = w a 1 ˙ = A ν 4 π 2 k 2 a 1 + 2 π k v a 2 a 2 ˙ = B 2 π v a 1 ν 4 π 2 a 2 d 1 ˙ = c 4 π 2 d 1 + b 4 π 2 a 1 + v 2 π d 2 d 2 ˙ = c 4 π 2 d 2 + b 4 π 2 a 1 v 2 π d 1
where ( a 1 ( t ) , a 2 ( t ) ) describe the temperature T ( t , x ) , ( d 1 ( t ) , d 2 ( t ) ) describe the solute concentrations S ( t , x ) and w ( t ) , a n d v ( t ) denotes the acceleration and the velocity of the fluid.
We prove in Appendix A.4 that the asymptotic behavior of the temperature T ( t , x ) and the solute concentration S ( t , x ) are given by the coefficients of Fourier expansions a k ( t ) and d k ( t ) , respectively, with k = ± 1 for a circular circuit. Moreover, the equation for k = 1 is the conjugate of k = 1 , and this is enough to consider k = 1 . Thus, if b 1 is the coefficient of the Fourier expansions of the heat flux, we obtain Equation (A72).
Finally, after of change of variables (see Remark A4), in order to simplify the notations, we can omit the subscripts since we only have one coefficient; a ( t ) and d ( t ) denote the temperature and solute concentrations, respectively. We consider the real and imaginary parts of these coefficients, given by
a ( t ) = a 1 ( t ) + i a 2 ( t ) , d ( t ) = d 1 ( t ) + i d 2 ( t ) and b 1 = A + i B .
Notice that if we assume the viscoelastic parameter ε , and the Soret diffusion coefficients b and c are zero, we recover the model analyzed in [14] for a Newtonian fluid. Moreover, if we include the Soret effect, but ε = 0 , our model reduces to that in Ref. [15].
Hereafter, we consider the full model (19) (that we refer to as Model M3), a simplified model without the Soret effect (Model M1), and the viscoelastic Lorenz model (Model M2). Although M1 and M2 are intrinsically equivalent, we study this using the traditional parametrization to preserve the familiarity with the Lorenz model in the chaos literature.
Also, as the general model (19) supersedes the other two, and it has a ground physical interpretation, we entitle Section 3.1 and Section 3.3 explicitly in terms of the physical assumptions underlying each model.

3.1. M1 ( ( v ( t ) , T ( t , x ) ) : Thermodiffusion and Viscoelasticity in a Closed-Loop Thermosyphon without Solute S(t,x)

In this section, we consider a fluid without solute and we analyze and discuss the behavior of the v ( t ) and the temperature T ( t , x ) ( a 1 ( t ) , and a 2 ( t ) ) of the fluid, given by the solutions of the following system:
ε w ˙ = 2 a 1 G ( v ) v w v ˙ = w a 1 ˙ = A ν 4 π 2 k 2 a 1 + 2 π k v a 2 a 2 ˙ = B 2 π v a 1 ν 4 π 2 a 2
We describe the results of numerical experiments that were obtained, considering various parameters of the viscoelastic coefficients.
In the absence of viscoelasticity, the Lorenz-type behavior of the system was observed, see [14]. Therefore, for the numerical simulations, we take the same values of the parameters and the initial conditions as in [14] to explore the effect of viscoelasticity on the chaotic behavior of the system. That is, A = 0 , B = 30 , and ν = 0.025 , and the initial conditions, are fixed as w ( 0 ) = , v ( 0 ) = 0 ,   a 1 ( 0 ) = 0.1 , and a 2 ( 0 ) = 1 .
We study the behavior of the system upon varying the viscoelastic coefficient ε . To assess the effect of viscoelasticity on the system response, we generated a bifurcation diagram to visualize the dynamic behavior of the system. The parameter ε serves as a control variable that represents the degree of viscoelasticity in the system. By varying ε and observing the corresponding peaks of the fluid velocity at large times, we can gain valuable insights into the system’s response under different viscoelastic conditions, see Figure 2.
Through the analysis of the bifurcation diagram, we can observe that if the parameter value is small, chaos remains as expected. Then, between epsilon 10 2 and 10 1 , there seems to be a region of periodicity that transitions into stable behavior between 10 1 and 10. However, the surprising finding is that near 10, the velocity behavior becomes more complex.
Given the multidimensional nature of the system, we plot the velocity in temporal graphs or the impact of the relevant modes of temperature on velocity in phase-space diagrams. During our analysis, we carefully observe significant qualitative changes in the system’s behavior. For low enough values of ε , we recover, as expected, Lorenz-type behaviors observed in [14]; see Figure 3. In fact, the system presents chaotic behavior when the value of the viscoelastic component is between 10 7 and 10 3 .
By increasing the viscoelastic coefficient, a significant alteration in the system’s behavior becomes apparent. As the coefficient increases within the range of 10 3 to 10 2 , the system’s behavior gradually undergoes a transformation, eventually transitioning into a periodic pattern. Observe the emergence of periodic behavior in the velocity at 10 2 , Figure 4, and observe that the phase diagram has simplified for ε = 10 2.4 , Figure 5.
Beyond 10 2 , the system begins to stabilize. For instance, for ε = 10 1 , after a transition, the system has stable behavior. Notice that, in this case, equilibrium is reached. This behavior holds until ε = 10 . However, it is important to note that the transition becomes more complex as ε is larger; see Figure 6.
Corroborating what was observed in the phase diagram, as the value of ε approaches 10, the system shows a tendency to exhibit more complex behaviors. This observation becomes evident when ε varies between 10 and 12. Surprisingly, within this range, the system reverts back to displaying chaotic behavior; see Figure 7.
Finally, when ε goes beyond 12, the system stabilizes. First, quasi-periodic behavior is observed, and then, as ε is large enough, periodic behavior is found; see Figure 8.
In summary, the effect of the viscoelastic parameter ε on chaos is far from trivial. For large values, it tends to remove chaos but, as shown in Figure 2, there are regions where chaos reappears. To study this more systematically, in the discussion section, we display the maximum value of the Lyapunov exponent for different values of ε .

3.2. M2 General Lorenz System: Considering the Viscoelastic Fluid

For the thermosyphon model without viscoelasticity, ε = 0 in (20), equivalence with the classical Lorenz equations was shown, see [14]. The authors performed some changes of variables to prove that the three-dimensional model of the thermosyphon corresponded exactly to the Lorenz system in the case of linear friction. Analogously, performing the same changes of variables as in [14], we obtain from (20) the following system:
ε w = σ ( y x ) w x = w y = x · ( ρ z ) y z = x y β z .
Observe that if the parameter ε is zero in the previous model, we recover the classical Lorenz system:
x = σ ( y x ) y = x · ( ρ z ) y z = x y β z .
We are interested in analyzing how the parameter ε variation affects the Lorenz system’s chaotic regime. Therefore, taking the typical parameters that cause chaos in the classical Lorenz system, we analyze and discuss the system’s behavior for different parameters of the viscoelastic coefficients.
Therefore, we consider σ = 10 , β = 8 3 , and ρ = 28 , and initial conditions are fixed as w ( 0 ) = 0 , x ( 0 ) = 1 , y ( 0 ) = 1 and z ( 0 ) = 0 .
Following the analysis conducted in the previous model, we first present the bifurcation diagram to understand the system’s behavior, as we vary the parameter ε , see Figure 9.
Similar to what occurred in the previous model, starting from chaos, when the viscoelasticity is sufficiently small, chaotic behavior persists. However, when ε reaches a value of 10 1 , a predictable behavior emerges. It is worth noting that near ε = 1 , small regions with more complex behaviors appear, which require further detailed analyses.
To accomplish this, we will plot the phase diagram of ( x , y , z ) for notable parameter values where the behavior changes. The chaotic behavior of the system is observed for small enough ε . Specifically, the classical Lorenz-type behavior remains from 10 7 to 10 2 ; see Figure 10.
At values below ε = 0.1 , the system exhibits chaotic behavior. However, when ε = 0.1 is reached, the system transitions towards periodic behavior. Remarkably, this periodic behavior persists until ε = 1 , but at ε = 0.1 , the system begins to change from periodic behavior. The system preserves this periodic behavior until ε = 1 ; see Figure 11.
Similar to the previous model, we have identified parameter values of epsilon that induce chaos from the initial periodic behavior of the system. In particular, for parameter ε values between 1.2 and 1.5 , the system reverts back to chaotic behavior Figure 12.
After ε = 1.5 the behavior changes drastically again. For ε = 1.6 , the solutions present a periodic behavior. Finally, the behavior becomes stable for larger values of the parameter ε .
As mentioned above, the Lorenz model is a paradigmatic case of a simple model displaying a route to chaos. To illustrate this, we show four bifurcation diagrams for the Lorenz parameter ρ in Figure 13. Note that, for small values of the viscoelastic parameter ε , the two bifurcation diagrams are almost identical (so the chaotic patterns are preserved even though ε 0 is a singular perturbation). Interestingly, larger values stretch the bifurcation pattern and create new stability gaps.
This similarity between bifurcation diagrams is best quantified in Figure 18 (red line). For values of ε 10 1.4 , the Lyapunov exponent is the same as the original Lorenz system’s. It is also worth noting that the appearance of chaos is intermittent for larger values, as seen in Figure 9.

3.3. M3 ( v ( t ) , T ( t , x ) , S ( t , x ) ) : General Model (19)

In this subsection, we analyze the influence of viscoelasticity in the general case where a binary fluid is considered, this is when we consider a solute (as water and antifreeze), and we study the evolution of the velocity v ( t ) , the temperature T ( t , x ) of the fluid (given by a 1 ( t ) , a 2 ( t ) ) and the solute concentrated S ( t , x ) (given by d 1 ( t ) , d 2 ( t ) ) , see (19).
We are going to proceed as in the previous cases. We will assume the values calculated [15] in order to obtain chaos, and then, we discuss the behavior of the new system depending on the value of the viscoelasticity ε . Therefore, we take A = 0 , B = 30 , ν = 0.025 , c = 0.001 and b = 10 3 . Moreover, we consider the initial conditions w ( 0 ) = 0 , v ( 0 ) = 0 , a 1 ( 0 ) = 0.1 , a 2 ( 0 ) = 1 , d 1 ( 0 ) = 0.01 , d 2 ( 0 ) = 1 .
First, to understand the system’s behavior globally, we plotted the bifurcation velocity diagram for large time values; see Figure 14.
Repeating the pattern observed in previous models, for small parameter values, chaotic behavior is preserved, and as the parameter representing viscoelasticity increases, chaos disappears. However, what is interesting is that for values greater than 10 in the diagram, more complex behaviors are once again observed.
To further validate the conclusions derived from the bifurcation diagram, we will graph the velocity and the temperature phase diagram for specific values of the parameter epsilon that are deemed relevant. First, we confirm the intuitive idea that while the viscoelastic coefficient is small, the behavior remains chaotic. In particular, from ε = 10 7 to ε = 10 2 chaotic behavior of the system is observed; see Figure 15.
However, once ε exceeds 10 2 , the system stabilizes rapidly, transitioning from chaotic behavior to equilibrium. Specifically, from ε = 10 2 to ε = 10 1 , the system demonstrates periodic behavior; see Figure 16. At ε = 10 1.2 , the velocity initially undergoes a complex transition but eventually converges to equilibrium. This stabilizing effect continues to manifest with varying degrees of complexity until reaching ε = 10 .
For values close to (but greater than) ε = 10 , the behavior of the solutions becomes more complex. First, we obtain periodic behavior, ε = 12 , then quasiperiodic behavior, and chaos from ε = 14 ; see Figure 17. Beyond ε = 15 , the system tends to stabilize. Sometimes the system presents periodic behavior, ε = 15 or ε = 20 , but other times, it converges to equilibrium. For instance, ε = 21 or 25. As discussed before, the maximum Lyapunov exponent allows to explain this non-monotonic effect of the viscoelastic parameters, as shown in Figure 18.

4. Discussion

Through our study of three chaotic systems, it is clear that the interaction of viscoelasticity can lead to intricate nonlinear and chaotic behaviors. As a result, drawing conclusions that are not specific to the studied system can be challenging. However, our article provides general insights into all three models based on a thorough analysis of numerical simulations.
For example, while chaos has been observed in ecological systems, our research reveals that this is not contradictory but rather a natural outcome of the Lotka–Volterra model’s inherent structure.
Traditionally, elasticity is linked to phenomena such as oscillations or resonance. However, our study demonstrates that the interplay between elasticity and other model features is more complex, resulting in complex interactions that preserve, modify, and generate novel chaotic regimes.
The intricacies resulting from the interaction of viscoelasticity are evident in the analysis of the three chaotic systems, see Table 1. As such, generalizing results beyond the specific studied systems can pose challenges. Nevertheless, our article provides comprehensive insights into all three models through detailed numerical simulations.
To summarize the behavior of models M1–M3, Figure 18 showcases the maximum Lyapunov exponent, which exhibits a fascinating pattern in the three systems as the viscoelasticity-affecting parameter increases. As viscoelasticity increases, we observe transitions from chaos to predictable long-term behaviors, like equilibria, periodicity, or quasi-periodicity. However, increasing the viscoelastic parameter can reintroduce chaos, highlighting the captivating nature of these systems, wherein changes in the viscoelastic parameter lead to transitions between ordered and chaotic dynamics.

5. Conclusions

This research delves into the role of elasticity in chaotic systems by examining four different models from various fields, such as mechanics, immunology, ecology, and rheology. Through numerical simulations, we discovered that elasticity has a significant impact on the chaotic behavior of these systems, which we refer to as elastic-Lorenz equations. Our study provides insight into the importance of elasticity in chaotic systems and sparks new areas of research.
The implications of our findings are far-reaching. For example, in immunology, where chaos plays a significant role in the behavior of the immune system, understanding the effects of elasticity on chaotic dynamics could be crucial in treating diseases. Our research also has potential applications in rheology, which studies the flow of matter and is vital in many industrial processes. Additionally, our discoveries shed light on complex systems, such as financial markets and ecosystems, where chaotic behavior is prevalent.
Furthermore, our research suggests potential extensions. We used the basic Maxwell model of viscoelasticity, but exploring linear and non-linear generalizations of the constitutive equation could reveal new and exciting phenomena. Another potential extension is the use of time-fractional models, which have successfully modeled viscoelastic materials. As our elasticity-driven models contain both the first and second derivatives of velocity, we hope that further developments in the field will result from our findings.

Author Contributions

Á.J.-C., M.C. and M.V.-P. conceived the research and wrote the manuscript. Á.J.-C. derived the inertial manifold equations for the thermosyphon model. M.V.-P. performed the numerical simulations. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially supported by grants PID2019-106339GB-I00 funded by MCIN/ AEI/10.13039/501100011033, Spain, and PID2019-103860GB-I00 funded by MCIN/AEI/Spain; by Gr5/08 Grupo 920894 BSCH-UCM from Grupo de Investigación CADEDIF and Grupo de Dinámica No Lineal (U.P. Comillas), Spain.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Mathematical Derivation of the Thermosyphon Model

For the sake of completeness, we summarize succinctly the method described in Ref. [24].
(i)
Starting from Equations (12), (13), (15) and (16), we project them along a unitary vector tangent of the pipe (hereafter, we will refer to the coordinate along that tangent simply as x).
(ii)
Given the sufficiently small size of the thermosyphon’s cross-sectional area, denoted as A, we can safely disregard variations in the transverse direction of the pipe, ensuring that the process of averaging does not result in any appreciable loss of precision.
(iii)
Finally, since the fluid is incompressible, the average velocity remains constant along the entire length of the pipe.
To emphasize the role of elasticity, we differentiate Equation (12) to eliminate the stress component σ ˜ (after multiplying by μ / E ). Thus,
ρ μ E 2 v t 2 + μ E p t + ρ v t + p ρ g = μ v + non - linear terms
Based on Refs. [24,27], it is assumed that the non-linear terms and viscous term, 2 v , can be expressed through the wall law force, F , which is dependent on the average velocity (see details below).
According to the incompressibility hypothesis stated in Equation (13), the average velocity v remains constant throughout the thermosyphon and is solely dependent on time, t, rather than the spatial coordinate along the pipe, x. Additionally, due to the projection and averaging procedures, temperature, T, and solute concentration, S, are only functions of time, t, and position along the loop, x. By projecting in the direction of the pipe, we obtain Equation (17).
In Equations (12)–(14), the velocity is linked with temperature and diffusion through the body force ρ g , where the density ρ depends on temperature (hotter fluids have lower density) and solute concentration (higher solute concentrations lead to higher densities). Mathematically, it is possible to assume a linear correlation between density and temperature. Similarly, in this work, we also assume a linear relationship between density and solute concentration. Thus, in the end,
ρ = ρ 0 ( 1 α T T + α S S ) ,
with α T and α S , being the so-called dilatation and thermophoretic coefficients, respectively.
In this work, we use the Boussinesq approximation [24,27], which assumes that the density, which is proportional to the dynamical terms d 2 v / d t 2 and d v / d t , remains approximately constant and equal to ρ 0 . It also assumes that the effects of temperature and solute concentration on inertia can be neglected. Therefore, by integrating Equation (17) over the length of the pipe L, resulting in zero integrated pressure gradients and a constant gravitational term ρ 0 g , we obtain the following equation:
ρ μ E d 2 v d t 2 + ρ d v d t + 0 L F s d x = g 0 L ( α T T ( t , x ) α S S ( t , x ) ) f ( x ) d x ,
Here, the function f ( x ) denotes a mathematical expression that relates cos ϕ to the coordinate x. Following Rodriguez-Bernal et al.’s study in 1995, we make the assumption that the local dependence of the wall law projection on velocity v can be expressed as 0 L F s d x G ( v ) v . The function G ( v ) exhibits specific general characteristics that are elucidated in the subsequent discussion.
To complete the closure of Equations (A3) and (15), it is necessary to supplement them with appropriate “constitutive equations” governing the behavior of temperature and solute concentrations. The primary mechanisms considered for temperature evolution are as follows:
  • Convection accounted for in the second term on the left-hand side of Equation (15).
  • Thermal diffusion (conduction), expressed as J T = ν T . Thus, from Equation (15), · J T = ν 2 T x 2 .
Similarly, for solute concentration, the following physical mechanisms are taken into consideration:
  • Convection, incorporated in the second term on the left-hand side of Equation (16).
  • Solute diffusion, described by J S , 1 = c S .
  • Onsager coupling (thermodiffusion or the Soret effect), represented by J S , 2 = b T .
Hence, from Equation (16), · J S , 1 + J S , 2 = c 2 S x 2 b 2 T x 2 .
Finally, we proceed to non-dimensionalize the resulting equations, whereby lengths are scaled by the loop length (with the entire loop represented by the interval [ 0 , 1 ] ), temperatures are scaled by ( g α T ) 1 , and solute concentration by ( g α S ) 1 , to arrive at the main focus of our work, summarized in the system of Equation 18.

Appendix A.1. Well-Posedness and Boundedness: Global Attractor

For completeness, we include the details of well-posedness and boundedness for the present model as it requires specific work, although the techniques used are similar to those in Refs. [28,29].
In this context, we enforce the condition that the thermal diffusivity ν 0 maintains physical consistency. The system of Equation (18) represents an ODE/PDE system governing the evolution of the observables v, T, and S. The solute concentration S ( t , x ) adheres to a solute diffusivity c > 0 and a positive Soret coefficient b, similar to previous models for binary fluids, like in [6,15,30].
In subsequent discussions, we utilize the notation · d x = 0 1 · d x for integration along a closed circuit path, which can be associated with integration over periodic functions with a period of 1. The function f in the equation characterizes the loop’s geometry and gravitational forces [24,27], where f ( x ) d x = 0 , since the integral of f ( x ) = cos ϕ ( x ) over a closed loop yields zero.
The parameter ε in Equation (18) represents the non-dimensional form of μ / E , with time units. Roughly, ε determines the (non-dimensional) timescale for the material to transition from an elastic to a fluid-like behavior.
We assume that the friction law G ( v ) at the inner wall of the loop is positive and remains bounded away from zero. Previous works have considered G ( v ) as a constant for linear friction [27] (Stokes flow) or | v | for quadratic friction [25,26]. A more general form of G ( v ) is given by g ˜ ( R e ) | v | , where R e = ρ v L / μ is the Reynolds number, and g ˜ is a function of R e . Here, we consider a general function of the velocity, assumed to be large [31,32]. The functions G, f, and l incorporate relevant physical constants of the model, such as the cross-sectional area D, the length of the loop L, and the Prandtl, Rayleigh, or Reynolds numbers. We require G and l to be continuous functions satisfying G ( v ) G 0 > 0 and l ( v ) l 0 > 0 , where G 0 and l 0 are positive constants. Additionally, we impose further constraints on G in Equation (A23), mandating G to be a nonlinear friction function, such that there exists a constant h 0 0 , satisfying the following:
lim sup s | G ( s ) | G ( s ) = 0 and lim sup s | s G ( s ) | G ( s ) h 0 .
We will introduce some function spaces that will be used to study the existence of solutions of (18). Let Ω = ( 0 , 1 ) and consider the spaces
L p e r 2 ( Ω ) = { u L l o c 2 ( I R ) , u ( x + 1 ) = u ( x ) a . e . x I R } , H p e r m ( Ω ) = H l o c m ( I R ) L p e r 2 ( Ω )
where m I N { 0 } , and u L l o c 2 ( I R ) (or H l o c m ( I R ) ) if for every open set ω I R one has u L l o c 2 ( ω ) (or H l o c m ( ω ) , respectively). Finally, we consider functions with zero average, and we denote by the following:
L ˙ p e r 2 ( 0 , 1 ) = { u L l o c 2 ( I R ) , u ( x + 1 ) = u ( x ) a . e . , u ( x ) d x = 0 } ,
H ˙ p e r m ( 0 , 1 ) = H l o c m ( I R ) L ˙ p e r 2 ( 0 , 1 ) .

Appendix A.2. Existence and Uniqueness of Solutions

In this section, we demonstrate the presence and singularity of solutions for the thermosyphon model (18), where f and h belong to L ˙ p e r 2 ( 0 , 1 ) , T 0 is in H ˙ p e r 1 ( 0 , 1 ) , and S 0 is in L ˙ p e r 2 ( 0 , 1 ) . Here, L ˙ p e r 2 ( 0 , 1 ) and H ˙ p e r 1 ( 0 , 1 ) are defined by (A5). It should be noted that the term dot indicates functions that have zero averages and not the time derivatives of the functions. To begin with, we choose this framework by observing that for ν > 0 , integrating the temperature equation along the loop, while considering the periodicity of T, we have the following:
T x d x = 2 T x 2 d x = 0 , and d d t T ( t , x ) d x = ( h ( t , x ) .
Next, if h ( t , x ) = h ( x ) , we have T ( t , x ) d x = T 0 ( x ) d x + t h ( x ) d x . Therefore, the temperature is unbounded, as t , unless h ( x ) d x = 0 . However, taking τ ( t , x ) = T ( t , x ) T ( t , x ) d x and h = h h ( x ) d x reduces to the case T ( t , x ) d x = T 0 ( x ) d x + t h ( x ) d x = 0 , since τ ( t , x ) would satisfy
τ t + v τ x = h * ( t , x ) + ν 2 τ x 2 , τ ( 0 , x ) = τ 0 ( x ) = T 0 ( x ) T 0 ( x ) d x .
Additionally, if we integrate the equation for the solute concentration along the loop and take into account the periodicity of S, we obtain the conditions S x d x = 2 S x 2 d x = 0 and d d t ( S ( t , x ) d x ) = 0 . Since S ( t , x ) d x is constant, it means that the solute S ( t , x ) = S 0 ( x ) for all t.
Let σ ( t , x ) = S ( t , x ) S 0 ( x ) d x (not to be confused with the stress component σ ˜ ). Then, from the third equation of system (18), σ satisfies the equation:
σ t + v σ x = c 2 σ x 2 b 2 τ x 2 , σ ( 0 , x ) = σ 0 ( x ) .
Finally, since f ( x ) d x = 0 , we have ( T ( t , x ) S ( t , x ) ) f ( x ) d x = ( τ ( t , x ) σ ( t , x ) ) f ( x ) d x , and the equation for v is
ε d 2 v d t 2 + d v d t + G ( v ) v = ( τ ( t , x ) σ ( t , x ) ) f ( x ) d x , v ( 0 ) = v 0 , d v d t ( 0 ) = w 0 .
Therefore, ( v , τ , σ ) satisfies system (18) with τ 0 , σ 0 , h * replacing T 0 , S 0 , h , respectively, and f ( x ) d x = τ 0 ( x ) d x = σ 0 ( x ) d x = h * d x = 0 and T ( t , x ) d x = S ( t , x ) d x = 0 for all t 0 . From now on, we consider all the functions of system (18) to have zero average.
Furthermore, since ν , c > 0 , the operators ν A = ν 2 x 2 and c A = c 2 x 2 , along with periodic boundary conditions, are unbounded, self-adjoint operators with a compact resolvent in L p e r 2 ( 0 , 1 ) , which are positive when restricted to the space of zero average functions in L ˙ p e r 2 ( 0 , 1 ) . Moreover, the equations for the temperature T and the solute concentration S in (18) are parabolic types for ν , c > 0 .
We express system (18) as the following evolution system for acceleration, velocity, temperature, and solute concentrations:
d w d t + 1 ε w = 1 ε G ( v ) v + 1 ε ( T ( t , x ) S ( t , x ) ) f ( x ) d x , w ( 0 ) = w 0 d v d t w = 0 , v ( 0 ) = v 0 T t ν 2 T x 2 = v T x + h ( t , x ) , T ( 0 , x ) = T 0 ( x ) S t c 2 S x 2 = v S x b 2 T x 2 , S ( 0 , x ) = S 0 ( x ) .
That is:
d d t w v T S + 1 ε 0 0 0 1 0 0 0 0 0 ν 2 x 2 0 0 0 0 c 2 x 2 w v T S = F 1 ( w , v , T , S ) F 2 ( w , v , T , S ) F 3 ( w , v , T , S ) F 4 ( w , v , T , S )
with
F 1 ( w , v , T , S ) = 1 ε G ( v ) v + 1 ε ( T ( t , x ) S ( t , x ) ) f ( x ) d x ,
F 2 ( w , v , T , S ) = 0 ,
F 3 ( w , v , T , S ) = v T x + h ( t , x ) ,
F 4 ( w , v , T , S ) = v S x b 2 T x 2
and initial data w v T S ( 0 ) = w 0 v 0 T 0 S 0 .
The operator B = 1 ε 0 0 0 1 0 0 0 0 0 ν 2 x 2 0 0 0 0 c 2 x 2 is a sectorial operator in Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) with domain D ( B ) = I R 2 × H ˙ p e r 3 ( 0 , 1 ) × H ˙ p e r 2 ( 0 , 1 ) and has a compact resolvent, see Equation (A5).
Using the results and techniques in the sectorial operator of [33] to prove the existence of solutions of the system, we have Theorem A1.
Theorem A1. 
We assume that H ( r ) = r G ( r ) is locally Lipschitz, f , h ( t ) L ˙ p e r 2 ( 0 , 1 ) , G ( v ) G 0 > 0 . Then, given ( w 0 , v 0 , T 0 , S 0 ) Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) , there exists a unique solution of (18), satisfying
( w , v , T , S ) C ( [ 0 , ) , Y ) C ( 0 , , I R 2 × H ˙ p e r 3 ( 0 , 1 ) × H ˙ p e r 2 ( 0 , 1 ) ) ,
d w d t , d v d t , T t , S t C ( 0 , , I R 2 × H ˙ p e r 3 δ ( 0 , 1 ) × H ˙ p e r 2 δ ( 0 , 1 ) ) ,
for every δ > 0 . In particular, (A7) defines a nonlinear semigroup, S * ( t ) in Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) , with S * ( t ) ( w 0 , v 0 , T 0 , S 0 ) = ( w ( t ) , v ( t ) , T ( t , x ) , S ( t , x ) ) .
Proof. 
In order to prove this, we cover several steps:
Step (i). We prove the local existence and regularity. This follows easily from the variation of the constants formula as presented in [33]. In order to prove this, we write the system as (A8), and we have:
U t + B U = F ( U ) , w i t h U = w v T S ,
B = 1 ε 0 0 0 1 0 0 0 0 0 ν 2 x 2 0 0 0 0 c 2 x 2 ,
F = F 1 F 2 F 3 F 4
where the operator B is a sectorial operator in Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) with domain D ( B ) = I R 2 × H ˙ p e r 3 ( 0 , 1 ) × H ˙ p e r 2 ( 0 , 1 ) , and has a compact resolvent. In this context, operator A = 2 x 2 must be understood in the variational sense, i.e., for every T , φ H ˙ p e r 1 ( 0 , 1 ) ,
A ( T ) , φ = T x φ x d x
and L ˙ p e r 2 ( 0 , 1 ) coincides with the fractional space of exponent 1 2 , as in [33]. We denote H ˙ p e r 1 ( 0 , 1 ) as the dual space and . as the norm on the space L ˙ p e r 2 ( 0 , 1 ) . If we prove that the nonlinearity F : Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) Y 1 2 = I R 2 × L ˙ p e r 2 ( 0 , 1 ) × H ˙ p e r 1 ( 0 , 1 ) is well defined, Lipschitz, and bounded on bounded sets, we obtain the local existence for the initial data in Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) .
Using H ( v ) G ( v ) v as a locally Lipschitz function, with f , h ( t ) L ˙ p e r 2 ( 0 , 1 ) , we can follow the approach in [28], and we will prove the nonlinear terms, F 1 4 in Equations (A9)–(A12), satisfying F 1 : I R 2 × L ˙ p e r 2 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) I R , F 2 : I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) I R , F 3 : I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) L ˙ p e r 2 ( 0 , 1 ) and F 4 : I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) H ˙ p e r 1 ( 0 , 1 ) ; that is, F : Y Y 1 2 is well defined, Lipschitz, and bounded on bounded sets.
For F 3 , in this case, we have that
F 3 ( w 1 , v 1 , T 1 , S 1 ) F 3 ( w 2 , v 2 , T 2 , S 2 ) = v 1 ( T 1 ) x + v 2 ( T 2 ) x
( v 2 v 1 ) ( T 1 ) x + v 2 ( ( T 2 T 1 ) ) x
v 2 v 1 T 1 ( t ) H ˙ p e r 1 + v 2 T 2 T 1 ( t ) H ˙ p e r 1
C v 2 v 1 + T 2 T 1 ( t ) H ˙ p e r 1 .
In this way, F 3 is locally Lipschitz and bounded on bounded sets.
Applying the techniques of the variation of constants formula from [33], we obtain the unique local solution ( w , v , T , S ) C ( [ 0 , t * ] , Y ) (with a suitable t * > 0 ) of (A7), which is given by
w ( t ) = w 0 e 1 ε t 1 ε 0 t e 1 ε ( t r ) H ( r ) d r + + 1 ε 0 t ( T ( r , x ) S ( r , x ) ) f ( x ) d x e 1 ε ( t r ) d r
with H ( r ) = G ( v ( r ) ) v ( r ) .
v ( t ) = v 0 + 0 t w ( r ) d r
T ( t , x ) = e ν A t T 0 ( x ) + 0 t e ν A ( t r ) h ( x , r ) ] d r 0 t e ν A ( t r ) v ( r ) T ( r , x ) x d r ,
S ( t , x ) = e c A t S 0 ( x ) + 0 t e c A ( t r ) [ v ( r ) S x ( r ) b 2 T x 2 ( r ) ] d r .
where ( w , v , T , S ) C ( [ 0 , t * ] , Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) ) and using the results from [33], (the smoothing effect of the equations with the bootstrapping method), we obtain the regularity of solutions.
Step (ii). To prove global existence, we must show that the solutions are bounded in the Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) norm on finite time intervals, and using the nonlinearity of F , maps bounded on bounded sets, we can conclude.
Part I: We study the norm of temperature, T.
To prove that the norm of T is bounded on finite time, we multiply the equation for the temperature by T in L ˙ p e r 2 ( 0 , 1 ) . Then, integrating by parts, we obtain:
1 2 d d t T 2 + ν T x 2 = h ( t , x ) T ( t , x ) d x
since T ( t , x ) T x d x = 1 2 x ( T 2 ) d x = 0 .
By using Cauchy–Schwartz and Young inequalities, as well as the Poincaré inequality for functions of zero average, since T ( t , x ) d x = 0 , , and as π 2 is the first nonzero eigenvalue of A = 2 x 2 in L ˙ p e r 2 ( 0 , 1 ) , we obtain
1 2 d d t T 2 + ν π 2 T 2 1 4 δ h ( t ) 2 + δ T 2 ,
for every δ > 0 . Now, taking δ = ν π 2 2 , we have
d d t T 2 + ν π 2 T 2 h ( t ) 2 ν π 2 ,
and conclude that the norm of T in L ˙ p e r 2 ( 0 , 1 ) remains bounded in finite time.
Now, we will prove that the norm T x remains bounded in finite time intervals. For this, multiply the second equation of (A7) by 2 T x 2 in L ˙ p e r 2 ( 0 , 1 ) . By integrating by parts, applying the Young inequality, and taking into account that
T x 2 T x 2 d x = 1 2 ( T x ) 2 x = 0 ,
since T x is periodic, we obtain
1 2 d d t T x 2 + ν 2 T x 2 2 C δ h ( t ) 2 + δ 2 T x 2 2
for every δ > 0 and C δ = 1 4 δ . Thus, taking δ = ν 2 , and applying the Poincaré inequality for functions with zero average, we obtain
d d t T x 2 + ν π 2 T x 2 h ( t ) 2 ν ,
which proves that the norm of T in H ˙ p e r 1 ( 0 , 1 ) remains bounded in finite time.
Part II: We study the norm of solute concentration S.
In this step, we show that the norm of S in L ˙ p e r 2 ( 0 , 1 ) does not blow up in finite time. Multiplying the fourth equation of (A7) by S, integrating by parts, applying the Young inequality and again taking into account that S ( t , x ) S x d x = 1 2 S 2 x d x = 0 , since S is also periodic, we have
1 2 d d t S 2 + ( c δ ) S x 2 b 2 C δ T x 2
for every δ > 0   with   C δ = 1 4 δ . Thus, taking δ = c 2 , and using (A19) with the Poincaré inequality for functions with zero average, we obtain
d d t S 2 + c π 2 S 2 b 2 c T x 2 k 1
with k 1 > 0 . Therefore, S ( t ) remains bounded in finite time.
Finally, since T and S are bounded in finite time, this implies that | w ( t ) | and | v ( t ) | remain bounded in finite time. Hence, we obtain a global solution in the nonlinear semigroup in Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) .

Appendix A.3. Boundedness of the Solutions: Global Attractor

In this section, we adapt the results and techniques in Refs. [15,28] with Refs. [6,14] for a fluid with one component, to prove the existence of the global attractor for a binary fluid for the semigroup defined in the space Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) , in this case, when we consider this transfer law with temperature diffusion.
To obtain the asymptotic bounds on the solutions, as t , we consider the prescribed flux h ( t , x ) , satisfying that there exists h 0 > 0 , such that
h 0 = lim sup t h ( t , x ) ,   with   h 0 = h   if   h ( t , x ) = h ( x ) .
Moreover, we consider the friction function G, as in [6,14,28], satisfying the hypotheses of the previous section; there exits a constant g 0 0 , such that
lim sup s | G ( s ) | G ( s ) = 0   and   lim sup s | s G ( s ) | G ( s ) g 0 .
Using the l’Hopital’s lemma proved in [32], we have the following lemma proved in [28].
Lemma A1. 
If we assume that G ( r ) and H ( r ) r G ( r ) satisfy the hypotheses of Theorem A1, with (A23), then:
lim sup t H ˜ ( t ) 1 ε 0 t e 1 ε ( t r ) H ˜ ( r ) d r G ˜ ( t ) H 0
with H 0 = ( 1 + g 0 ) ε being a positive constant, such that H 0 0 if ε 0 , and G ˜ ( r ) = G ( v ( r ) ) and H ˜ ( r ) = v ( r ) G ˜ ( r ) .
Remark A1. 
We note that the conditions (A23) are satisfied for all friction functions G considered in the previous works, i.e., the thermosyphon models where G is constant, or linear or quadratic laws. Moreover, the conditions (A23) are true for G ( s ) A | s | n , a s s .
Theorem A2. 
Under the above notation and hypotheses of Theorem A1, if we assume that G satisfies (A24) for a constant H 0 0 , and h ( t , x ) satisfies (A22), then
(i) 
lim sup t T ( t ) h 0 ν π 2 ,
(ii) 
lim sup t T x ( t ) h 0 ν π 2 ,
(iii) 
lim sup t S ( t ) b h 0 c ν π 2 a n d lim sup t T ( t ) S ( t ) h ν π 2 ( 1 + b c ) .
(iv) 
Moreover, if c = ν + b , then we also have
lim sup t T ( t ) S ( t ) h 0 c π 2   w i t h   h c π 2 h ν π 2 ( 1 + b c ) i f ν 2 c c + 1 .
(v) 
lim sup t | v ( t ) | ( b + c ) h 0 f c ν π 2 G 0 + H 0   a n d
i f c = ν + b , t h e n lim sup t | v ( t ) | H 0 + h 0 f c π 2 G 0 .
(vi) 
lim sup t | w ( t ) | G 0 * H 0 + 1 + G 0 * G 0 ( b + c ) h 0 f c ν π 2   a n d
i f c = ν + b , t h e n lim sup t | w ( t ) | G 0 * H 0 + 1 + G 0 * G 0 h 0 f c π 2 .
Here, G 0 * = lim sup t G ( v ( t ) when v satisfies (A29).
In particular, we have a global compact and connected attractor A in Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) .
Proof. 
(i)
From (A17), using Gronwall’s lemma, for every t t 0 , we have
T 2 h 0 2 ν 2 π 4 + T 0 2 h 0 2 ν 2 π 4 + e ν π 2 t ,
where h 0 is given by (A22), i.e.,
h 0 = lim sup t h ( t , x ) = h if h ( t , x ) = h ( x ) ,
and we have (A25).
(ii)
From (A19), and working as above, for every t t 0 ,
T x 2 h 0 2 ν 2 π 2 + T 0 x 2 h 0 2 ν 2 π 4 + e ν π 2 t
then we obtain (A26).
(iii)
Using Gronwall’s lemma, from (A21) with (A26), we have (A27).
(iv)
From system (A7), if c = ν + b , T S satisfies
( T S ) t + v ( T S ) x = h ( t , x ) + c 2 ( T S ) x 2 ,
this is the same kind of equation for T. Therefore, working as in (i), we have
T S 2 h 0 2 c 2 π 4 + T 0 S 0 2 h 0 2 π 4 + e c π 2 t ,
and we conclude (A28).
(v)
Working as ([28]), from (A7), we have
d w d t + 1 ε w = 1 ε G ( v ) v + 1 ε ( T ( t , x ) S ( t , x ) ) f ( x ) d x
and w ( t ) = d v d t satisfies
d v d s = w ( 0 ) e 1 ε s 1 ε 0 s e 1 ε ( s r ) H ˜ ( r ) d r + + 1 ε 0 s ( T ( r , x ) S ( r , x ) ) f ( x ) d x e 1 ε ( s r ) d r
where H ˜ ( r ) = H ( v ( r ) ) = v ( r ) G ( v ( r ) ) and G ˜ ( s ) = G ( v ( s ) ) . We rewrite (A35) as
d v d s + G ˜ ( s ) v = w ( 0 ) e 1 ε s + I 1 ( s ) + I 2 ( s ) ,
with
I 1 ( s ) = 1 ε 0 s ( T ( r , x ) S ( r , x ) ) f ( x ) d x e 1 ε ( s r ) d r   and I 2 ( s ) = H ˜ ( s ) 1 ε 0 s e 1 ε ( s r ) H ˜ ( r ) .
For any δ > 0 , there exits t 0 > 0 , such that δ ( s ) = w ( 0 ) e 1 ε < δ for any s t 0 , and integrating (A36) with t t 0 , we obtain
| v ( t ) | | v ( t 0 ) | e t 0 t G ˜ ( s ) d s + e t 0 t G ˜ ( s ) d s t 0 t e t 0 s G ˜ ( r ) d r ( δ + | I 1 ( s ) | + | I 2 ( s ) | )
Using L’Hopital’s lemma proved in [32], we obtain the following two results:
lim sup t e t 0 t G ˜ ( s ) d s t 0 t e t 0 s G ˜ ( r ) d r ( | I 1 ( s ) | + | I 2 ( s ) | + δ ) =
= lim sup t t 0 t e t 0 s G ˜ ( r ) d r ( | I 1 ( s ) | + | I 2 ( s ) | + δ ) d s e t 0 t G ˜ ( s ) d s lim sup t | I 1 ( t ) | + | I 2 ( t ) | + δ G ˜ ( t )
for any δ > 0 and
lim sup t | I 1 ( t ) | lim sup t 0 t e r ε | ( T ( t , x ) S ( t , x ) ) f ( x ) d x | ε e t ε lim sup t | ( T ( t , x ) S ( t , x ) ) f ( x ) d x |
with
lim sup t | ( T ( t , x ) S ( t , x ) ) f ( x ) d x | f lim sup t ( T + S ) ;
and from (A38) with (A24), we conclude for any δ ,
lim sup t | v ( t ) | f lim sup t ( T S ) G 0 + H 0 + δ .
Thus, using T S T + S with the above results, we obtain (A29).
(vi) From (A34) with Gronwall’s lemma, we have
| w ( t ) | | w ( t 0 ) | e 1 ε t + 1 ε t 0 t e 1 ε ( t r ) G ( r ) | v ( r ) | + | ( T ( r , x ) S ( r , x ) ) f ( x ) d x | d r
where G ˜ ( r ) = G ( v ( r ) ) . Consequently, for any δ > 0 , there exits t 0 , such that for any t t 0
1 ε t 0 t e 1 ε ( t r ) G ( v ( r ) ) | v ( r ) | + | ( T ( r , x ) S ( r , x ) ) f ( x ) d x | d r
[ δ + lim sup t G ( v ( t ) ) | v ( t ) | + | ( T ( t , x ) S ( t , x ) ) f ( x ) d x | ( 1 e 1 ε ( t t 0 ) )
this is
lim sup t | w ( t ) | lim sup t G ( v ( t ) ) | v ( t ) | + | ( T ( t , x ) S ( t , x ) ) f ( x ) d x | + δ ,
for any δ > 0 , and using the above results, we have (A30).
Finally, since the sectorial operator B, as defined in Theorem A1, has a compact resolvent, the rest follows from [34] [Theorems 4.2.2 and 3.4.8]. □

Appendix A.4. Asymptotic Behavior: Reduction to Finite-Dimensional Systems

Finally, in this section, we aim to analyze the long-term behavior of system (A7) by examining its asymptotic properties. To achieve this, we utilize the Fourier expansions of the functions involved in the system, namely the temperature and solute concentrations. Our focus lies on the Fourier coefficients of the loop’s geometric function, denoted as f ( x ) , and the prescribed flux across the loop wall, denoted as h = h ( x ) . These coefficients play a crucial role in studying the system’s asymptotic behavior. We establish the system’s asymptotic behavior (A7) through the appropriate Fourier coefficients associated with f and h.
To investigate the system’s asymptotic behavior, we consider the coefficients a k ( t ) for temperature and d k ( t ) for the solute concentration, where k K J represent the relevant modes. Thus, we obtain a finite system of differential equations:
d w d t = 2 R e ( k ( K J ) + a k ( t ) c k ) ε 2 R e ( k ( K J ) + d k ( t ) c k ) ε G ( v ) v ε w ε , d v d t = w , d ( a k ) d t = l ( v ) b k l ( v ) a k ( t ) 4 ν π 2 k 2 a k ( t ) 2 π k v a k ( t ) i , d ( d k ) d t = 4 c π 2 k 2 d k ( t ) + 4 b π 2 a k ( t ) 2 π k v d k ( t ) i .
Furthermore, it is worth noting that the Fourier expansion of any function g belonging to H ˙ p e r m ( 0 , 1 ) , where m 0 , can be expressed as g ( x ) = k I Z * a k e 2 π k i x , with I Z * = I Z \ { 0 } . Moreover, we have the following expression for the norm of g in H ˙ p e r m ( 0 , 1 ) :
g H ˙ p e r m ( 0 , 1 ) = ( 2 π ) m k I Z * k 2 m | a k | 2 1 2 .
Assume that T H ˙ p e r 1 ( 0 , 1 )   and   f , h , S L ˙ p e r 2 ( 0 , 1 ) are given by the following Fourier series expansions:
h ( x ) = k I Z * b k e 2 π k i x   and   f ( x ) = k I Z * c k e 2 π k i x   with   I Z * = I Z \ { 0 }
T ( t , x ) = k I Z * a k ( t ) e 2 π k i x   and   S ( t , x ) = k I Z * d k ( t ) e 2 π k i x
with the initial data T 0 H ˙ p e r 1 ( 0 , 1 ) are given by T 0 ( x ) = k I Z * a k 0 e 2 π k i x and S 0 L ˙ p e r 2 ( 0 , 1 ) is given by S 0 ( x ) = k I Z * d k 0 e 2 π k i x . Since all functions involved are real and periodic, we have a ¯ k = a k , b ¯ k = b k , c ¯ k = c k and d ¯ k = d k .
Proposition A1. 
Under the above notation and hypotheses of Theorem A1, we consider h , f L ˙ p e r 2 ( 0 , 1 ) given by (A46) and the initial data T 0 H ˙ p e r 1 ( 0 , 1 ) given by T 0 ( x ) = k I Z * a k 0 e 2 π k i x and S 0 L ˙ p e r 2 ( 0 , 1 ) given by S 0 ( x ) = k I Z * d k 0 e 2 π k i x . Let ( w , v , T , S ) be the solution of system (18) given by Theorem A1, we then have the following:
(i) 
Coefficients a k ( t ) and d k ( t ) in (A47) satisfy the following equations:
d ( a k ) d t + 2 k π v i + 4 ν k 2 π 2 a k ( t ) = b k , a k ( 0 ) = a k 0 , k I Z * d ( d k ) d t + 2 k π v i + 4 c k 2 π 2 d k ( t ) = 4 b π 2 k 2 a k ( t ) , d k ( 0 ) = d k 0 , k I Z * .
(ii) 
The equation for velocity is
ε d 2 v d t 2 + d v d t + G ( v ) v = k I Z * ( a k ( t ) d k ( t ) ) c ¯ k .
Proof. 
By using the following Fourier expansion for a generic function u ( t , x ) of u ( t , x ) = k I Z * u k ( t ) e 2 π k i x   with   I Z * = I Z \ { 0 } , the Fourier coefficients u k ( t ) satisfy the following relations:
u ( t , x ) t = k I Z * d u k d t e 2 π k i x ,
u ( t , x ) x = k I Z * 2 π k i u k ( t ) e 2 π k i x ,
2 u ( t , x ) x 2 = k I Z * 4 π 2 k 2 u k ( t ) e 2 π k i x .
Consider the model (18) and the Fourier series expansions of all functions, depending on the spatial variable x, i.e., T ( t , x ) = k I Z * a k ( t ) e 2 π k i x ,
With S ( t , x ) = k I Z * d k ( t ) e 2 π k i x given by (A47), h ( x ) = k I Z * b k e 2 π k i x and f ( x ) = k I Z * c k e 2 π k i x given by (A46), with the expansions of initial data for temperature T 0 ( x ) = k I Z * a k 0 e 2 π k i x and for solute S 0 ( x ) = k I Z * d k 0 e 2 π k i x , we can easily find that the coefficients for temperature a k ( t ) and solute concentration d k ( t ) are the solutions of (A48).
Moreover, it is sufficient to note that
( T ( t , x ) S ( t , x ) ) f ( x ) d x = k I Z * ( a k ( t ) d k ( t ) ) c ¯ k ,
since all the functions involved are real and periodic, we have for all k I Z * = I Z \ { 0 } , a ¯ k = a k , b ¯ k = b k , c ¯ k = c k   and   d ¯ k = d k ; this allows us to conclude that (18) is equivalent to infinite systems of ODEs consisting of (A48) coupled with
ε d 2 v d t 2 + d v d t + G ( v ) v = k I Z * ( a k ( t ) d k ( t ) ) c k , v ( 0 ) = v 0 , d v d t ( 0 ) = w 0 .
Remark A2. 
It is worth mentioning that system (18) is synonymous with system (A7) concerning acceleration, velocity, temperature, and solute concentration. Moreover, based on the proposition above, it can be equated to the following infinite system of ordinary differential Equation (A53):
d w d t + 1 ε w = 1 ε G ( v ) v + 1 ε k I Z * ( a k ( t ) d k ( t ) ) c k , w ( 0 ) = w 0 d v d t = w , v ( 0 ) = v 0 d ( a k ) d t + 2 k π v i + 4 ν k 2 π 2 a k ( t ) = b k , a k ( 0 ) = a k 0 , k I Z * d ( d k ) d t + 2 k π v i + 4 c k 2 π 2 d k ( t ) = 4 b k 2 π 2 a k ( t ) , d k ( 0 ) = d k 0 , k I Z * .
The system of Equation (A53) captures two significant characteristics: (i) the interaction between the modes occurs through the velocity, while diffusion acts as a linear damping term, and (ii) the evolution of temperature impacts the evolution of solute concentration.
In the subsequent analysis, we will utilize this explicit equation for the Fourier modes of temperature and solute concentrations to examine the asymptotic behavior of the system and derive explicit low-dimensional models. Similar explicit constructions were presented in various previous studies, such as [28], and by Bloch and Titi in [35] for a nonlinear beam equation, where nonlinearity arises solely from the appearance of the L 2 norm of the unknown. Stuart also provided a related construction in [36] for a nonlocal reaction–diffusion equation.
In the following section, we will establish the boundedness of these coefficients, which enhances, to some extent, the boundedness of temperature and solute concentrations. This, in turn, proves the existence of the inertial manifold for system (A7), employing similar techniques as in Refs. [6,14,15,28,29].

Appendix A.5. Inertial Manifold

We consider the general case ν > 0 with the heat transfer law across the loop wall given by a prescribed function h ( x ) , and use inertial manifold techniques, in the spirit of the non-diffusion case of [37]. In this case, the existence of an inertial manifold does not rely on the existence of significant gaps in the spectrum of the elliptic operator but on the invariance of particular sets of Fourier modes.
Proposition A2. 
Under the above notation and hypotheses of Theorem A2, with initial conditions ( w 0 , v 0 , T 0 , S 0 ) Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) , for every solution of system (A7), ( w , v , T , S ) , and for every k I Z * , and recalling the expansions
h ( x ) = k I Z * b k e 2 π k i x   a n d   f ( x ) = k I Z * c k e 2 π k i x   w i t h   I Z * = I Z \ { 0 }
T ( t , x ) = k I Z * a k ( t ) e 2 π k i x   a n d   S ( t , x ) = k I Z * d k ( t ) e 2 π k i x
with the initial data T 0 H ˙ p e r 1 ( 0 , 1 ) given by T 0 ( x ) = k I Z * a k 0 e 2 π k i x and S 0 L ˙ p e r 2 ( 0 , 1 ) given by S 0 ( x ) = k I Z * d k 0 e 2 π k i x , we have
(i) 
lim sup t | a k ( t ) | | b k | 4 π 2 ν k 2 , i . e . , lim sup t k 2 | a k ( t ) | | b k | 4 π 2 ν ,
lim sup t | d k ( t ) | b | b k | 4 π 2 ν k 2 c , i . e . , lim sup t k 2 | d k ( t ) | b | b k | 4 π 2 ν c ,
lim sup t | v ( t ) | 1 4 π 2 ν I 0 G 0 1 + b c + H 0 ,   w i t h   I 0 = k I Z * | b k | | c k | k 2
and G 0 is a positive constant, such that G ( v ) G 0 ,
lim sup t | w ( t ) | G 0 * H 0 + 1 + G 0 * G 0 1 + b c 1 4 π 2 ν I 0 ,   w i t h   G 0 * = lim sup t G ( v ( t ) ) .
(ii) 
lim sup t T ( t ) h 4 π 2 ν , lim sup t T ( t ) H ˙ p e r 1 h 2 π ν   a n d   lim sup t T ( t ) H ˙ p e r 2 h ν ,
lim sup t S ( t ) b c h 4 π 2 ν , lim sup t S ( t ) H ˙ p e r 1 b c h 2 π ν   a n d   lim sup t S ( t ) H ˙ p e r 2 b c h ν ,
lim sup t | v ( t ) | f h G 0 1 + b c 1 4 π 2 ν + H 0
lim sup t | w ( t ) | G 0 * H 0 + 1 + G 0 * G 0 1 + b c 1 4 π 2 ν f h .
In particular, if h H ˙ p e r m ( 0 , 1 ) , we have the global compact and connected attractor A [ M , M ] × [ N , N ] × C × C I R 2 × H ˙ p e r m + 2 ( 0 , 1 ) × H ˙ p e r m + 2 ( 0 , 1 ) , where M , N are the upper bounds for acceleration and velocity, as given in (A61) and (A60), respectively, and T , S C = { R ( x ) H ˙ p e r m + 2 ( 0 , 1 ) , R ( x ) = k I Z * r k e 2 π k i x , | r k | d | b k | } , where d = 1 4 π 2 ν max { 1 , b c } and A is a compact set in I R 2 × H ˙ p e r m + 2 ( 0 , 1 ) × H ˙ p e r m + 2 ( 0 , 1 ) .
Proof. 
(i) From (A48), we have
a k ( t ) = a k 0 e 4 ν π 2 k 2 t e 0 t [ 2 π k v i ] + b k 0 t e 4 ν π 2 k 2 ( t s ) e s t [ 2 π k v i ] d s
and taking into account that
| e 0 t 2 π k v i | = | e s t 2 π k v i | = 1 ,
we obtain:
| a k ( t ) | | b k | 4 ν π 2 k 2 + | a k 0 | | b k | 4 ν π 2 k 2 + e 4 ν π 2 k 2 t
and we have (A54), i.e., lim sup t | a k ( t ) | | b k | 4 ν π 2 k 2 .
From (A48), we have
d k ( t ) = d k 0 e 4 c π 2 k 2 t e 0 t 2 π k v i + 4 b π 2 k 2 0 t a k ( s ) e 4 c π 2 k 2 ( t s ) e s t 2 π k v i d s ,
therefore, from (A65), we have
| d k ( t ) | | d k 0 | e 4 c π 2 k 2 t + b c lim sup t | a k ( t ) | ( 1 e 4 c π 2 k 2 t ) .
Finally, taking into account that (A66) with (A54), we obtain
| d k ( t ) | b | b k | c 4 π 2 k 2 ν + | d k 0 | b | b k | c 4 π 2 k 2 ν + e 4 c π 2 k 2 t ,
and we have (A55), i.e., lim sup t | d k ( t ) | b | b k | c 4 π 2 k 2 ν .
From (iii) in Theorem A2, with
( T ( t , x ) S ( t , x ) ) f ( x ) d x = k I Z * a k ( t ) c k k I Z * d k ( t ) c k
and using (A54) and (A55), we have
lim sup t | ( T ( t , x ) S ( t , x ) ) f ( x ) d x | 1 + b c 1 4 π 2 ν I 0 ,
where I 0 = k I Z * | b k ( t ) | | c k ( t ) | k 2 . Hence, from (A40), we obtain (A56), namely
lim sup t | v ( t ) | I 0 G 0 1 4 π 2 ν 1 + b c + H 0
and using (A43), we obtain (A57), i.e.,
lim sup t | w ( t ) | G 0 * H 0 + 1 + G 0 * G 0 1 + b c 1 4 π 2 ν I 0 ,   with   G 0 * = lim sup t G ( v ( t ) ) .
(ii) Using Theorem A2 and taking into account Equation (A45) with I 0 k I Z * | b k ( t ) | | c k ( t ) | h f , we find that for any solution of (A7), in conjunction with (A58), (A59), (A60), and (A61), given that the sectorial operator B defined above (in Section 2.1.1) has a compact resolvent from [Theorem 4.2.2 and 3.4.8] in Ref. [30], the system has a global, compact, and connected attractor, A , in Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) .
Next, we show that A [ M , M ] × [ N , N ] × C × C , where C = { R ( x ) H ˙ p e r m + 2 ( 0 , 1 ) , R ( x ) = k I Z * r k e 2 π k i x , | r k | d | b k | } , , with d = 1 4 π 2 ν max { 1 , b c } .
From Equations (A54) and (A55), for any ( w ( t ) , v ( t ) , T ( t , x ) , S ( t , x ) ) A , we have k 2 | a k | d | b k | and k 2 | d k | d | b k | ; therefore, T H ˙ p e r m + 2 d h H ˙ p e r m and S H ˙ p e r m + 2 d h H ˙ p e r m , i.e., if h H ˙ p e r m ( 0 , 1 ) then T , S H ˙ p e r m + 2 ( 0 , 1 ) , and we have ( w ( t ) , v ( t ) , T ( t , x ) , S ( t , x ) ) [ M , M ] × [ N , N ] × C × C ; that is, we have A [ M , M ] × [ N , N ] × C × C .
Finally, we show that C is compact in H ˙ p e r m + 2 ( 0 , 1 ) .
Certainly, given any sequence { R n } in C , we can derive a subsequence denoted as { R n } , which exhibits weak convergence to a function R. Moreover, for any k I Z * , the Fourier coefficients satisfy r k n r k as n , where r k represents the kth Fourier coefficient of R. Consequently, we have k 2 | r k | d | b k | for every integer N 0 .
R n R m + 2 2 | k | = 1 N 0 | k | 2 ( m + 2 ) | r k n r k | 2 + C | k | = N 0 + 1 | k | 2 m | b k | 2 ,
since there exists a positive constant C, such that
| k | = N 0 + 1 | k | 2 m | k | 4 | r k n r k | 2 C | k | = N 0 + 1 | k | 2 m | b k | 2 ;
where . m + 2 denotes the norm in H ˙ p e r m + 2 ( 0 , 1 ) . Hence, the first term goes to zero, such as n , and the second can be made arbitrarily small, such as N 0 , since h H ˙ p e r m ( 0 , 1 ) with h m 2 = C k I Z * | k | 2 m | b k | 2 < . Consequently, R C and R n R in H ˙ p e r m + 2 ( 0 , 1 ) , and the result is proved. □
Now, we will prove that there exists an inertial manifold M (see a definition in Ref. [38]) for the semigroup S * ( t ) in the phase space Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) , i.e., a submanifold of Y , such that
(i)
S * ( t ) M M for every t 0 ,
(ii)
there exists δ > 0 , satisfying that for every bounded set, B Y , there exists C ( B ) 0 , such that d i s t ( S ( t ) , M ) C ( B ) e δ t , t 0 ; see, for example, [38,39].
Assume that h H ˙ p e r 1 ( 0 , 1 ) with
h = k K b k e 2 π k i x
with b k 0 for every k K I Z * with 0 K , since h ( x ) d x = 0 . We denote by V 1 and V 0 the closure of the subspaces of H ˙ p e r 1 ( 0 , 1 ) and L ˙ p e r 2 ( 0 , 1 ) , respectively, generated by { e 2 π k i x , k K } .
Theorem A3. 
Assume that h H ˙ p e r 1 ( 0 , 1 ) and f L ˙ p e r 2 ( 0 , 1 ) . Then, the set M = I R 2 × V 1 × V 0 is an inertial manifold for the flow of S * ( t ) ( w 0 , v 0 , T 0 , S 0 ) = ( w ( t ) , v ( t ) , T ( t ) , S ( t ) ) in the space Y = I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 2 ( 0 , 1 ) . Moreover, if K is a finite set, the dimension of M is 2 | K | + 2 , where | K | is the number of elements in K.
Proof. 
Step (i). First, we show that M is invariant. Note that if k K , then b k = 0 ; therefore, if a k 0 = 0 , from (A62), we have that a k ( t ) = 0 for every t, i.e., T ( t , x ) = k K a k ( t ) e 2 π k i x , and if d k 0 = 0 , using a k ( t ) = 0 , from (A65), we have d k ( t ) = 0 for every t, i.e., S ( t , x ) = k K d k ( t ) e 2 π k i x . Therefore, if ( w 0 , v 0 , T 0 , S 0 ) M , then ( w ( t ) , v ( t ) , T ( t ) , S ( t ) ) M for every t , i.e., M is invariant.
Step (ii). From previous assertions, ( T ( t , x ) S ( t , x ) ) f ( x ) d x = k K a k ( t ) · c k k K d k ( t ) · c k , the flow on M is given by
d w d t + 1 ε w + 1 ε G ( v ) v = 1 ε k K a k ( t ) · c k 1 ε k K d k ( t ) · c k d v d t = w d ( a k ) d t + 2 π k v i + 4 ν π 2 k 2 a k ( t ) = b k , k K d ( d k ) d t + 2 π k v i + 4 c π 2 k 2 d k ( t ) = 4 b π 2 k 2 a k ( t ) , k K a k = d k = 0 , k K
Now, we consider the following decomposition in H ˙ p e r 1 ( 0 , 1 ) , T = T 1 + T 2 , where T 1 is the projection of T on V 1 and T 2 is the projection of T on the subspace generated by { e 2 π k i x , k I Z * \ K } i.e., T 1 = k K a k e 2 π k i x and T 2 = k I Z * \ K a k e 2 π k i x = T T 1 .
Analogously, we consider the decomposition S = S 1 + S 2 in L ˙ p e r 2 ( 0 , 1 ) , where S 1 is the projection of S on V 0 , i.e., S 1 = k K d k e 2 π k i x and S 2 = S S 1 . Then, given ( w 0 , v 0 , T 0 , S 0 ) Y , we decompose T 0 = T 0 1 + T 0 2 , S 0 = S 0 1 + S 0 2 , and T ( t ) = T 1 ( t ) + T 2 ( t ) , S ( t ) = S 1 ( t ) + S 2 ( t ) and we consider ( w ( t ) , v ( t ) , T 1 ( t ) , S 1 ( t ) ) M and
( w ( t ) , v ( t ) , T ( t ) , S ( t ) ) ( w ( t ) , v ( t ) , T 1 ( t ) , S 1 ( t ) ) = ( 0 , 0 , T 2 ( t ) , S 2 ( t ) ) .
From (A64), and taking into account that b k = 0 for k I Z * , we derive | a k ( t ) | | a k 0 | e 4 ν π 2 k 2 t ; moreover, with 4 ν π 2 k 2 t 4 ν π 2 t for every k I Z * with (A45), it implies that T 2 ( t ) H ˙ p e r 1 T 0 2 H ˙ p e r 1 e 4 ν π 2 t i.e., T 2 ( t ) 0 in H ˙ p e r 1 ( 0 , 1 ) if t .
Moreover, we have S 2 ( t ) = k I Z * \ K d k ( t ) e 2 π k i x ; therefore,
S 2 ( t ) L ˙ p e r 2 2 = k I Z * \ K | d k ( t ) | 2 .
Since b k = 0 for k I Z * \ K , from (A67), we have
| d k ( t ) | 2 | d k 0 | 2 e 8 c π 2 k 2 t .
Thus,
S 2 ( t ) L ˙ p e r 2 2 4 e 8 c π 2 t S 20 L ˙ p e r 2 2 .
Therefore, T 2 ( t ) H ˙ p e r 1 and S 2 ( t ) L ˙ p e r 2 0 as t with the exponential decay rate e 4 π 2 δ t , where δ = min { ν , c } . Thus, M attracts ( w ( t ) , v ( t ) , T ( t ) , S ( t ) ) with the exponential rate e 4 π 2 δ t in I R 2 × H ˙ p e r 1 ( 0 , 1 ) × L ˙ p e r 1 ( 0 , 1 ) . □
Remark A3. 
If T 0 , S 0 H ˙ p e r m ( 0 , 1 ) , from | a k ( t ) | | a k 0 | e 4 δ π 2 t and taking into account of (A45), we have T 2 ( t ) H ˙ p e r m ( 0 , 1 ) e 4 δ π 2 t T 0 2 H ˙ p e r m ( 0 , 1 ) ; Furthermore, we note that, by working as above, we have:
S 2 ( t ) H ˙ p e r m ( 0 , 1 ) 2 4 e 8 δ π 2 t S 20 H ˙ p e r m ( 0 , 1 ) 2
and the invariant M , attracts the solutions ( w ( t ) , v ( t ) , T ( t ) , S ( t ) ) in
I R 2 × H ˙ p e r m ( 0 , 1 ) × H ˙ p e r m ( 0 , 1 )
with exponential rate e 4 π 2 δ t .

Appendix A.6. The Explicit Reduced Subsystem

Under the hypotheses and notation of Theorem A3, we suppose that
f ( x ) = k J c k e 2 π k i x ,
with c k 0 for every k J I Z and c k = 0 if k J . . On the inertial manifold
( T ( t , x ) S ( t , x ) ) f ( x ) d x = k K a k ( t ) · c k k K d k ( t ) · c k = K J ( a k ( t ) d k ( t ) ) c k .
Hence, the behaviors of velocity v and acceleration w are influenced by the Fourier coefficients of T and S within the set K J . To solve the complete system, we first determine the coefficients of T and S belonging to K J , and then proceed to solve the equations for the coefficients of T and S, where k K J (i.e., | K \ ( K J ) | ).
It is worth noting that 0 K J . Since K = K and J = J , the set K J contains an even number of elements, denoted as 2 n 0 . Therefore, the number of positive elements in K J , denoted as ( K J ) + , is n 0 .
Corollary A1. 
Under the assumptions and notation of Theorem A3, if we assume that the set K J is finite with | K J | = 2 n 0 , then the asymptotic behavior of system (A7) can be described by a system of N = 4 n 0 + 2 coupled equations in I Z R N . These equations determine the variables ( w , v , a k , d k ) for k K J , along with a family of | K \ ( K J ) | linear non-autonomous equations.
Proof. 
On the inertial manifold:
( T ( t , x ) S ( t , x ) ) f ( x ) d x = k K ( a k d k ) ( t ) c k = k K J a k ( t ) c k k K J d k ( t ) c k .
Consequently, the system’s dynamics rely on the coefficients within K J . Additionally, the equations for a k and d k are conjugates of the equations for a k and d k , respectively. This implies that k K J a k ( t ) c k = 2 R e k ( K J ) + a k ( t ) c k and k K J d k ( t ) c k = 2 R e k ( K J ) + d k ( t ) c k .
From this, taking real and imaginary parts of a k , ( a 1 k , a 2 k ) and d k , ( d 1 k , d 2 k ) , k ( K J ) + in (A53) with n 0 = | ( K J ) + | , we conclude. □
Remark A4. 
Taking the real and imaginary parts of the coefficients of the temperature, a k ( t ) , the heat flux at the wall of the loop, b k , the geometry of the circuit, c k , and the solute concentration, d k ( t ) , k ( K J ) + , as
a k ( t ) = a 1 k ( t ) + i a 2 k ( t ) , b k = b 1 k + i b 2 k , c k = c 1 k + i c 2 k , d k ( t ) = d 1 k ( t ) + i d 2 k ( t )
the asymptotic behavior of system (A7) is given by a reduced explicit system of ODEs in I R N , with N = 4 n 0 + 2 , given by
d w d t + 1 ε w + 1 ε G ( v ) v ( t ) = 1 ε 2 k ( k J ) + [ a 2 k ( t ) c 2 k a 1 k ( t ) c 1 k ] 1 ε 2 k ( k J ) + [ d 2 k ( t ) c 2 k d 1 k ( t ) c 1 k ] d v d t = w d ( a 1 k ) d t + [ 4 π 2 k 2 ν a 1 k ( t ) 2 π k v ( t ) a 2 k ( t ) ] = b 1 k , k ( K J ) + d ( a 2 k ) d t + [ 2 π k v ( t ) a 1 k ( t ) + 4 π 2 k 2 ν a 2 k ( t ) ] = b 2 k , k ( K J ) + d ( d 1 k ) d t + [ 4 c π 2 k 2 d 1 k ( t ) 2 π k v ( t ) d 2 k ( t ) ] = 4 b π 2 k 2 a 1 k ( t ) , k ( K J ) + d ( d 2 k ) d t + [ 4 c π 2 k 2 d 2 k ( t ) + 2 π k v ( t ) d 1 k ( t ) ] = 4 b π 2 k 2 a 2 k ( t ) , k ( K J ) +
Thus, we reduced the asymptotic behavior of the initial system (A7) to the dynamics of the reduced explicit system (A72). It is worth noting that, from the above analysis, it is possible to design the geometry of the circuit and/or the external heating source, by properly choosing the functions f and/or the ambient temperature, h, so that the resulting system has an arbitrary number of equations of the form N = 4 n + 2 .
Note that K and J may be infinite sets, but their intersection is finite. For instance, for a circular circuit, we have f ( x ) a sin ( x ) + b cos ( x ) , i.e., J = { ± 1 } , and then K J is either { ± 1 } or the empty set.
If we take, for the sake of simplicity, a circular geometry, then J = { ± 1 } and K J = { ± 1 } . Also, if we take k = 1 and omit the equation for k , the conjugate of k, we have the following transformed set of equations:
d w d t = 2 R e ( a 1 ( t ) c 1 ) ε 2 R e ( d 1 ( t ) c 1 ) ε G ( v ) v ε w ε , d v d t = w , d ( a 1 ) d t = l ( v ) b 1 l ( v ) a 1 ( t ) 4 ν π 2 a 1 ( t ) 2 π v a 1 ( t ) i , d ( d 1 ) d t = 4 c π 2 d 1 ( t ) + 4 b π 2 a 1 ( t ) 2 π v d 1 ( t ) i
where the variables of interest are w ( t ) , representing the fluid acceleration, v ( t ) , the fluid velocity, a 1 ( t ) , the Fourier mode of temperature, and the Fourier mode of solute concentration, d 1 ( t ) . To simplify the notation, we can omit the subscripts as we only have a single coefficient.
To reduce the number of independent parameters, we introduce a change of variables: a 1 c 1 a and d 1 c 1 d . Additionally, we express the equations in terms of their real and imaginary components using the following approach:
a ( t ) = a 1 ( t ) + i a 2 ( t ) ,
d ( t ) = d 1 ( t ) + i d 2 ( t )
b 1 = A + i B
with a 1 ( t ) , a 2 ( t ) , d 1 ( t ) , d 2 ( t ) , A, B I R . Finally, we integrate the system to arrive at Equation (19) in Section 3.

References

  1. Lorenz, E.N. Deterministic nonperiodic flow. J. Atmos. Sci. 1963, 20, 130–141. [Google Scholar] [CrossRef]
  2. Strogatz, S.H. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering; CRC Press: Boca Raton, FL, USA, 2018. [Google Scholar]
  3. Ott, E. Chaos in Dynamical Systems; Cambridge University Press: Cambridge, UK, 2002. [Google Scholar]
  4. Shukla, J. Predictability in the midst of chaos: A scientific basis for climate forecasting. Science 1998, 282, 728–731. [Google Scholar] [CrossRef] [PubMed]
  5. Boffetta, G.; Cencini, M.; Falcioni, M.; Vulpiani, A. Predictability: A way to characterize complexity. Phys. Rep. 2002, 356, 367–474. [Google Scholar] [CrossRef] [Green Version]
  6. Jiménez-Casas, A.; Rodríguez-Bernal, A. Finite-Dimensional Asymptotic Behavior in a Thermosyphon Including the Soret Effect. Math. Methods Appl. Sci. 1999, 22, 117–137. [Google Scholar] [CrossRef]
  7. Loiseau, J.C. Data-driven modeling of the chaotic thermal convection in an annular thermosyphon. Theor. Comput. Fluid Dyn. 2020, 34, 339–365. [Google Scholar] [CrossRef]
  8. Malcai, O.; Biham, O.; Richmond, P.; Solomon, S. Theoretical analysis and simulations of the generalized Lotka-Volterra model. Phys. Rev. E 2002, 66, 031102. [Google Scholar] [CrossRef] [Green Version]
  9. García-Algarra, J.; Galeano, J.; Pastor, J.M.; Iriondo, J.M.; Ramasco, J.J. Rethinking the logistic approach for population dynamics of mutualistic interactions. J. Theor. Biol. 2014, 363, 332–343. [Google Scholar] [CrossRef] [Green Version]
  10. Currie, J.; Castro, M.; Lythe, G.; Palmer, E.; Molina-París, C. A stochastic T cell response criterion. J. R. Soc. Interface 2012, 9, 2856–2870. [Google Scholar] [CrossRef] [Green Version]
  11. Castro, M.; Lythe, G.; Molina-Paris, C.; Ribeiro, R.M. Mathematics in modern immunology. Interface Focus 2016, 6, 20150093. [Google Scholar] [CrossRef] [Green Version]
  12. Bravo-Gutiérrez, M.E.; Castro, M.; Hernández-Machado, A.; Poiré, E.C. Controlling Viscoelastic Flow in Microchannels with Slip. Langmuir 2011, 27, 2075–2079. [Google Scholar] [CrossRef]
  13. Jiménez Casas, Á.; Castro Ponce, M. A Thermosyphon Model with a Viscoelastic Binary Fluid. Electron. J. Differ. Equ. 2015, 22, 54–61. [Google Scholar]
  14. Rodríguez-Bernal, A.; Van Vleck, E.S. Diffusion Induced Chaos in a Closed Loop Thermosyphon. SIAM J. Appl. Math. 1998, 58, 1072–1093. [Google Scholar] [CrossRef] [Green Version]
  15. Jiménez-Casas, Á.; Ovejero, A.M.L. Numerical Analysis of a Closed-Loop Thermosyphon Including the Soret Effect. Appl. Math. Comput. 2001, 124, 289–318. [Google Scholar] [CrossRef]
  16. Storey, I.; Bourmistrova, A.; Subic, A. Smooth control over jerk with displacement constraint. Proc. Inst. Mech. Eng. Part J. Mech. Eng. Sci. 2012, 226, 2656–2673. [Google Scholar] [CrossRef]
  17. Lu, Y.S.; Lin, Y.Y. Smooth motion control of rigid robotic manipulators with constraints on high-order kinematic variables. Mechatronics 2018, 49, 11–25. [Google Scholar] [CrossRef]
  18. Park, J.H.; Waickman, A.T.; Reynolds, J.; Castro, M.; Molina-París, C. IL7 receptor signaling in T cells: A mathematical modeling perspective. Wiley Interdiscip. Rev. Syst. Biol. Med. 2019, 11, e1447. [Google Scholar] [CrossRef]
  19. Arias, C.F.; Herrero, M.A.; Acosta, F.J.; Fernandez-Arias, C. Population mechanics: A mathematical framework to study T cell homeostasis. Sci. Rep. 2017, 7, 9511. [Google Scholar] [CrossRef]
  20. Castro, M.; van Santen, H.M.; Férez, M.; Alarcón, B.; Lythe, G.; Molina-París, C. Receptor pre-clustering and T cell responses: Insights into molecular mechanisms. Front. Immunol. 2014, 5, 132. [Google Scholar] [CrossRef] [Green Version]
  21. Samardzija, N.; Greller, L.D. Explosive route to chaos through a fractal torus in a generalized Lotka-Volterra model. Bull. Math. Biol. 1988, 50, 465–491. [Google Scholar] [CrossRef]
  22. Morrison, F.A. Understanding Rheology; Oxford University Press: New York, NY, USA, 2001; Volume 1. [Google Scholar]
  23. Thurston, G.B.; Henderson, N.M. Effects of Flow Geometry on Blood Viscoelasticity. Biorheology 2006, 43, 729–746. [Google Scholar]
  24. Welander, P. On the Oscillatory Instability of a Differentially Heated Fluid Loop. J. Fluid Mech. 1967, 29, 17–30. [Google Scholar] [CrossRef]
  25. Liñan, A. Analytical Description of Chaotic Oscillations in a Toroidal Thermosyphon. In Fluid Physics, Lecture Notes of Summer Schools; World Scientific: London, UK, 1994; p. 17. [Google Scholar]
  26. Herrero, M.A.; Velazquez, J.J.L. Stability Analysis of a Closed Thermosyphon. Eur. J. Appl. Math. 1990, 1, 1–23. [Google Scholar] [CrossRef]
  27. Keller, J.B. Periodic Oscillations in a Model of Thermal Convection. J. Fluid Mech. 1966, 26, 599–606. [Google Scholar] [CrossRef]
  28. Yasappan, J.; Jiménez-Casas, Á.; Castro, M. Stabilizing Interplay between Thermodiffusion and Viscoelasticity in a Closed-Loop Thermosyphon. Discret. Contin. Dyn. Syst. B 2015, 20, 3267–3299. [Google Scholar] [CrossRef] [Green Version]
  29. Jiménez-Casas, Á. Asymptotic Behaviour of a Closed-Loop Thermosyphon with Linear Friction and Viscoelastic Binary Fluid. Math. Methods Appl. Sci. 2019, 42, 6791–6809. [Google Scholar] [CrossRef]
  30. Hart, J.E. A Model of Flow in a Closed-Loop Thermosyphon Including the Soret Effect. J. Heat Transf. 1985, 107, 840–849. [Google Scholar] [CrossRef]
  31. Velázquez, J.J.L. On the Dynamics of a Closed Thermosyphon. SIAM J. Appl. Math. 1994, 54, 1561–1593. [Google Scholar] [CrossRef]
  32. Rodríguez-Bernal, A.; Van Vleck, E.S. Complex Oscillations in a Closed Thermosyphon. Int. J. Bifurc. Chaos 1998, 8, 41–56. [Google Scholar] [CrossRef]
  33. Henry, D. Geometric Theory of Semilinear Parabolic Equations; Springer: Berlin/Heidelberg, Germany, 2006; Volume 840. [Google Scholar]
  34. Hale, J.K. Asymptotic Behavior of Dissipative Systems; Number 25; American Mathematical Society: Providence, RI, USA, 2010. [Google Scholar]
  35. Bloch, A.M.; Titi, E.S. On the dynamics of rotating elastic beams. In New Trends in Systems Theory; Springer: Berlin/Heidelberg, Germany, 1991; pp. 128–135. [Google Scholar]
  36. Stuart, A. Perturbation theory for infinite dimensional dynamical systems. In Advances in Numerical Analysis; Oxford University Press: Oxford, UK, 1995. [Google Scholar]
  37. Rodriguez-Bernal, A. Attractors and Inertial Manifolds for the Dynamics of a Closed Thermosiphon. J. Math. Anal. Appl. 1995, 193, 942–965. [Google Scholar] [CrossRef] [Green Version]
  38. Foias, C.; Sell, G.R.; Temam, R. Inertial Manifolds for Nonlinear Evolutionary Equations. J. Differ. Equ. 1988, 73, 309–353. [Google Scholar] [CrossRef] [Green Version]
  39. Rodriguez Bernal, A. Inertial Manifolds for Dissipative Semiflows in Banach Spaces. Appl. Anal. 1990, 37, 95–141. [Google Scholar] [CrossRef]
Figure 1. A cartoon of a linear-jerk-controlled chaotic waterwheel. Water inflow fills the water cups, creating a torque that keeps the wheel rotating. At the center of the waterwheel, a controlled spring is designed to control the jerk and avoid (in principle) large pulls.
Figure 1. A cartoon of a linear-jerk-controlled chaotic waterwheel. Water inflow fills the water cups, creating a torque that keeps the wheel rotating. At the center of the waterwheel, a controlled spring is designed to control the jerk and avoid (in principle) large pulls.
Mathematics 11 03099 g001
Figure 2. Bifurcation diagram Model M1. We only vary the viscoelastic parameter ε to emphasize the role of elasticity in the chaotic/non-chaotic nature of the solutions. Despite the fact that elasticity is associated with oscillatory motion, we show that new patterns of chaotic behavior arise for large values of ε .
Figure 2. Bifurcation diagram Model M1. We only vary the viscoelastic parameter ε to emphasize the role of elasticity in the chaotic/non-chaotic nature of the solutions. Despite the fact that elasticity is associated with oscillatory motion, we show that new patterns of chaotic behavior arise for large values of ε .
Mathematics 11 03099 g002
Figure 3. Temperature velocity phase diagram ( a 1 , a 2 , v ) for Model M1.
Figure 3. Temperature velocity phase diagram ( a 1 , a 2 , v ) for Model M1.
Mathematics 11 03099 g003
Figure 4. Example of velocity evolution for ε = 10 2 for Model M1.
Figure 4. Example of velocity evolution for ε = 10 2 for Model M1.
Mathematics 11 03099 g004
Figure 5. Temperature velocity phase diagram ( a 1 , a 2 , v ) for Model M1.
Figure 5. Temperature velocity phase diagram ( a 1 , a 2 , v ) for Model M1.
Mathematics 11 03099 g005
Figure 6. Comparison of four velocity time-courses for different values of ε and fixing the rest of the parameters of Model M1. Note how small changes in ε between panels (c,d) also induce transient chaotic behavior followed by stability at a fixed point.
Figure 6. Comparison of four velocity time-courses for different values of ε and fixing the rest of the parameters of Model M1. Note how small changes in ε between panels (c,d) also induce transient chaotic behavior followed by stability at a fixed point.
Mathematics 11 03099 g006
Figure 7. Velocity evolution on time.
Figure 7. Velocity evolution on time.
Mathematics 11 03099 g007
Figure 8. Velocity evolution on time.
Figure 8. Velocity evolution on time.
Mathematics 11 03099 g008
Figure 9. Bifurcation diagram, Model M2, compared with that in Figure 2. Again, elasticity does not mean periodic behavior; rather, new complex regions of chaotic-stable transitions emerge.
Figure 9. Bifurcation diagram, Model M2, compared with that in Figure 2. Again, elasticity does not mean periodic behavior; rather, new complex regions of chaotic-stable transitions emerge.
Mathematics 11 03099 g009
Figure 10. Phase diagram ( x , y , z ) for Model M2.
Figure 10. Phase diagram ( x , y , z ) for Model M2.
Mathematics 11 03099 g010
Figure 11. Phase diagram ( x , y , z ) for Model M2.
Figure 11. Phase diagram ( x , y , z ) for Model M2.
Mathematics 11 03099 g011
Figure 12. Phase diagram ( x , y , z ) for Model M2.
Figure 12. Phase diagram ( x , y , z ) for Model M2.
Mathematics 11 03099 g012
Figure 13. Comparison of bifurcation diagrams for the Lorenz parameter ρ with and without elasticity. The Lorenz system, (22), is represented with blue dots, and the system (21) is represented with green dots.
Figure 13. Comparison of bifurcation diagrams for the Lorenz parameter ρ with and without elasticity. The Lorenz system, (22), is represented with blue dots, and the system (21) is represented with green dots.
Mathematics 11 03099 g013
Figure 14. Bifurcation diagram for Model M3. As discussed for models M1 and M2, the role of elasticity is not to force oscillations (or resonance) but to change the chaotic behavior dramatically.
Figure 14. Bifurcation diagram for Model M3. As discussed for models M1 and M2, the role of elasticity is not to force oscillations (or resonance) but to change the chaotic behavior dramatically.
Mathematics 11 03099 g014
Figure 15. Velocity evolution for Model M3.
Figure 15. Velocity evolution for Model M3.
Mathematics 11 03099 g015
Figure 16. Velocity evolution for ε = 10 1 for Model M3.
Figure 16. Velocity evolution for ε = 10 1 for Model M3.
Mathematics 11 03099 g016
Figure 17. Temperature phase diagram for ε = 14 .
Figure 17. Temperature phase diagram for ε = 14 .
Mathematics 11 03099 g017
Figure 18. Maximum Lyapunov Exponent. This diagram summarizes the main idea of the work: elasticity is a non-trivial contributor to the emergence or disappearance of chaotic behavior.
Figure 18. Maximum Lyapunov Exponent. This diagram summarizes the main idea of the work: elasticity is a non-trivial contributor to the emergence or disappearance of chaotic behavior.
Mathematics 11 03099 g018
Table 1. Categorization of Models Based on the viscoelastic parameter ε . For the sake of simplicity, we group stable and oscillatory behaviors under the label non-chaotic.
Table 1. Categorization of Models Based on the viscoelastic parameter ε . For the sake of simplicity, we group stable and oscillatory behaviors under the label non-chaotic.
ModelChaosNon-ChaoticChaosNon-Chaotic
M1 10 7 < ε < 10 3 10 3 < ε < 10 10 < ε < 12 12 < ε
M2 10 7 < ε < 10 2 10 1 < ε < 1 1.2 < ε < 1.5 ε > 1.5
M3 10 7 < ε < 10 2 10 2 < ε < 12 14 < ε < 15 ε > 15
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Jiménez-Casas, Á.; Castro, M.; Villanueva-Pesqueira, M. The Role of Elasticity on Chaotic Dynamics: Insights from Mechanics, Immunology, Ecology, and Rheology. Mathematics 2023, 11, 3099. https://doi.org/10.3390/math11143099

AMA Style

Jiménez-Casas Á, Castro M, Villanueva-Pesqueira M. The Role of Elasticity on Chaotic Dynamics: Insights from Mechanics, Immunology, Ecology, and Rheology. Mathematics. 2023; 11(14):3099. https://doi.org/10.3390/math11143099

Chicago/Turabian Style

Jiménez-Casas, Ángela, Mario Castro, and Manuel Villanueva-Pesqueira. 2023. "The Role of Elasticity on Chaotic Dynamics: Insights from Mechanics, Immunology, Ecology, and Rheology" Mathematics 11, no. 14: 3099. https://doi.org/10.3390/math11143099

APA Style

Jiménez-Casas, Á., Castro, M., & Villanueva-Pesqueira, M. (2023). The Role of Elasticity on Chaotic Dynamics: Insights from Mechanics, Immunology, Ecology, and Rheology. Mathematics, 11(14), 3099. https://doi.org/10.3390/math11143099

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop