Marc Heddebaut

Souheir Mili

David Sodoyer

Eduardo Jacob

Marina Aguado

Christian Pinedo Zamalloa

Igor Lopez

Virginie Deniau
email: virginie.deniau@ifsttar.fr.

Towards

Towards a resilient railway communication network against electromagnetic attacks

Keywords: Railway, radio, attack, jammer, architecture, resilience

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Résumé

En raison de son déploiement très vaste et de sa relative facilité d'accès, l'infrastructure ferroviaire apparait susceptible de subir d'éventuelles attaques. Celles-ci pourraient avoir des conséquences significatives sur l'économie ainsi qu'en termes de sécurité. Ces dernières années, des efforts considérables ont été déployés en vue de développer un système de gestion du trafic ferroviaire paneuropéen interopérable. Ce systèmes sont de ce fait bien documentés et leurs caractéristiques sont également largement diffusées et disponibles. Des communications radio de différentes natures sont largement employées pour gérer ce trafic ferroviaire. Dans ce contexte, le projet européen «SECRET» vise à développer des solutions innovantes afin d'améliorer la résilience électromagnétique des systèmes radio ferroviaires. Cet article explore différents aspects de recherche du projet. Une première partie analyse l'infrastructure ferroviaire et les risques d'attaques électromagnétiques. La section suivante explore certaines des conditions qui peuvent inhiber la radio sol-train. Afin de contrer ces systèmes de brouillage, une architecture de communication résiliente est analysée en dernière section.

Mots-clé: Transport ferroviaire ; radio ; attaque ; brouilleur ; architecture ; résilience

Introduction

The terrorism threat to European citizens has been constantly elevated for several years. This threat also concerns rail transportation and its associated infrastructure which provide mass transport. Both rail operation and its infrastructure are recognized as critical priorities because of the economic and security impacts of terrorist attacks (loss of service, destruction of vehicles, and destruction of infrastructure…). Extended consequences on the surrounding businesses are also expected, as well as the impaired reputation of the railway as a safe and secure transport system. Indeed, the railway is an attractive target for security attacks, because of its familiarity, ease of access and openness. The attacks carried out on the French railway infrastructure in November 2008 involving metal bars on the catenaries which caused delays for 160 TGVs, Eurostar and Thalys (coordinated attacks that targeted four different rail lines on the French Railway network) showed that more subtle actions than bombings can be carried out by terrorists: the terrorists can base their action methods on the vulnerabilities of technologies employed. Equally, electromagnetic (EM) terrorism is based on a similar approach, as it consists in failing equipment or devices which serve the efficiency and safety of railway transport systems.

The European railway network constitutes a mass transport system which includes a large number of telecommunication, command-control, electronic and computer systems and subsystems potentially vulnerable to intentional electromagnetic interferences (IEMI). Generally, electromagnetic interferences (EMI) can act on wireless communications, on wired networks and on electronic and computer components and systems. They can have three main effects:

-Destruction of the electronic components; jamming information transmitted and, thus, impeding communication between two components; modifying the information transmitted and, thus, enabling or disabling certain functions unexpectedly.

Considering the current evolution of the European rail network which is based on the deployment of the GSM-Railway communication system to transmit the signalling information between the control-centres and the trains, the first feared effect is the jamming of communication signals. Indeed, the jamming of wireless communication signals can be easily implemented in comparison to other techniques aiming to modifying the carried information. However, today, the real consequences of such attacks on the whole railway network are not wellknown.

Consequently, the European project SECRET (SECurity of Railways against Electromagnetic aTtacks) aims to assess the real risks concerning EM attacks on rail networks, to identify areas for strengthening and to develop a detection and management system for EM attacks that is integrated into the rail infrastructure, making it an architecture resilient to any EM attack. Obviously, to be in agreement with the current evolution of the European railway network, the architecture has to be adapted to interoperability needs associated with the current harmonization process of the European railway network.

Knowing that, within the railway, information exchanges between trains and the infrastructure are necessarily provided by wireless links, it is, therefore, impossible to protect such information from EM attacks with conventional shields, filters or other topological solutions. Thus, it becomes necessary to provide a dynamic solution based on the detection of the EM attack and on the redundancy of communication solutions to ensure that the correct information reaches the recipient (infrastructure or train) in an adequate amount of time.

This research article is structured as follows. Section 2 underpins the potential risk that EM attacks may introduce in railway operation: our research motivation. Section 3 determines the different railway radio communication subject to electromagnetic attacks and focuses our research on a specific target scenario: the GSM-R communication link. Next, we identify current jamming systems and measure their effectiveness to disturb GSM-R communication. Based on real measurements, we also provide a necessary baseline of EM behaviour characterization for normal GSM-R radio operation in the railway scenario. Once, we identify the target communication to be protected, the attack scenario and the effect of the EM attack in the GSM-R communication link, section 4 describes the resilient railway communication network able to react in real time to this attack scenario. The proposed architecture is grounded on two major concepts: a Detection System, which will use of the detection techniques of GSM-R attacks detailed in the section 3 among others, and a Multipath Communication System, which will strengthen the connections between trains and trackside.

EM attacks and railway network

EM attacks consist of intentionally generated EM signals in order to induce EM disturbances on electronic, computer or communication components. A definition of Intentional Electromagnetic Interferences (IEMI) is given by the International Electrotechnical Commission (IEC) as "intentional malicious generation of electromagnetic energy introducing noise or signals into electric and electronic systems thus, disrupting, confusing, or damaging these systems for terrorist or criminal purposes."

In general, we can assume that EM attacks involving sufficiently high-level power interference can produce permanent damages to electronic equipment, and that lower power level EM interferences can mainly impact the quality of wireless transmissions.

In the case of EM attacks on railways, a large number of components of railway infrastructure can be affected: control-command systems (track switches, barrier crossing, traffic signal light...), spot communication systems for detection (track circuits, eurobalise, euroloop), continuous communication systems (local train radio systems, GSM-Railway). Consequently, the railway infrastructure is an effective potential victim of electromagnetism terrorism. According to the systems or subsystems attacked and to the nature of the intentional interferences, the consequences on the individuals, the structures and the materials can significantly vary.

Moreover, the deployment of ERTMS (European Rail Traffic Management System) not only homogenizes the technologies to manage the trains over the European territory, but also the vulnerability points to EM interferences. Some examples can be set up in the context of the harmonization of systems and rules due to interoperability requirements for all operational domains in the rail transportation in Europe (operation, control, management, maintenance…). These include the strategy to reduce the number of control centres of the track switches in Europe that will rely on interlocking remote controlled systems in order to activate switches. Thus, if a terrorist is designing an intentional EM emissions device capable of disrupting management systems for rail infrastructure in Berlin, for instance, the same device will have the same attack capacity in all European cities. This will cause immediate economic consequences and possibly more…. Harmonization thus facilitates the implementation of organized and simultaneous EM attacks.

The easiness of implementation of an EM attack will notably depend on the accessibility of the required devices to generate a given EM attack signal. Today, with the proliferation of the telecommunication applications, it is really easy to obtain equipment providing relatively low-level power interference able to jam wireless communication signals. Consequently, in SECRET, the priority is given to the study of the potential impacts of low power interference sources, such as jammers, on wireless railway communication systems.

Meanwhile, the technologies and frequencies employed in the railway field can be similar to technologies and frequencies used for applications available to the general public. Indeed, the railway no longer develops technology "owners" but adapts general public technologies. This increases the vulnerability of the railway because it is easy to obtain emission devices capable of disrupting rail technologies. With a relatively low basic knowledge of electronics and the performance of electronic components and antennas available on the open market, these emission devices can be combined with amplifiers to increase the capacity of EM attacks.

Given the potential vulnerability of the railway network and the ease of implementation of such EM attacks, we chose to work on an integrated solution to ERTMS offering the capability of detecting EM attack situations and engaging appropriate responses to maintain the security and capability of the railway network.

Detection of EM attack signals

Critical railway radio targets

A preliminary step of this study was to determine the potential radio interface targets for electromagnetic threats and malicious attacks. They are mainly related to the use of the GSM-Railway (GSM-R) for train to ground communication, the use of GNSS/GPS satellite navigation received signals for the train odometry system and the Eurobalise/Balise Transmission Module for spot communication. TETRA, Terrestrial Trunked Radio, which was also considered many years ago for ground to train radio communication, is also a potential target. TETRA is currently used by some railway operators for railway and depot communication and also for railway security applications. This paper will concentrate on GSM-R.

GSM-R presentation

GSM-R is a standard communication platform for railways. It is a strategic communication system focused on the interoperability between European railway infrastructures. By the end of 2016, 56 countries in five continents should have operational GSM-R networks. Its specifications are widely disseminated. Europe is leading GSM-R implementation in 11 countries (Austria, Belgium, Czech Republic, France, Germany, Italy, the Netherlands, Norway, Spain, Sweden and Switzerland). In 2011, GSM-R deployment provided coverage of about 30% of the European railway network with 68,000 km in operation; 156,000 km are planned to be covered (70% of the European railway network). The main objectives of GSM-R are to replace analogical radio communication and to provide a unique data transmission solution for ERTMS/ETCS (European Rail Traffic Management System / European Train Control System) level 2 and level 3. GSM-R is an evolution of public GSM dedicated to railway application. Therefore, GSM-R has similar characteristics to those of public GSM system. GSM-R functionality is built on standards and recommendations supported by mainly two organizations, ETSI (European Telecommunications Standards Institute) and UIC (Union Internationale des Chemins de fer, International Union of Railways). ETSI technical committee RT (Railway Telecommunications) is responsible for those aspects of Global System for Mobile communications standardization which are specific to Railway (GSM-R) and Private Mobile Radio (PMR) operations. UIC specifies through its EIRENE (European Integrated Railway Radio Enhanced Network) project, a digital radio standard for the European railways. It forms part of the specification for technical interoperability.

GSM-R performances are guaranteed for high-speed conditions (up to 500 km/h). A common European frequency band is allocated to GSM-R below the frequencies of the public GSM standard. The allocated frequency bands are for the uplink: 876 MHz -880 MHz and for the downlink: 921 MHz -925 MHz. As represented in figure 1, GSM-R has its own cellular network installed along the track designed to provide a minimum level of received power > -95 dBm anywhere along the track. This is a mandatory requirement coming from the EIRENE specifications. This represents a limited level of power with the potential of being jammed. -Type 'D' devices 'Direct Receive and Transmit Jammers': These devices operate as a small independent base station. The jammer is predominantly in receiving mode and will choose intelligently interaction and blocking the cell phone if it is within close proximity of the jammer. -Type "E" devices 'EMI Shield -Passive Jamming': These jamming solutions consist in suppressing electromagnetic signals by using the properties of the Faraday cages. The Faraday cage essentially blocks, or greatly attenuates, all electromagnetic signals entering or leaving the cage.

In this presentation, we shall only consider type 'A' jammers. The next step of this study was to determine the effectiveness of these particular jammers.

Jamming effectiveness

To evaluate the effectiveness of any particular type 'A' jammer, a dedicated test bench was set up. It is presented in figure 1. A GSM-R communication is established in laboratory using a base transceiver station (BTS) emulator (CMU) and GSM-R train mobile equipment. Using combiners, a supplementary radio signal is injected in the radiofrequency link generating different jamming waveforms, at different power levels. A spectrum analyzer is used to visualize these waveform frequency occupancies [START_REF] Mili | Recognition process of jamming signals superimposed on GSM-R radiocommunications[END_REF]. In the first column, a constant -38 dBm GSM-R level is delivered to the input of the mobile GSM-R receiver.

CMU: GSM-R base

This represents a comfortable level of operation, as compared to the -95 dBm minimum level requested by the EIRENE specifications. The second column indicates the variable continuous wave jamming power injected at the input of the GSM-R mobile receiver. This jamming signal is centred in the GSM-R communication channel. The third column shows the resulting measured bit error rate (BER) in the case of a pure sinusoidal jammer. The last column indicates the results obtained with FM modulated jamming signals. We obtain that:

-When the received GSM-R power is of the same order of magnitude as the power received from the jammer then, the communication is perturbed or interrupted.

-A 6 dB difference (or more) between the two power levels greatly reduces the impact of the jammer.

-A pure sine wave centered on the GSM-R channel has less impact than a wider bandwidth FM modulated jamming signal.

We conclude that limited power jammers can effectively affect GSM-R communication.

The next step of the study was to evaluate "normal" EM environments in order to be able to then detect "jammed" electromagnetic conditions. Therefore, it was decided to perform electromagnetic environment measurements in different critical railway environments i.e. in stations, along railway tracks and in trains.

Figure 3 represents one typical result of spectral measurements performed in a railway station (on the left side, a 3D representation: frequency, time, and received power and, on the right side, the equivalent waterfall representation). The results presented cover a bandwidth of 10 MHz. As stated before, the GSM-R band occupies the 4 MHz central part of this 10 MHz band. A single GSM-R beacon channel, providing strong radio coverage of the station, is received in the GSM-R band; Line of Sight (LOS) of the corresponding GSM-R BTS antenna is nearly possible from this measurement location. Another, weaker, GSM-R BTS signal is also discernible higher in frequency, corresponding to a farther GSM-R BTS. The upper part of the band is fairly busy due to the number of cellular phone users in the station and its surroundings. The lower part of the monitored band is free of activity at the particular time of this recording.

After determining "normal" electromagnetic ambient conditions, sensors shall be installed in order to detect potential attacks. They will use possible discrepancies between the observed electromagnetic environments and the recorded "normal" electromagnetic environments. Sensor outputs deliver data to an optimized communication architecture providing the capability to react against these attacks.

A resilient communication architecture against EM attacks

This section describes the proposed resilient railway communication architecture able to detect and react against an EM or jamming attack similar to those described in section 3 and, consequently, able to manage the disruption of the train-trackside communication.

As previously introduced in section 3, ERTMS Levels 2 and 3 require a permanent communication path between trackside and the train so that the Radio Block Centre (RBC) can signal the train. All the possible communication paths between the trackside and the train rely on wireless communication technologies, at least in a certain section of the communication path. Up to four wireless technologies are defined in the current specifications of ERTMS: GSM-R, Eurobalise, Euroloop and Radio-infill. The usage and importance of the previous technologies depends on the ERTMS Level deployed in the track.

Facing the impossibility to absolutely guarantee wireless communication between trackside and the train during EM attacks, due to the open-air interface inherent features, the aim of the proposed resilient architecture is to provide a flexible architecture that will provide multiple ways to make more resilient trackside-train communications depending on, on the one hand, on the ERTMS deployment and, on the other hand, on the security level the management staff wishes.

The architecture is based on two main concepts: detection of EM attacks and a more resilient wireless communication path. These two main concepts are the basis of the resilient architecture that consists of two subsystems: the Detection System (DS) and the Multipath Communication System (MCS) presented in figure 5. Firstly, based on the data sent by the sensors, the DS allows the detection in real time of the EM attacks that occur in the railway infrastructure. Thus, the management staff can be provided with accurate information which allows them to anticipate a loss scenario of connectivity in certain sections of the trackside exposed to an EM attack or, in the worst case scenario, resolve quickly if the loss of communication with a train is due to an EM attack or due to a different reason (lack of power, hardware failure, …). The availability of this information for the railway management staff is crucial to perform the suitable actions in order to restore communications and normal operation. Given the information, management staff could deploy manual reactive actions or even, depending on the railway management system, automatic reactive actions. Secondly, the MCS aims to provide multiple paths of communication between the RBC and the train in order to maintain the ETCS session of ERTMS Level 2 or Level 3, even if one or several paths fail. This implies the requirement of deploying an end-to-end multipath solution between the RBC and the on-board ERTMS located in the train.

The resilience level provided by MCS will depend on the dissimilarity level of the established paths, especially the wireless sections, which are the most sensitive to EM attack. Thus, the use of different wireless access technologies or, even better, frequencies would be desirable.

Finally, deploying the DS and the MCS together is of great interest because both subsystems are complementary. On the one hand, the DS might provide useful information to the MCS in order to change the active interface in case an attack that affects the current active interface is being performed. Thus, the MCS might change the communication paths based on the information provided by the DS. On the other hand, the DS could profit from the more resilient communication provided by the MCS to increase the resiliency of the communication between the trackside detection subsystem and the on-board detection subsystem located in the train.

The Detection System (DS) architecture

The main aim of the Detection System (DS) is to detect EM attacks that are performed in the protected infrastructure and, thus, be able to react to overcome them. The DS presents a geographically distributed architecture composed of one Central Health/Attack Manager (CHAM) that would be generally located at the command centre, and multiple detection subsystems distributed along the track and inside trains. Figure 6 presents this architecture.

The CHAM is the main element of the DS. It consolidates the information provided by its detection subsystem with regard to the information EM attacks and, thus, it provides a complete view of the EM state of the railway infrastructure. The CHAM also eases the centralized management of the detection subsystems. The detection subsystems are autonomous systems, governed by a local Health/Attack Manager (HAM), which is responsible for the detection of EM attacks in a specific geographical area and for carrying out the pertinent actions when EM attacks are detected. These subsystems may be deployed along the track or in the trains.

Inside each detection subsystem, sensors are deployed to feed with information regarding EM attacks the HAM that governs the local detection subsystem. Different types of sensors may be deployed; and the information provided by the sensors and the processing requirement of that information varies depending on the type of sensor used. For example, it can be the radio frequency signal received by an electromagnetic sensor or even information from upper layers like the latency, bit error rate or number of retransmissions.

In order to deploy sensors for detecting EM attacks inside GSM-R frequency bands with the processing algorithms proposed in section 3, two categories of electromagnetic sensor can be used. A first category considers sensors or inputs that provide access to information inside the GSM-R receiver situated in the train or at ground level, at different stages along the receiving chain. This allows us to collect information before/after different filters in the GSM-R receiver. This approach requires the modification of the existing equipment in order to collect the necessary analogue data. A second category considers a sensor independent of the GSM-R equipment, receiving the same radiofrequency signals as those effectively received by the railway equipment, by sharing the electromagnetic signals received by the train or ground BTS antenna.

Regarding the actions performed when an EM attack is discovered, they may depend on the current configuration of the detection subsystem:  Discard the EM attack.  Report the information about the EM attack to the CHAM in order to inform the railway management staff.  In case of availability of a Multipath Communication System (MCS), change the communications configuration of the MCS to adapt communications to the on-going EM attack.

Multipath Communication System (MCS) architecture

The aim of the Multipath Communication System (MCS) is to offer redundant communication paths between the Radio Block Centre (RBC) and the train. Figure 7 presents the architecture of this multipath communication system. The current ERTMS specification is based on circuit-switching technology which involves the use of a reserved channel. This supposes an inefficient use of the limited spectrum resources of GSM-R that would be increased in the case of performing a multipath technology on it. In addition, the deployment of a multipath communication system involves a modification of communication protocols or the application layer. In the case of an application level proposal, this would require a substantial modification of the ETCS equipment, which have already been certified by the industry to fulfil the strict requirements of the SIL4 standard. Thus, all modifications should focus on communication protocols. Nowadays, the industry is considering the migration of the ERTMS from GSM-R towards new wireless access technologies due to the end of support of GSM technology in 2020 and the inefficient use of the spectrum of this technology [START_REF] Ruesche | The European Switch. Packet-Switched Approach to a Train Control System[END_REF], [START_REF] Sniady | An overview of GSM-R technology and its shortcomings[END_REF]. All new proposals, such as GPRS and LTE, are based on packet switching technologies (UIC, 2012), (Tentea, 2011) and TCP/IP protocol stack. This migration to TCP/IP protocols allows us to apply different multipath strategies, which are already used in IP networks, and the packet-switching technologies provide a more efficient resource usage.

Taking all this into account, the MCS proposal has been designed to cover the scenarios of railway communications that will be deployed in the near future. And, thus, the MCS will only consider TCP/IP as a communication protocol stack. In order to take full advantage of MCS, trains need to be multihomed or, in other words trains need to have two or more IP interfaces. This brings into focus the necessity to study protocols [START_REF] Gladisch | Survey on Mobility and Multihoming in Future Internet[END_REF] like HIP, SHIM6, MPTCP or, SCTP. The availability of Software Defined Networking technologies which make a clear separation between data path and signalling control planes and the availability of several protocols, such as OpenFlow and ForCES, present the possibility to use a technology agnostic mechanism to control the MCS. It is important to remark that the MCS is unaware of the wireless technologies used, which provides a great flexibility (WiMAX (Aguado et al., 2013), LTE [START_REF] Tingting | A high-speed railway mobile communication system based on LTE[END_REF], Tetra-IP…). Depending on the number and type of interfaces the communication will be more reliable against EM attacks. Indeed, in order to provide a better resilience against EM attacks it would be advisable to deploy multiple interfaces using different technologies that ideally use different range of frequencies. However, the use of multiple interfaces of the same technology may be also useful as a fallback solution against eventual equipment failures. Thanks to the MCS, the ETCS application of the RBC and train are unaware of the multiple paths since underlying networking protocols manage them transparently to offer a more robust and efficient data transport service to the upper layers of the network stack. This is a great advantage, because no change in the ETCS application protocol is required.

Conclusion

The railway is an attractive target for security attacks, because of its familiarity, ease of access and openness. Among the potential security attacks, electromagnetic terrorism aims at failing equipment or devices which serve the efficiency and safety of the railway transport system. This paper has presented some work performed on the SECurity of Railways against Electromagnetic aTtacks FP7 project aiming to assess the real risks concerning EM attacks on rail networks. It has concentrated on the ground to train communication by presenting the impact of a low-power jammer on this communication. Elements regarding the development of a detection system able to discriminate jamming signals from "normal" electromagnetic railway environmental conditions were also discussed. Finally, using this input, the resilient architecture discussed offers an improved resilience against EM attacks in order to protect wireless communications between trackside and trains principally. The proposed architecture provides a dynamic protection solution combining resilient communication architecture with a resilient health and attack management subsystem.

Figure 1 :

 1 Figure 1: Cells deployed along a GSM-R network

Figure 3 :

 3 Figure 3: GSM-R downlink frequencies -span = 10 MHz (on a station platform).

Figure 5 :

 5 Figure 5: Detection System (DS) and Multipath Communication System (MCS) for the resilient architecture

Figure 6 :

 6 Figure 6: Architecture of Detection System (DS)

Figure 7 :

 7 Figure 7: Architecture of the Multipath Communication System

Acknowledgements

The research leading to these results has received funding from the European Community's Framework Program FP7/2007-2013 under grant agreement n°"285136".


Guillaume Maîtrejean
email: guillaume.maitrejean@ujf-grenoble.fr

Amine Ammar
email: amine.ammar@ensam.eu

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Miroslav Grmela
email: miroslav.grmela@polymtl.ca

Deterministic solution of the kinetic theory model of colloidal suspensions of structureless particles

Keywords: Colloidal suspensions, Smoluchowski equation, Proper generalized decomposition

A direct modeling of colloidal suspensions consists of calculating trajectories of all suspended objects. Due to the large time computing and the large cost involved in such calculations, we consider in this paper another route. Colloidal suspensions are described on a mesoscopic level by a distribution function whose time evolution is governed by a Fokker-Plancklike equation. The difficulty encountered on this route is the high dimensionality of the space in which the distribution function is defined. A novel strategy is used to solve numerically the Fokker-Planck equation circumventing the curse of dimensionality issue. Rheological and morphological predictions of the model that includes both direct and hydrodynamic interactions are presented in different flows.

Introduction

Colloidal systems are encountered in many natural and industrial processes such as blood, foodstuffs, or paints where the presence of particles modifies in a non-linear way the behavior of the fluid. These suspended particles are subjected to different kind of forces coming from the fluid drag, inter-particles potentials, particlessolvent interaction, Brownian forces, etc. Thus, usual descriptions lie in discrete description of each particle and usual models proceed within the Brownian dynamics (BD) or Stokesian dynamics (SD) [START_REF] Brady | Stokesian dynamics[END_REF] framework. When the solvent is considered as continuum and the Reynolds number of particles is small, SD is able to account for interactions between colloidal suspension mediated by the solvent, the socalled HI. As soon as the position and eventually the momentum of all the colloidal particles are known, the virial stress could be computed without any major difficulty. However, the tracking of all the particles in the population implies significant computing time in order to limit the inherent statistical noise that stochastic simulations induce [START_REF] Brady | Computer simulation of viscous suspensions[END_REF]).

An alternative to the BD and SD lies in deriving a continuous and deterministic kinetic theory description. This route leads to the so-called Fokker-Planck equation that characterizes the microstructure through the phase space distribution function (x, R, V). On a more macroscopic level, the Smoluchowski description replaces the Fokker-Planck's formalism by a pair correlation function g(x, R) related to (x, R, V) by

g(x, R) = dV (x, R, V) (1)
where x is the position vector, R the vector connecting two particles, and V their relative velocity.

The difference between the Fokker-Planck and the Smoluchowski formulation lies in the elimination in g(x, R) of fast processes related to inertia (motion of the suspended particles through the fluid). By the simplifying assumptions made for the definition of the Smoluchowski equation, where only binary hydrodynamic interactions are considered, we may expect that the suspensions are rather dilute as the higher-order correlations are neglected. Furthermore, we would like to emphasize that colloidal systems have been widely studied (see, for example, reviews from [START_REF] Stickel | Fluid mechanics and rheology of dense suspensions[END_REF] and more recently [START_REF] Morris | A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow[END_REF]) from a theoretical [START_REF] Batchelor | Effect of Brownian-motion on bulk stress in a suspension of spherical-particles[END_REF][START_REF] Brady | Microstructure of strongly sheared suspensions and its impact on rheology and diffusion[END_REF][START_REF] Lionberger | A Smoluchowski theory with simple approximations for hydrodynamic interactions in concentrated dispersions[END_REF][START_REF] Stickel | A constitutive model for microstructure and total stress in particulate suspensions[END_REF][START_REF] Wilson | An analytic form for the pair distribution function and rheology of a dilute suspension of rough spheres in plane strain flow[END_REF] or experimental [START_REF] Gao | Direct investigation of anisotropic suspension structure in pressure-driven flow[END_REF]) point of view as well as numerically [START_REF] Bergenholtz | The non-Newtonian rheology of dilute colloidal suspensions[END_REF][START_REF] Foss | Structure, diffusion and rheology of Brownian suspensions by Stokesian dynamics simulation[END_REF][START_REF] Phung | Stokesian dynamics simulation of Brownian suspensions[END_REF]). However, the major part of this work addressed concentrated system which cannot be accurately reproduced by Smoluchowski equation as mentioned previously. Nonetheless as we will see in this paper, our approach gives relevant results when comparing to theory. The Smoluchowski formulation is defined on both physical space x and configuration space R; thus, a strategy for solving high-dimensional problems is needed. For that reason, in the literature, the problem is reformulated in a more macroscopic description where g(x, R) is replaced by its second moment m(x, R) as it has been suggested by Phan-Thien in [START_REF] Phan-Thien | Constitutive equation for concentrated suspensions in Newtonian liquids[END_REF]. Let n be the unit vector

n = R R (2)
where R = |R|, the second moment writes:

m αβ = R 2 4π S n α n β g dn (3)
with S the surface of sphere of radius R. However, the equation that governs the evolution of the second-order moment involves the fourth order and so on. Thus, a closure relation relating the fourth-order moment as a function of the lower-order moments is compulsory. This relation constitutes the so-called closure relation.

The contribution of this paper is to solve numerically the pair correlation function without reducing it to its second moment avoiding by the way the use of any closure relation needed to solve the equation governing the evolution of the second moment m(x, R) (see [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF]. We shall not contribute in any way to the formulation of micro-hydrodynamic suspensions model: the pair correlation function is taken as a micro-macro model and is solved using advanced model reduction techniques based on a separated representation, the so-called Proper Generalized Decomposition (PGD) that works efficiently in high dimensional spaces by circumventing the curse of dimensionality that such models imply (cf. "Tensor form of the discrete problem" section and "Appendix 2").

Governing equations of the model

Smoluchowski's equation

Following [START_REF] Smoluchowski | Über Brownsche Molekularbewegung unter Einwirkung äußerer Kräfte und den Zusammenhang mit der verallgemeinerten Diffusionsgleichung[END_REF], we describe the microstructure of the suspension by the pair correlation function g(x, R).

The suspensions under consideration are assumed to be semi-dilute, isothermal, monodisperse, and incompressible. The suspended particles are rigid spheres, all identical, with radius a L, where L is the characteristic length of the fluid container. Both a and L are assumed to be constants. The pair correlation function is highly influenced by the flow kinematics, i.e., influenced by both the type and the strength of the flow. Thus, we allow the suspensions to be spatially inhomogeneous in the sense that the pair correlation function g is a function of R and also of x, but we assume that the density n p of the suspended spheres remains a constant independent of both x and the time t. It is well-known (see, e.g., [START_REF] Morris | A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow[END_REF] and [START_REF] Miller | Suspension flow modeling for general geometries[END_REF]) that in flows with varying shear rate, the spheres migrate and thus n p becomes a function of both x and t. It is therefore necessary to adopt n p (x, t) as an additional state variable joining the velocity v(x, t) and the pair correlation function g(x, R, t). We intend to discuss rigid sphere suspensions in such extended setting in a future paper. In this paper, we consider flows with constant shear rates (in "Steady-state solution in a simple shear flow" section) for which the extended setting is not necessary, but we also investigate spatially varying shear rates (in "Computing the pair correlation function in a simple recirculating flow" section) for which the setting with n p = const. provides only an approximative model in which the migration of the spheres is ignored. How much will the results obtained in "Computing the pair correlation function in a simple recirculating flow" section change if the migration is taken into account remains an open question that will be addressed in future works.

The Smoluchowski equation [START_REF] Smoluchowski | Über Brownsche Molekularbewegung unter Einwirkung äußerer Kräfte und den Zusammenhang mit der verallgemeinerten Diffusionsgleichung[END_REF]) (see also [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF] and references cited therein) writes:

∂g ∂t = - ∂ ∂ x α (gv α )- ∂ ∂ R α gR β ∂v α ∂ R β + ∂ ∂ R γ gλ γ αβ D αβ + ∂ ∂ R α g αβ ∂ ∂ R β (U + k B T ln g) (4)
By v(x), we denote the overall velocity of the suspension and D denotes the rate of strain tensor,

D αβ = 1 2 ∂v α ∂ x β + ∂v β ∂ x α (5)
and the tensors and λ are kinetic coefficients which are functions of R (see "Appendix 1"). The total interparticle potential U(R) writes:

U = U HC + U nl (6)
with

U HC (R) = 0, for R > 2a ∞, for R ≤ 2a. (7)
and

U nl (R) = R 12 , for R ≥ 2a. (8)
where is a positive constant. g eq denotes the equilibrium correlation function as

g eq (R) = exp U nl (R) k B T (9)
Equation 4 is subjected to two boundary conditions:

-When R → ∞, g(x, R) → n 2 p . -When R = 2a, the flux through the boundary vanishes.

According to [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF], the dimensionless form of the equation governing the evolution of the correlation function writes:

∂g ∂t = -Pe ∂ ∂ x α (gv α ) -Pe ∂ ∂ R α R β ∂v α ∂ x β g + Pe ∂ ∂ R γ gλ γ αβ D αβ + ∂ ∂ R α g eq αβ ∂ ∂ R β ∂g ∂g eq ( 10
)
where

x → x/a, R → R/a g → g n 2 p , U nl → 1 k B T U nl
and the Péclet number defined as

Pe = γ 6πη 0 a 3 k B T
with γ the strain rate and η 0 a positive coefficient. g(R) = g(R)g eq (R) denotes the deviation from equilibrium of the correlation function.

For more details concerning definition of kinetic coefficients, please refer to "Appendix 1" as well as [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF] and the references therein.

Expression of the stress tensor

The stress tensor groups three contributions:

σ = σ hydr + σ DI + σ HI ( 11
)
σ hydr is the contribution from the Newtonian fluid in which particles are suspended. σ DI comes from the direct interactions including hard core and long-range potential. σ HI is a contribution arising due to the hydrodynamic interactions.

The resulting dimensionless form of the extra stress tensor reads [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF])

σ αβ = -1 + 2 5 φ D αβ - 27 8π φ 2 Pe × 2 S 2 dn n α n β - 1 3 δ αβ × [1 -μ 1 (2)] g(2n) -dR 3(μ 1 -μ 2 ) + R ∂μ 1 ∂ R + (1 -μ 1 )R ∂U nl ∂ R × n α n β - 1 3 δ αβ g(R) ( 12
)
where φ = (4/3)n p πa 3 is the volume fraction of particles in the suspension and S 2 refers to integration over the unit sphere S 2 .

Tensor form of the discrete problem

The evolution equation of the pair correlation function involves both the physical and the configurational coordinates which are up to three dimensions each. If we account for the time, the total dimension of the problem is up to seven. In such high-dimensional spaces, the classical Finite Element (FE) method cannot be applied successfully because mesh-based discretization techniques suffer the so-called curse of dimensionality.

The separated representation involved in the PGD [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part II: Transient simulation using space-time separated representations[END_REF][START_REF] Ammar | Reduced numerical modeling of flows involving liquid-crystalline polymers[END_REF] allows circumventing or at least alleviating the curse of dimensionality encountered in the solution of the equation governing the evolution of the pair correlation function. The main idea of the PGD method is to find a solution of the linear system defined by

AG = B
where A, G, and B are written under separated form as a sum of tensorial products defined on a subspace of the initial domain. Without loss of generality, when choosing to separate the problem onto D subspaces, A, G, and B write, respectively

A = n A i=1 A i 1 ⊗ A i 2 ⊗ ... ⊗ A i D G = n i=1 G i 1 ⊗ G i 2 ⊗ ... ⊗ G i D B = n B i=1 B i 1 ⊗ B i 2 ⊗ ... ⊗ B i D
where G is the discrete tensor form of g. The PGD method is an a priori method, opposed to the a posteriori method (as the well-known Proper Orthogonal Decomposition method (POD)) and does not need information concerning the solution to compute the reduced approximation basis. Furthermore, PGD method can be used with any partial differential equation, written in continuous or discrete form, as long as the operator A and the right-hand side term B may be written under separated form. For the pair correlation function equation, numerous separation strategies are feasible, for example:

-the functions may be defined on subspaces of dimension 1 only, such as

G = n i=1 G i x ⊗ G i y ⊗ G i z ⊗ G i u ⊗ G i v ⊗ G i w ⊗ G i t (13)
where time coordinate is t, space coordinates are x = (x, y, z) T , and configurational coordinates R = (u, v, w) T , -it is possible to separate the solution over physical, configurational, and time spaces

G = n i=1 G i x ⊗ G i R ⊗ G i t ( 14
)
-or isolate the physical and configurational spaces with a classical time discretization

G t = n i=1 G i x ⊗ G i R (15)
where G t corresponds to the solution G at time t.

We chose to separate the problem on two subspaces, the physical and configurational spaces, and to treat the time space using classical discretization (an implicit Euler scheme) because it is the most convenient way to deal with the influence induced by the presence of the particles on the fluid evolution.

Let introduce the tensor form related to the Eq. 10 which can be rewritten as:

∂g ∂t + v • ∂g ∂ x + E 0 g + ∂ ∂ R (E 1 g) - ∂ ∂ R E 2 ∂ ∂ R g g eq = F (16)
where

E 0 (x, R) = Pe ∂v γ ∂ x γ + ∂ R β ∂ R α ∂v α ∂ x β -Pe ∂λ γ αβ ∂ R γ D αβ + λ γ αβ ∂ D αβ ∂ R γ E 1 (x, R) = PeR β ∂v α ∂ x β -Pe λ γ αβ D αβ E 2 (x) = g eq λ αβ F(x, R) = Pe ∂ R β ∂ R α ∂v α ∂ x β g eq + R β ∂v α ∂ x β ∂g eq ∂ R α
Some of the above coefficients involve both the physical and the configurational coordinates, but their dependence on the physical space lies only through the velocity gradient. To rewrite the previous coefficients in a separate form, we chose to express their contribution in the velocity gradient canonical basis (see "Appendix 3" for more details concerning the procedure).

Thus E 0 (x, R) can be rewritten as:

E 0 (x, R) = ∂v α ∂ x β (x) E αβ 0 (R) (17)
We follow the same procedure for the other coefficients E 1 (x, R) and F(x, R) leading to the weak form of Eq. 16:

×I g * ∂g ∂t + v • ∂g ∂ x + E 0 g + ∂ ∂ R (E 1 g) - ∂ ∂ R E 2 ∂ ∂ R g g eq dx dR = ×I g * F dx dR (18)
where the variables x ∈ and R ∈ I are, respectively, the position vectors of the physical and configurational space. In this formulation, the second derivative term is rewritten as follows:

×I g * ∂ ∂ R E 2 ∂ ∂ R g g eq dx dR = × g * E 2 ∂ ∂ R g g eq dx dR - ×I ∂g * ∂ R E 2 ∂ ∂ R g g eq dx dR (19)
where is the integration over the boundary of I. The problem is subjected to two boundary conditions for R → 2 and R → ∞. For the far-field limit, R → ∞, we assume that the pair correlation function tends to the uniform distribution, i.e., g = 1. At the contact surface of a pair (R → 2), we adopt the boundary condition introduced by [START_REF] Batchelor | Effect of Brownian-motion on bulk stress in a suspension of spherical-particles[END_REF]: a zero flux condition. It means that the flux vanishes on R = 2:

n α Pe -R β ∂v α ∂ x β + λ γ αβ D αβ g + E 2 ∂ ∂ R β g g eq = 0 (20)
with n α = R α /R. Thus, on the boundary R = 2:

n α E 2 ∂ ∂ R β g g eq = n α ⎡ ⎢ ⎢ ⎢ ⎣ Pe R β ∂v α ∂ x β -λ γ αβ D αβ E 1 g ⎤ ⎥ ⎥ ⎥ ⎦ (21)
Thus, the second derivative term writes:

×I g * ∂ ∂ R E 2 ∂ ∂ R g g eq dx dR = × g * (n • E 1 ) g dx dR - ×I ∂g * ∂ R E 2 1 g eq ∂g ∂ R + g ∂ ∂ R 1 g eq dx dR (22)
Finally, the weak form reads:

×I g * ∂g ∂t + g * v • ∂g ∂ x + g * E 0 g + g * E 1 ∂g ∂ R + ∂g * ∂ R E 2 ∂g ∂ R + ∂g * ∂ R E 3 g dx dR - × g * (n • E 1 ) g dx dR = ×I g * F dx dR (23) with E 0 (x, R) = Pe ∂v γ ∂ x γ + ∂ R β ∂ R α ∂v α ∂ x β -Pe ∂λ γ αβ ∂ R γ D αβ + λ γ αβ ∂ D αβ ∂ R γ (24) E 1 (x, R) = Pe R β ∂v α ∂ x β -λ γ αβ D αβ (25) E 2 (R) = αβ (26) E 3 (R) = αβ g eq ∂ ∂ R 1 g eq (27) F(x, R) = Pe ∂ R β ∂ R α ∂v α ∂ x β g eq + R β ∂ ∂ R α ∂v α ∂ x β g eq +R β ∂v α ∂ x β ∂g eq ∂ R α ( 28
)
Equation 23 is discretized using a standard Finite Element approximation in each coordinate. In homogeneous (rheometric) flows, no specific strategy is needed to solve the pair correlation function. However, in complex flows (i.e., when the velocity gradient is no longer homogeneous), we use reduced modeling to alleviate difficulty related to the high-dimensional domain of definition of the pair correlation function. For this purpose, we introduce N and M that represent the vectors containing the shape functions associated, respectively, with each coordinate x and R.

By defining the following matrices:

A 1 x = Nv dN T dx dx A 1 R = I MM T dR A 2 x = NN T dx A 2 R = I dM dR E 2 dM T dR dR A 3 x = NN T dx A 3 R = I dM dR E 3 M T dR A 3(α-1)+β+3 x = N dv α dx β N T dx A 3(α-1)+β+3 R = I ME αβ 0 M T + ME αβ 1 dM T dR dR -M n • E αβ 1 M T dR B 3(α-1)+β x = dv α dx β N dx B 3(α-1)+β R = I F αβ M dR
the discrete form results:

AG = B ( 29
)
with

A = 12 i=1 A i x ⊗ A i R ( 30
)
and

B = 9 i=1 B i x ⊗ B i R (31)
Note that this equation is a tensorial equation. The unknown field G contains the tensor product of the two functions given in Eq. 15. Thus, the resulting problem becomes non-linear and an appropriate linearization technique must be applied (see "Appendix 2").

Numerical results

The pair correlation function is defined both on the physical space and configurational space, up to three dimensions each. In what follows, we choose to represent the configurational space either in two or three dimensions. By doing this, we can compare the results obtained with the ones given in [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF] and then evaluate the impact of the closure introduced in [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF]. In this case, the physical space disappears and Eq. 18 is solved, which is a classical Finite Element formulation which does not require any special numerical strategy. According to Zmievski's results (see [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF], the suspension reach an equilibrium state for R > 5. Thus, we choose to represent the configuration space by a hollow disk and a hollow sphere, respectively, in two and three dimensions, with an internal radius of 2 and the external one of 5.

Steady

By solving the pair correlation function instead of reducing it to the second-order tensor m, it is clear that we retain more information on the microstructure. We can see on Figs. 1 and 2 the orientation of the microstructure revealed by the angular dependence of g. We observe that the anisotropy decreases as R increases, i.e., g tends to 1, which corresponds with the prescribed boundary condition g | R=5 = 1. On the other hand, for small R, the anisotropy of the microstructure is significant and orientated in the direction of the velocity gradient in perfect agreement with the Zmievski 's results [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF]. Notice in Figs. 1 and2 the magnitude of the predicted anisotropy which is quite different when considering 2D or 3D configuration spaces whereas the distribution on the plane of shear is similar.

N 1 = Pe -1 (σ 22 -σ 11 ) ( 33
)
and

N 2 = Pe -1 (σ 33 -σ 22 ) ( 34
)
with σ αβ the components of the dimensionless extrastress tensor, we can notice the error induced by the closure relation used by [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF] for calculating the microstructure tensor.

From a physical point of view, the thinning behavior observed for the first normal stress coefficient is in accordance with experimental measurements (Kolli Fig. 8 Evolution of g along a trajectory for FE calculation et al. 2002;[START_REF] Ohl | The characterization of the steady-state shear and normal stress functions of highly concentrated suspensions formulated with viscoelastic liquids[END_REF][START_REF] So | Microstructure evolution and theological responses of hard sphere suspensions[END_REF], as well as the negative value of the second stress coefficient. However, for large Pe, thickening behavior which has been observed experimentally [START_REF] Kolli | Transient normal stress response in a concentrated suspension of spherical particles[END_REF][START_REF] Ohl | The characterization of the steady-state shear and normal stress functions of highly concentrated suspensions formulated with viscoelastic liquids[END_REF][START_REF] So | Microstructure evolution and theological responses of hard sphere suspensions[END_REF]) and numerically [START_REF] Brady | Computer simulation of viscous suspensions[END_REF][START_REF] Foss | Structure, diffusion and rheology of Brownian suspensions by Stokesian dynamics simulation[END_REF] in concentrated suspensions cannot be reproduced by the model here considered.

Computing the pair correlation function in a simple recirculating flow

We compute the evolution of the pair correlation function in the recirculating flow depicted in Fig. 7. To evaluate the accuracy of the PGD solution, we compare it with the one computed by applying the Finite Element method.

The flow is defined by the velocity field

v = -y x 2 + y 2 x x 2 + y 2 . ( 35
)
Because Finite Element calculations fail when considering the whole physical and conformational domains, we adopt the following strategy: g is computed only along a streamline (which is a perfect circle in the flow here addressed) taking into account that the equation governing the pair correlation function g is purely advective in the physical space. The steady-state solution is then computed from an isotropic state after several rotations (cf. Fig. 8).

Note that the FE solution is computed at a single position (in the present case, it suffices because the problem symmetry) whereas the PGD solution is computed for the whole physical space. The problem is discretized by considering 3,905 and 1,680 degrees of freedom, respectively, for the physical (two-dimensional) and conformational (twodimensional) domains. Figures 9 and 10 depict the results obtained, respectively, with the FE and PGD method after 50 enrichments (that is, the PGD solution consists in 50 modes).

As it can be noticed, the results shown in Figs. 9 and 10 are very close, and we need to reduce the pair correlation function to its second-order moment tensor in order to appreciate the difference between FE and PGD solutions. Figure 11 shows the trace of the second-order moment of the pair correlation function computed at position (-0.51, 0.33) as a function of the number of terms involved in the separated representation of the pair correlation function. The solution computed by applying the PGD method even with few enrichments (of around ten terms in the separated representation) is close to the reference FEM solution even if 40 enrichments are required to achieve the desired accuracy.

Figure 12 shows the five first functions defined on each subspace involved in the separated representation of the solution.

In this geometry, the fluid velocity field is assumed known, and we look for the steady state of the pair correlation function, allowing us to perform some numerical test illustrating the PGD capabilities. First, we represent in Fig. 13 the error of the PGD solution as where Sol FE and Sol PGD are, respectively, the solutions obtained with FE and PGD method.

The influence of the refinement on the convergence in not really significant as the error reached after 80 iterations remains almost unchanged, although the finer mesh has twice more nodes than the coarser one. The optimality of the separated representation, that is, the number of modes that it involves, constitutes a work in progress [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF].

When plotting the calculation time of the FE and PGD method as a function of the total number of degrees of freedom (see Fig. 14), we can notice the efficiency of the PGD for treating high-dimensional problem. These results have been obtained on a standard laptop with Intel Core i7 @ 2-GHz Processor with 6 GB of RAM.

Solving the pair correlation function in complex flows

For illustrating the performances of the PGD for addressing complex flows, in what follows we consider the driven cavity flow problem. We consider a suspension characterized by φ = 0.05 with long-range particle interactions ( = 10 6 in Eq. 8).

In addition to the equation governing the pair correlation function that we considered in the previous sections, now we should consider also the ones related to the fluid flow kinematics that involve the mass and momentum balance equations:

∂v α ∂t = -Pe ∂ ∂ x β v β u α + I αβ p + 1 Re σ αβ (37) ∇ • v = 0 ( 38
)
where u = v/ρ, σ is defined in Eq. 11 and Re the Reynolds number defined by Re = ρa 2 γ η 0

The physical (2D) and the conformational (2D) domains are discretized by employing 1,976 and 1,680 nodes. At t < 0, the fluid is at rest, i.e., v = 0, and at t = 0, we impose a velocity on the upper side defined by:

v(x) = 16 h 4 v max x 2 (h -x) 2 (39)
with the width of the container h = 1 and v max = 1.

We perform the calculation for ten time steps with t = 0.01. The resolution of the non-linear equation 37 is performed using a Newton-Raphson procedure, and 40 terms of the separated representation in the PGD framework were calculated at each time step (this number was enough for reaching the desired accuracy). For sake of simplicity in what follows, the word iteration refers to time step.

The fluid velocity field is depicted in Figs. 15 and16, respectively, at iterations 1 and 10.

Figure 17 shows the five first functions defined on each subspace involved in the separated representation of the solution G at iteration 10.

We evaluate the evolution of the pair correlation function for the three different nodes in the physical space depicted in Fig. 18 with the corresponding velocity gradient: As we can notice in Figs. 19, 20, and 21, the more is the velocity gradient, the higher is the anisotropy of the microstructure. For the node 1 that is located in a zone exhibiting large shear rates, flow-induced anisotropy is noticed from the first time step, whereas the distribution in the other nodes remains essentially isotropic.

Figures 22 and 23 depict the effects of the suspension on the stress at iteration 1, 5, and 10.

Conclusion

The most direct approach to investigations of colloidal suspensions is the approach adopted in stochastic simulations. Smoluchowski's mesoscopic theory provides an alternative setting. The microstructure in the Smoluchowski theory is characterized by the pair correlation function. In this paper, the PGD was used for solving the Smoluchowski equation. Because the PGD numerical complexity scales linearly with the dimension of the space instead of the exponential growing characteristic of mesh-based discretizations, PGD allows circumventing the curse of dimensionality that models defined in high-dimensional spaces (like the ones related to the Smoluchowski equation) exhibit. Thus, by applying the PGD, we can avoid the solution of the equation governing the evolution of the second-order moment related to the pair correlation function and then the issues related to the closure approximation that such an approach involves.

In this paper, we proved the capabilities of simulating complex flows within the kinetic theory framework. However, we are conscious that the model here considered is not the more pertinent for addressing scenarios of real interest because such scenarios imply taking into account particle migration as well as other ignored effects essentially in concentrated regimes. These topic will be addressed in future works.

Appendix 1: Kinetic coefficients

The kinetic coefficients and λ derive from the Stokes problem (motion of rigid spheres in a Newtonian fluid) and appear in Eq. 4. They are defined in [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF] by:

αβ = k B T 3πη 0 a σ 1 (R/a) R α R β R 2 +σ 2 (R/a) δ αβ - R α R β R 2
(40)

λ αβγ = μ 1 (R/a) R α R β R 2 + μ 2 (R/a) δ αβ - R α R β R 2 R γ (41)
The dimensionless form of the above terms is given by

αβ → 3πη 0 a k B T αβ and λ αβγ → 1 a λ αβγ
In this paper, we strictly follow the approximations of the coefficient proposed by [START_REF] Zmievski | Nonlinear microstructure and rheology of semidilute colloidal suspensions of structureless particles[END_REF] (based on [START_REF] Batchelor | Effect of Brownian-motion on bulk stress in a suspension of spherical-particles[END_REF] and [START_REF] Batchelor | The hydrodynamic interaction of two small freely-moving spheres in a linear flow field[END_REF]):

σ 1 ≈ 1 - 3 2 a R + a 3 R 3 - 15 4 a 4 R 4 + 109.0 a 6 R 6 -236.0 a 7 R 7 (42) σ 2 ≈ 1 - 3 4 a R - 1 2 a 3 R 3 -10.336 a 6 R 6 (43) μ 1 ≈ 5 a 3 R 3 -8 a 5 R 5 + 25 a 6 R 6 -
A G = B
where

A = n A j=1 A j 1 ⊗ A j 2 ⊗ ... ⊗ A j N B = n B j=1 B j 1 ⊗ B j 2 ⊗ ... ⊗ B j N
We are looking for the solution as a tensor product of functions

G = ∞ j=1 G j 1 ⊗ G j 2 ⊗ ... ⊗ G j N
For computing the involved functions, we are minimizing the residual

J (G) = A G -B 2
that could be rewritten

J (G) = G -A -1 B T A T A G -A -1 B = G -A -1 B A T A
Thus, by differentiation in relation to the field G, the new problem to be solved writes

∂GA T A G = ∂GA T B
The adaptive strategy consists of enriching progressively the separated representation. Thus, at iteration n G + 1, we look for:

G = n G j=1 G j 1 ⊗ G j 2 ⊗ ... ⊗ G j N G F + R 1 ⊗ R 2 ⊗ ... ⊗ R N G R
Because of the non-linearity, these functions are computed by using an alternating directions fixed point algorithm that proceeds by calculating R j assuming all the others are known. Thus, when we are looking for R j , the following test function is used

G * = R 1 ⊗ ... ⊗ R j-1 ⊗ R * j ⊗ R j+1 ⊗ ... ⊗ R N
The system to be solved then writes

G * A T AG R + G * A T AG F = G * A T B that results in n F i=1 n A k =1 n A k=1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ A k j A k j G i j N h = 1 h = j R T h A k h A k h G i h ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ + n A k =1 n A k=1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ A k j A k j R j N h = 1 h = j R T h A k h A k h R h ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = n A k =1 n B k=1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ A k j B k j N h = 1 h = j R T h A k h B k h ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
Appendix 3: Writing coefficient under separate form Some of the coefficients of the pair correlation function are defined both on physical and conformational space, and following the PGD algorithm (see "Appendix 2"), it is needed to rewrite these coefficients in a separated form. Let A(x, R) be such a coefficient :

A(x, R) = i B i (x)C i (R) (46)
However, it is not easy to write the pair correlation function coefficients in a separate form. As the only dependence of these coefficients on the physical coordinates is through the velocity gradient, we could express it in the velocity gradient canonical basis. In that basis, the velocity gradient writes (in 2D for the sake of simplicity):

∂v ∂ x = ∂v 1 ∂ x 1 ∂v 1 ∂ x 2 ∂v 2 ∂ x 1 ∂v 2 ∂ x 2 = ∂v 1 ∂ x 1 1 0 0 0 + ∂v 1 ∂ x 2 0 1 0 0 + ∂v 2 ∂ x 1 0 0 1 0 + ∂v 2 ∂ x 2 0 0 0 1 (47)
A(x, R) can be rewritten as:

A(x, R) = ∂v 1 ∂ x 1 (x)A 11 (R) + ∂v 1 ∂ x 2 (x)A 12 (R) + ∂v 2 ∂ x 1 (x)A 21 (R) + ∂v 2 ∂ x 2 (x)A 22 (R) = 2 i=1 2 j=1 ∂v i ∂ x j (x)A ij (R) (48)
where all the terms A ij (R) are determined numerically.

Fig. 1

 1 Fig. 1 Steady-state pair correlation function in simple shear flow in a 2-D configuration space. The Péclet number is set to 1 and = 0

  v vs u.(c) w vs u.(d) w vs v.

Fig. 2 (Fig. 3

 23 Fig. 2 Steady-state pair correlation function in simple shear flow in a 3-D configuration space with vector R = (u v w) T . The Péclet number is set to 1 and = 0

Fig. 4 Fig. 5 Fig. 6 Fig. 7

 4567 Fig. 4 g(R) on the shear plane for the results obtained by Morris and Katyal (2002) (a) and with the Smoluchowski model (b) with a 3D configurational space

Fig. 9

 9 Fig. 9 Configurational space obtained with FE method

Fig. 10

 10 Fig. 10 Configurational space obtained with PGD method after 50 iterations

Fig. 11

 11 Fig. 11 Trace of the second-order moment of the pair correlation function computed at position (-0.51, 0.33) as a function of the number of terms involved in the separated representation of the pair correlation function

Fig

  Fig. 12 Five first functions involved in the PGD representation defined on the physical (top) and configurational space (bottom)

Fig. 13

 13 Fig. 13 Error of the PGD method as a function of the number of terms m involved in the separated representation of the pair correlation function and for three different refinements of the configurational space

Fig. 14

 14 Fig. 14 Calculation time for both FE and PGD methods as a function of the total number of degrees of freedom used for describing the problem

Fig. 15

 15 Fig. 15 Fluid velocity field and associated streamline after the first iteration

Fig. 16

 16 Fig. 16 Fluid velocity field and associated streamline at iteration 10

FigFig. 18

 18 Fig. 17Five first functions involved in the PGD separated representation defined on the physical (top) and configurational space (bottom)

  Iteration 1. (b) Iteration 5. (c) Iteration 10.

Fig. 19 (

 19 Fig. 19 Configurational space for node 1 at iterations 1, 5, and 10

Fig. 20 (

 20 Fig. 20 Configurational space for node 2 at iterations 1, 5, and 10

Fig. 21

 21 Fig. 21 Configurational space for node 3 at iterations 1, 5, and 10

Fig. 22 (

 22 Fig. 22 (σ 22 -σ 11 ) at iteration 1, 5, and 10

Fig. 23 σ

 23 Fig. 23 σ 12 at iteration 1, 5, and 10

  This appendix describes the general PGD framework for solving a model defined in a space of dimension N whose discrete form reads:

				423.712	a 7 R 7 + 907.424	a 8 R 8
						(44)
	μ 2 ≈ -	16 3	a 5 R 5 -15.317	a 6 R 6		(45)
	Appendix 2: PGD algorithm