John Carran

Rachel Waschbüsch

Philippe Savignac

Carbanionic Displacement Reactions at Phosphorus : Diethyl (2-Pyridyl)methylphosphonate Synthesis *

Keywords: 2-picoline, LDA, diethyl chlorophosphate, Wittig-Horner reaction

We describe the formation and reactions of the α-lithiated (2-pyridyl)methylphosphonate 1. The lithiated α-picoline is obtained by metallation in THF at low temperature with LDA (lithium diisopropylamide) (2 equiv.). This is then condensed with diethyl chlorophosphate. The thus-formed stable carbanion 1 reacts with various electrophiles: H 2 O, D 2 O, TMSCl, RX, and aldehydes in a Wittig-Horner reaction.

Introduction

The potential value of phosphonates in both biological and industrial applications has stimulated much interest in their synthesis. Currently the formation of the carbon-phosphorus bond can be accomplished according to four main procedures. [START_REF] Engel | Synthesis of Carbon-Phosphorus bonds[END_REF] The first procedure, the Michaelis-Arbuzov reaction, 2 involves a trivalent phosphorus atom which is converted to a pentavalent species by displacement reactions at aliphatic carbon with a trivalent phosphorus ester. The second, the Michaelis-Becker reaction, 3 involves treatment of an aliphatic carbon with a dialkyl phosphite anion. The third, the Kinnear-Perren reaction, 4 involves alkylation of phosphorus halides mediated by Lewis acid, and the fourth involves the conjugate addition of trivalent phosphorus reagents at unsaturated carbon. 5 In connection with our interest in phosphorylated carbanions, we have directed our attention towards a fifth process in which treatment of a pentavalent phosphorus species by an organometallic reagent is involved. This reaction, which readily generates carbon-phosphorus bonds, occurs by carbanionic displacement of halides or ester groups at a relatively electrophilic phosphorus center. [START_REF] Engel | Synthesis of Carbon-Phosphorus bonds[END_REF] While there are a number of reports in the literature relating to the four previous procedures, the carbanionic displacement reactions at phosphorus are largely neglected. The synthetic value of this approach lies in the facile formation of stabilized carbanionic reagents generated by reaction of strong bases with a large variety of readily available starting materials and in the subsequent condensation reaction with trialkyl-or dialkylchlorophosphates. The method described herein involves the use of the stabilized anion of 2-methylpyridine (α-picoline) in reaction with diethyl chlorophosphate. Because of their use in heterocyclic chemistry, pyridines and substituted pyridines have attracted considerable attention. Pyridine derivatives, for example, are useful precursors for the preparation of polymers, pharmacological compounds, [START_REF] Kleeman | [END_REF] liquid crystals [START_REF] Barbera | [END_REF] and ligands for transition metal cations.

Results and discussion

Diethyl (2-pyridyl)methylphosphonate 2 was first obtained by Maruszewska-Wieczorkowska, Michalski and Skowronska from 2-picolyl chloride and triethylphosphite. 9 Later it was used in synthesis by Piechucki under phase-transfer catalysed Wittig-Horner reactions. 10 The previously reported procedure requires the 2-picolyl chloride free-base which is known as a very unstable and expensive starting material. By contrast we have developed a convenient route to 2 in a one-pot reaction which uses easily accessible and stable starting materials: αpicoline and diethyl chlorophosphate. The key intermediate, 2-lithiated picolyl carbanion, was obtained at low temperature by addition of 2-methylpyridine to LDA (2 equiv.) in THF solution. On standing at low temperature we observed the appearance of a deep red-orange colour. After about 30 min the carbanion is condensed with diethyl chlorophosphate in an exothermic reaction. As observed by 31 P NMR the chlorophosphate is completely consumed and the phosphorylated carbanion 1 is stable and can be kept at 0°C for a long time without decomposition (δ 31 P (THF) +40.6). On treatment at 0°C with a saturated solution of NH 4 Cl, we isolate the (2-pyridyl)methylphosphonate 2 in 83-85% yield after distillation. The preparation can be reproduced on large scale with the same success (Scheme 1).

Scheme 1

Attempts to synthesise the 3-and 4-methylphosphonate analogues of 2 with 3-methyl and 4methyl pyridine as starting materials failed to give the expected pyridylmethylphosphonates. Phosphorylation of these compounds under the same conditions as described above resulted

N H 3 C N CH (EtO) 2 P Li O N CH 2 (EtO) 2 P O 1) 2 LDA, THF -78°C 2) (EtO) 2 P(O)Cl H 2 O 1 2
in the formation of compounds 3a and 3b respectively. The increased acidity of the ortho proton of the starting material leads to ortho lithiation and thus to a 5-methyl-2-(diethyl phosphonyl)pyridine (3a) or 4-methyl-2-(diethyl phosphonyl)pyridine (3b) product respectively. Compound 3a is only formed in low yield (20%) and is obtained in a mixture of compounds whereas 3b is formed in 85% yield and is easily isolable in pure form by distillation.

The (2-pyridyl)methylphosphonate 2 may be readily converted into the 1,1-dideuteriated compound 5 under mild basic conditions. Thus, when 2 was metallated at low temperature with LDA (2 equiv.) in THF solution followed by the exothermic condensation of TMSCl at room temperature, the α-lithiated, α-silylated, (2-pyridyl)methylphosphonate 4 (δ 31 P (THF) +39.8) was obtained and readily converted with excess D 2 O at room temperature into the corresponding 1,1-dideuterio-1-(2'-pyridyl)methylphosphonate 5. We have already described a procedure for the conversion of activated methylphosphonates into alkylphosphonates via the introduction in the α-position of a transient trimethylsilyl group. 11 With activated phosphonates there is often a serious risk that the alkylphosphonate formed could undergo further alkylation. The trimethylsilyl moiety normally acts as a protecting group allowing the introduction of only one alkyl group thus avoiding the precautions needed in order to exclude the formation of dialkyl compounds resulting from acid-base equilibrium. The trimethylsilyl group is then easily removed by treatment with EtOH / LiOEt at 0°C. However, utilising this methodology we found that, although the desired monoalkyl product was formed, it only occured as a component of a mixture of monoalkyl-, dialkyl-, and starting product, and was thus difficult to isolate by distillation. The dialkyl compound 6 can be easily synthesised in high yield (85% for dimethyl) by the use of two equivalents of base in reaction with 2. We found that for this reaction we obtained best results when one equivalent of LiHMDS (lithium hexamethyldisilazane) and one equivalent of n-BuLi were used together with an excess of alkylating agent. (Scheme 2)

N CH 3 3 a N P(OEt) 2 O 3 b CH 3 (EtO) 2 P O Scheme 2
The α-phosphorylated carbanion 1 is ideally suited as a reagent in Wittig-Horner reactions. [START_REF] Johnson | Ylides and Imines of Phosphorus[END_REF] We have developed a particularly interesting and convenient one-pot procedure for converting α-picoline into the vinylpyridine derivatives [START_REF] Windsceif | [END_REF] 7 in good overall yield for the two steps starting from α-picoline without isolation of the intermediate pyridyl phosphonate 2 (Scheme 3, Table I).

Scheme 3

Experimental NMR spectra were recorded on a Bruker AC 200 spectrometer operating at 200 MHz for proton, 50.3 MHz for carbon and 81.01 MHz for phosphorus. 31 P downfield shifts (δ) are expressed with a positive sign, in ppm, relative to external 85% H 3 PO 4 in H 2 O. 1 H and [START_REF] Windsceif | [END_REF] C chemical shifts (δ) are reported in ppm relative to CDCl 3 as internal standard. Coupling constants (J ) are given in Hz. The following abbreviations are used: s, d, t, q, m for singlet, doublet, triplet, quadruplet and multiplet respectively. Organic solvents were purified by standard procedures. THF was distilled under an inert atmosphere from purple solutions of sodium:benzophenone ketyl. The synthesis of all compounds were carried out under dry nitrogen. To a solution of n-BuLi (66 ml of a 1.59 M soln. in hexane; 0.105 mol) in THF (100 ml) at -20°C is added diisopropylamine (11.6 g; 0.116 mol) dropwise. This solution is stirred for 10 min and then cooled to -78°C. A mixture of α-picoline (4.65 g; 50 mmol) in THF (10 ml) is then added dropwise to the freshly prepared LDA. A red-orange colour is observed. The reaction mixture is stirred at -78°C for 30 min before a solution of diethyl chlorophosphate (9.50 g; 55 mmol) in THF (15 ml) is slowly added via a dropping funnel. The mixture is stirred at -78°C for 10 min and then slowly allowed to warm to room temperature. A solution of NH 4 Cl in water (30 ml) is then added to the reaction mixture and the phases are separated after stirring for 10 min. The aqueous phase is then extracted with dichloromethane (3 x 50 ml) and the combined organic phases are dried over MgSO 4 prior to evaporation and distillation by bulb-to-bulb apparatus to yield the title compound as a yellow oil (9.

4-methyl-2-(diethylphosphonyl)pyridine 3b

To freshly prepared LDA (from 15.15 g (0.15 mol) diisopropylamine and 66.6 ml of BuLi) in THF (50 ml) at -78°C is added a solution of 4-picoline (4.65 g, 50 mmol) in THF (10 ml).

The resulting yellow solution is stirred at -78°C for ten minutes and then allowed to slowly warm to 0°C over a period of 1 h. The reaction mixture is then recooled to -78°C and a solution of diethyl chlorophosphate (9.5 g, 55 mmol)) in THF (10 ml) is added. An exothermic reaction is observed and the reaction mixture is allowed to slowly warm and is stirred for 1 h at room temperature. A saturated solution of ammonium chloride (50 ml) is added to the reaction mixture and the phases are then separated. The aqueous phase is first washed with diethyl ether (3 x 20 ml) and then dichloromethane (1 x 20 ml). The combined organic extracts are dried over MgSO 4 and are then filtered and solvents are evaporated. The resulting orange liquid is then distilled to give the title product as a mobile yellow oil (bp 175-185°C / 0.1 mm / Hg; 9.7 g, 85 %). 31 P-NMR (CDCl 3 ): δ= +11.62.

Diethyl 1,1-dideuterio-1-(2'-pyridyl)methylphosphonate 5

The same procedure as above is employed except that after addition of diethyl chlorophosphate the reaction mixture is warmed to room temperature and then a solution of TMSCl (1.1 equiv.) is added, the reaction is exothermic and the reaction mixture is then stirred at 20°C.

Diethyl 1-methyl-1-(2'-pyridyl)ethylphosphonate 6a

LiHMDS (4.4 mmol) is prepared at -78°C by the addition of a solution of HMDS (0.78 g, 1.1 equiv.) in THF(5 ml) to a solution of n-BuLi (3.03 ml of a 1.45 M soln, 4.4 mmol) in THF (20 ml) at -78°C. After 20 min a solution of diethyl (2-Pyridyl)methylphosphonate 2 (1 g, 4.4 mmol) in THF (10 ml) is added to the reaction mixture at -78°C. The now yellow solution is allowed to slowly warm to RT. The reaction mixture is then cooled to -78°C and n-BuLi (3.03 ml, 4.4 mmol) is added. The reaction mixture is cooled to -50°C and MeI (1.25 g, 2.2 equiv) in THF (5 ml) is slowly added. The mixture is allowed to slowly warm to RT during which time the yellow colour disappears. NH 4 Cl (20 ml) is slowly added to the reaction mixture and the phases are separated. The aqueous phase is washed with diethyl ether (3 x 10 ml) and the combined organic extracts are dried over MgSO 4 filtered and evaporated to yield a yellow oil which on distillation gives the title compound as a clear colourless liquid (0.91g, 85%, b. p. 150°C / 0.4 mm Hg). 31 P-NMR (CDCl 3 ): δ= 31.6. 1 H-NMR (CDCl 3 ): δ= 1.17 (t, 6H, 3 J HH = 7.2, OCH 2 CH 3 ), General Synthesis of Pyridylethenes 7a-e LDA (11 mmol) is formed as before and cooled to -78°C. α-Picoline (0.49 g, 5.2 mmol) in THF (10 ml) is added dropwise maintaining the temperature at -78°C, the solution is stirred for 30 min at this temperature. Diethyl chlorophosphate (0.98 g, 5.5 mmol) in THF (10 ml) is then added dropwise at -78°C, the reaction is stirred for 10 min and then allowed to warm to room temperature. A solution of aldehyde (5 mmol) in THF (10 ml) is then added at this temperature. After 10 min. stirring, H 2 O (20 ml) is added. The phases were separated and the aqueous phase extracted with diethyl ether (3 x 20 ml). The combined organic phases were dried over MgSO 4 and the solvents evaporated. Distillation by bulb-to-bulb apparatus yields the pyridylethenes 7a-e (see yields in Table 1). ). m/z (EI) 187 (M + , 28 %), 186 (100), 154 (10).

6H, 3 J

 3 PH = 16.6, PCMe 2 ), 3.8-4.02 (m, 4H, OCH 2 ), 7.05-7.15 (m, 1H, pyr-H), 7.58-7.64 (m, 2H, pyr-H), 8.52-8.57 (m, 1H, pyr-H). 13 C-NMR (CDCl 3 ) (JMOD): δ= 15.6 (OCH 2 CH 3 ), 20.27 (pyr-CH 3 ), 62.11 (OCH 2 ), 126.26 (pyr-CH), 128.39 (d, J=25, pyr-CH), 146.85 (d, J=13, pyr-C), 149.6 (d, J=19, pyr-CH), 151.15 (d, J=177, pyr-C).

( 2 '

 2 -pyridyl)-2-(3',4'-dioxymethylenephenyl)ethene (E) 7db.p. 235°C / 0.1 mmHg.[START_REF] Engel | Synthesis of Carbon-Phosphorus bonds[END_REF] H-NMR (CDCl 3 ): δ= 5.99 (s, 2H, OCH 2 O), 6.81 (d, 2H, 3 J HH = 8.0, Ar), 7.01 (d, 1H, 3 J HH = 16.1, H 2 ), 7.03 (dd, 1H, 3 J HH = 8.0, 4 J HH = 1.7, Ar), 7.13 (d, 1H, 4 J HH = 1.7, Ar), 7.14 (ddd, 1H, 3 J HH = 7.4, 4.9, 4 J HH = 1.0, pyr-H 5' ), 7.36 (d, 1H, 3 J HH = 8.0, pyr-H 3' ), 7.57 (d, 1H, 3 J HH = 16.1, H 1 ), 7.67 (td, 1H,3 J HH = 7.7, 4 J HH = 1.8, pyr-H 4' ), 8.58 (dm, 1H, 3 J HH = 4.9, pyr-H 6' ).[START_REF] Windsceif | [END_REF] C-NMR (CDCl 3 ): δ= 101.8 (s, OCH 2 O), 106.4 (s, Ar), 109.0 (s, C 2' Ar), 122.4 (s, pyr-C 3' ), 122.5 (s, pyr-C 5' ), 123.1 (s, Ar), 126.4 (s, C 2 ), 131.6 (s, Ar), 133.1 (s, C 1 ), 137.2 (s, pyr-C 4' ), 148.5 (s, Ar), 148.7 (s, Ar), 149.9 (s, pyr-C 6' ), 156.1 (s, pyr-C 2' ). NMR (CDCl 3 ): δ= 6.94 (d, 1H, 3 J HH = 15.8, H 2 ), 6.91 (dd, 1H, 3 J HH = 5.1, 3.6, thioph-H 4' ), 7.11 (ddm, 1H,3 J HH = 7.5, 4.8, pyr-H 5' ), 7.14 (d, 1H, 3 J HH = 3.6, thioph-H 3' ), 7.23 (d, 1H, 3 J HH = 5.1, thioph-H 5' ), 7.29 (d, 1H, 3 J HH = 7.5, pyr-H 3' ), 7.62 (td, 1H, 3 J HH = 7.5,4

  CDCl 3 ): δ= 24.9. 1 H-NMR (CDCl 3 ): δ= 1.20-1.27 (m, 6H, OCH 2 CH 3 ), 3.38 (d, 2H, 2 J PH = 22, PCH 2 ), 3.98-4.12 (m, 4H, OCH 2 ), 7.15-7.18 (m, 1H, pyr-H), 7.34-7.38 (m,

	31 P-NMR (1H, pyr-H), 7.58-7.62 (m, 1H, pyr-H), 8.50-8.53 (m, 1H, pyr-H).	13 C-NMR
	(CDCl 3 )(JMOD): δ= 15.72 (OCH 2 CH 3 ), 34.61 (d, 1 J PC = 135, PCH 2 ), 61.42 (OCH 2 ), 121.30
	(pyr-CH), 123.69 (pyr-CH), 135.75 (pyr-CH), 148.77 (pyr-CH), 152.18 (pyr-C).
	N
	H 3 C
	H 3 CO
	O
	O
	S

7 g, 85 %, bp 240-250°C / 20 mmHg).

  After 30 min, D 2 O (30 ml) is added at room temperature with vigorous stirring of the mixture. When deuteriolysis is complete (60 min) the reaction mixture is treated with an aqueous solution of NH 4 Cl. The aqueous phase is extracted with CH 2 Cl 2 (3 x 50 ml), the organic layer is dried (MgSO 4 ) and the solvents are evaporated off under reduced PC = 135, PCD 2 ), 61.76 (OCH 2 ), 121.56 (pyr-CH), 123.82 (pyr-CH), 136.22 (pyr-CH), 148.81 (pyr-CH), 152.10 (pyr-C).

	pressure. The crude product is then purified by bulb-to-bulb distillation to yield compound 5
	as a yellow oil. (9.6g, 83%, b.p. 230-240°C / 20 mmHg).

31 P-NMR (CDCl 3 ): δ= 25.4. 1 H-NMR (CDCl 3 ): δ= 1.20-1.27 (m, 6H, OCH 2 CH 3 ), 4.0-4.18 (m, 4H, OCH 2 ), 7.15-7.21 (m, 1H, pyr-H), 7.36-7.41 (m, 1H, pyr-H), 7.6-7.68 (m, 1H, pyr-H), 8.52-8.57 (m, 1H, pyr-H). 13 C-NMR (CDCl 3 ) (JMOD): δ= 15.79 (OCH 2 CH 3 ), 35.8 (d, 1 J

  CDCl 3 ): δ= 7.17 (ddd, 2H,3 J HH = 7.5, 4.8, 4 J HH = 1.2, pyr-H 5' ), 7.41 (dt, 2H,3 J HH = 7.8, 4 J HH = 1.0, pyr-H 3' ), 7.67 (td, 2H, 3 J HH = 7.7, 4 J HH = 1.8, pyr-H 4' ), 7.69 (s, 2H, H 1 , H 2 ), 8.62 (ddd, 2H, 3 J HH = 4.8, 4 J HH = 1.8, 0.9, pyr-H 6' ).[START_REF] Windsceif | [END_REF] C-NMR (CDCl 3 ): δ= 123.3 (s, pyr-C 3' ), 124.1 (s, pyr-C 5' ), 132.3 (s, C 1 , C 2 ), 137.4 (s, pyr-C 4' ), 150.4 (s, pyr-C 6' ), 155.6 (s, pyr-C 2' ). m/z (EI) 182 (M + , 31 %), 181 (100), 179 (2), 154(8).

	(2'-pyridyl)-2-(4'-methylphenyl)ethene (E) 7b
	b.p. 195-200°C / 0.2 mmHg.
	1 H-NMR (CDCl 3 ): δ= 2.37 (s, 3H, CH 3 ), 7.11 (ddd, 1H, 3 J HH = 7.5, 2.6, 4 J HH = 1.2, pyr-H 5' ),
	7.13 (d, 1H, 3 J HH = 16.0, H 2 ), 7.18 (d, 2H, 3 J HH = 8.8, Ar-H), 7.37 (dt, 1H, 3 J HH = 7.9, 4 J HH =
	0.9, pyr-H 3' ), 7.48 (d, 2H, 3 J HH = 8.8, Ar-H), 7.61 (d, 1H, 3 J HH = 16.0, H 1 ), 7.65 (td, 1H,
	3 J HH = 7.7, 4 J HH = 1.8, pyr-H 4' ), 8.60 (ddd, 1H, 3 J HH = 4.8, 4 J HH = 1.7, 0.8, pyr-H 6' ). 13 C-
	NMR (CDCl 3 ): δ= 21.8 (s, CH 3 ), 122.4 (s broad , pyr-C 3' , C 5' ), 127.1 (s, C 2 ), 127.6 (s, Ar),
	129.9 (s, Ar), 132.3 (s, C 1 ), 134.2 (s, Ar), 137.2 (s, pyr-C 4' ), 138.9 (s, Ar), 149.8 (s, pyr-C 6' ),
	156.1 (s, pyr-C 2' ). m/z (EI) 195 (M + , 25 %), 194 (100), 179 (4).
	(2'-pyridyl)-2-(4'-methoxyphenyl)ethene (E) 7c
	b.p. 245-250°C / 1 mmHg.
	1 H-NMR (CDCl 3 ): δ= 3.83 (s, 3H, OCH 3 ), 6.91 (d, 2H, 3 J HH = 8.8, Ar-H), 7.05 (d, 1H,
	(s, C 2 ), 129.1 (s, Ar), 130.0 (s, Ar), 133.1 (s, C 1 ), 137.3 (s, pyr-C 4' ), 150.0 (s, pyr-C 6' ), 156.4
	(s, pyr-C 2' ), 160.5 (s, Ar). m/z (EI) 211 (M + , 37 %), 210 (100), 195 (14), 167 (34).
	1,2-bis(2'-pyridyl)ethene (E) 7a
	b.p. 210-215°C / 0.2 mmHg.

[START_REF] Engel | Synthesis of Carbon-Phosphorus bonds[END_REF]

H-NMR (3 J HH = 16.0, H 2 ), 7.14 (ddd, 1H, 3 J HH = 7.4, 2.6, 4 J HH = 1.1, pyr-H 5' ), 7.37 (dt, 1H, 3 J HH = 7.9, 4 J HH = 0.9, pyr-H 3' ), 7.53 (d, 2H,

3

J HH = 8.8, Ar-H), 7.59 (d, 1H, 3 J HH = 16.0, H 1 ), 7.67 (td, 1H,

3

J HH = 7.7, 4 J HH = 1.8, pyr-H 4' ), 8.59 (ddd, 1H, 3 J HH = 4.8, 4 J HH = 1.7, 0.8, pyr-H 6' ).

[START_REF] Windsceif | [END_REF]

C-NMR (CDCl 3 ): δ= 55.9 (s, OCH 3 ), 114.8 (s, Ar), 122.3 (s, pyr-C 3' ), 122.4 (s, pyr-C 5' ), 126.2

  J HH = 1.8, pyr-H 4' ), 7.76 (d, 1H,3 J HH = 15.8, H 1 ), 8.56 (dm, 1H, 3 J HH = 4.8, pyr-H 6' ).[START_REF] Windsceif | [END_REF] C-NMR (CDCl 3 ): δ= 122.6 (s, pyr-C 3' ), 122.8 (s, pyr-C 5' ), 126.1 (s, thioph-C 3' ), 126.3 (s, thioph-C 4' ), 127.8 (s, thioph-C 5' ), 128.4 (s, C 1 , C 2 ), 137.2 (s, pyr-C 4' ), 142.8 (s, thioph-C 2' ), 150.3 (s, pyr-C 6' ), 155.8 (s, pyr-C2'

Acknowledgements

We thank the CNRS for grants for J. C. and R. W.. We are also grateful to Mr. M. Levard (URA CNRS 1307) of the Ecole Polytechnique for the mass spectra.


Simon Gsell
email: simon.gsell@univ-amu.fr

Umberto D'ortona

Julien Favier

Lattice-Boltzmann simulation of creeping generalized Newtonian flows: theory and guidelines

Keywords: lattice-Boltzmann, creeping flows, non-Newtonian flows

The accuracy of the lattice-Boltzmann (LB) method is related to the relaxation time controlling the flow viscosity. In particular, it is often recommended to avoid large fluid viscosities in order to satisfy the low-Knudsen-number assumption that is essential to recover hydrodynamic behavior at the macroscopic scale, which may in principle limit the possibility of simulating creeping flows and non-Newtonian flows involving important viscosity variations. Here it is shown, based on the continuous Boltzmann equations, that a two-relaxation-time (TRT) model can however recover the steady Navier-Stokes equations without any restriction on the fluid viscosity, provided that the Knudsen number is redefined as a function of both relaxation times. This effective Knudsen number is closely related to the previously-described parameter controlling numerical errors of the TRT model, providing a consistent theory at both the discrete and continuous levels. To simulate incompressible flows, the viscous incompressibility condition M a 2 /Re 1 also needs to be satisfied, where M a and Re are the Mach and Reynolds numbers. This concept is extended by defining a local incompressibility factor, allowing one to locally control the accuracy of the simulation for flows involving varying viscosities. These theoretical arguments are illustrated based on numerical simulations of the two-dimensional flow past a square cylinder. In the case of a Newtonian flow, the viscosity independence is confirmed for relaxation times up to 10 4 , and the ratio M a 2 /Re = 0.1 is small enough to ensure reliable incompressible simulations. The Herschel-Bulkley model is employed to introduce shear-dependent viscosities in the flow. The proposed numerical strategy allows to achieve major viscosity variations, avoiding the implementation of artificial viscosity cut-off in high-viscosity regions. Highly non-linear flows are simulated over ranges of the Bingham number Bn ∈ [1, 1000] and flow index n ∈ [0.2, 1.8], and successfully compared to prior numerical works based on Navier-Stokes solvers. This work provides a general framework to simulate complex creeping flows, as encountered in many biological and industrial systems, using the lattice-Boltzmann method.

Introduction

The lattice-Boltzmann (LB) method is a popular method to simulate fluid flows [START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF][START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. In contrast to Navier-Stokes methods which directly describe the dynamics of the macroscopic flow quantities (e.g. pressure, momentum), the LB method is a statistical approach based on gas kinetic theory. Originally derived from lattice-gaz models [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF][START_REF] Higuera | Lattice gas dynamics with enhanced collisions[END_REF], the LB method is a space-time discretization of the Boltzmann equation, which can recover the Navier-Stokes equations provided that the relevant conditions are satisfied [START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF][START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF].

In many textbooks, the ability of the LB method to accurately predict macroscopic flow behaviors is assumed to rely on a small Knudsen number assumption, especially because this number is used as the expansion parameter in the Chapman-Enskog analysis, allowing to derive the Navier-Stokes equations from the Boltzmann equation [START_REF] Chapman | The Mathematical Theory of Non-uniform Gases[END_REF]. The Knudsen number, denoted by Kn, is the ratio between the mean free path of gas particles and the typical length scale of the macroscopic problem, denoted by D in the following. It is often expressed using the von Kármán formulation, namely Kn ∼ M a/Re [START_REF] Karniadakis | Microflows and nanoflows: fundamentals and simulation[END_REF], where M a and Re are the Mach and Reynolds numbers of the physical system. The Mach number reads M a = u 0 /c s , where u 0 is the typical flow velocity of the problem and c s is the sound speed, which is set by the lattice in most LB implementations, unless a more generic method is used [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF]. As the condition M a 1 must be satisfied to achieve incompressible flows, the condition Kn 1 is simple to achieve in problems involving intermediate or large values of Re. In contrast, when Re 1, the Mach number may have to be decreased, by decreasing the flow velocity u 0 , in order to satisfy the low Knudsen assumption, resulting in higher computational costs. Due to this restriction, the lattice-Boltzmann method is not efficient in solving low Reynolds number (creeping) flows.

Alternatively, LB simulations of creeping flows can be performed by assuming the Stokes regime, where the flow is fully linear [START_REF] Talon | Assessment of the two relaxation time Lattice-Boltzmann scheme to simulate Stokes flow in porous media[END_REF]. In this regime, the Reynolds number is not a physical parameter and it can thus be freely varied to satisfy the low Knudsen assumption. However, this approach requires an a priori assumption on the effect of flow inertia, which may not be obvious depending on the physical configuration [START_REF] Lauga | Continuous breakdown of Purcells scallop theorem with inertia[END_REF]. In this work, a direct simulation approach is considered, where no assumption is made concerning the flow behavior, allowing one to perform simulations over a large range of Re including non-inertial and weakly-inertial flows.

In LB simulations, the fluid viscosity is controlled through a collision relaxation time, denoted by τ . Based on the kinematic viscosity ν = c 2 s (τ -1/2), the Reynolds number can be expressed as

Re = u 0 D c 2 s (τ -1/2) . (1)
Using the von Kármán formulation, the Knudsen number can thus be expressed as a function of the relaxation time, the sound speed and the macroscopic length scale,

Kn = c s (τ -1/2) D . (2)
It appears that the Knudsen restriction is closely connected to the relaxation time τ : to achieve small Knudsen numbers, small values of τ and thus small fluid viscosities have to be employed, if c s and D are fixed. Note that in some LB implementations, c s can be considered as a free parameter that is not connected to the fluid viscosity [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF][START_REF] Ginzburg | A free-surface lattice Boltzmann method for modelling the filling of expanding cavities by Bingham fluids[END_REF][START_REF] Talon | Assessment of the two relaxation time Lattice-Boltzmann scheme to simulate Stokes flow in porous media[END_REF][START_REF] Khirevich | Coarse-and fine-grid numerical behavior of MRT/TRT lattice-Boltzmann schemes in regular and random sphere packings[END_REF]. Even though such an approach may allow small variations of Kn through c s , the present work focuses on the dominant effect of τ , and c s is kept constant as in common LB implementations [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF][START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF]. The analysis and correction of the effect of the fluid viscosity (or relaxation time) on the accuracy and stability of the lattice-Boltzmann method have motivated a number of works during last decades. Several approaches have been proposed to achieve viscosity independence of the simulations, especially by introducing multiple-relaxation-time collision models [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability[END_REF] with suitable scaling of the relaxation parameters [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions[END_REF][START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF][START_REF] Khirevich | Coarse-and fine-grid numerical behavior of MRT/TRT lattice-Boltzmann schemes in regular and random sphere packings[END_REF]. The fluid viscosity may also greatly alter the accuracy of boundary conditions and immersed boundaries, which may require specific corrections [START_REF] Gsell | Explicit and viscosity-independent immersed-boundary scheme for the lattice Boltzmann method[END_REF]. In the two-relaxation-time (TRT) collision model [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions[END_REF][START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF], symmetric and anti-symmetric parts of the distribution functions are relaxed at two different rates. While the first rate controls the kinematic viscosity, the second rate is not directly related to macroscopic properties and it can be freely varied to adapt the numerical accuracy. D'Humières and Ginzburg [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF] showed that the numerical error can then become independent of the fluid viscosity if both relaxation rates satisfy a constant relation [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF]. This result should allow one to employ large fluid viscosities to simulate low Reynolds number flows efficiently (see equation [START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF]). However, the viscous relaxation time must in principle remain limited by the above mentioned Knudsen restriction (see equation [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]). The present work aims at clarifying this aspect.

In addition to its impact on low Re simulations, the limitation of the relaxation time also questions the reliability of the LB method for the simulation of non-Newtonian flows. Viscoplastic fluids, as shear-thinning, shear-thickenning and yield-stress (or Bingham) fluids, are often modeled using shear-dependent viscosities, an approach referred to as generalized Newtonian modelling [START_REF] Irgens | Rheology and Non-Newtonian Fluids[END_REF]. The LB method is a promising tool for the simulation of such flows, as it allows local computation of the flow shear rate using the distribution functions [START_REF] Krüger | Shear stress in lattice Boltzmann simulations[END_REF] and local variation of the fluid viscosity through the relaxation time [START_REF] Gabbanelli | Lattice Boltzmann method for non-Newtonian (power-law) fluids[END_REF][START_REF] Boyd | A second-order accurate lattice Boltzmann non-Newtonian flow model[END_REF][START_REF] Chai | Multiple-relaxation-time lattice Boltzmann model for generalized Newtonian fluid flows[END_REF][START_REF] Conrad | Accuracy of non-Newtonian lattice Boltzmann simulations[END_REF]. However, the resulting variation of the fluid viscosity raises the question of the relevant range of viscosities that can be covered while ensuring numerical accuracy. In particular, Bingham and shear-thinning fluids involve a major viscosity increase in low-shear regions of the flows. Even though the TRT lattice-Boltzmann model has already been successfully applied to Bingham flows in previous works [START_REF] Ginzburg | A free-surface lattice Boltzmann method for modelling the filling of expanding cavities by Bingham fluids[END_REF][START_REF] Talon | On the determination of a generalized darcy equation for yield-stress fluid in porous media using a Lattice-Boltzmann TRT scheme[END_REF], general guidelines to design such simulations using LBM are still missing.

In this work, a detailed analysis of the TRT lattice-Boltzmann equation is presented in order to emphasize the impact of the physical and numerical parameters on the macroscopic flow behavior in the creeping regime, with a special attention paid to the effect of the relaxation times. A focus is placed on the steady case, as unsteady diffusive phenomena can not be accurately predicted at large viscosities using the LB method (see §2.4). The analysis shows that the Navier-Stokes equations can be recovered without any restriction on the fluid viscosity, Mach number, Reynolds number or nominal Knudsen number Kn = M a/Re, provided that the physical Knudsen number is correctly redefined as a function of both relaxation times. The incompressible flow equations, however, can only be recovered if the viscous incompressibility factor T = M a 2 /Re is small. This result is verified numerically based on the two-dimensional creeping flow past a square cylinder. This configuration is chosen as it involves non-trivial external flow dynamics while allowing accurate and simple implementation of the no-slip boundary conditions using the bounce-back method. In the case of a viscoplastic fluid modelled by the Herschel-Bulkley law, a local definition of the T criterion allows one to control the viscosity range of the flow solution and to achieve considerable viscosity ratios between low-shear and high-shear regions of the flow. The numerical results are discussed and compared to prior works based Navier-Stokes models.

The present paper is organized as follows. First, the lattice-Boltzmann method is briefly introduced in §2. Mesoscopic equations are then analyzed in §3, and the macroscopic flow equations are derived and discussed. The method is applied to the simulation of a Herschel-Bulkley flow past a square cylinder. The physical setup and the numerical implementation are first described in §4; the numerical results are then analyzed in §5. Finally, the principal conclusions of this work are summarized in a last section, in §6.

Lattice-Boltzmann method

The main aspects of the lattice-Boltzmann model employed in this work are presented in the following. After a brief presentation of the discrete-velocity Boltzmann equation in §2.1, the two-relaxation-time collision model is described in §2.2 and the space-time discretization leading to the LB equation is introduced in §2.3. The ability of the present method to simulate unsteady flows is discussed in §2.4.

Boltzmann equation and velocity discretization

The flow dynamics is described through the particle distribution function f (x, ξ, t), representing the density of fluid particles moving with velocity ξ, at location x and time t. The dynamics of the distribution function is governed by the Boltzmann equation,

∂f ∂t + ξ • ∇f = Ω(f ), (3)
where Ω(f ) is the collision operator. The discretization of the Boltzmann equation in the velocity space, physical space, and time leads to the lattice-Boltzmann method.

The velocity space is discretized on a set of velocity vectors {c i , i = 0, ..., Q -1}, where Q is the number of discrete velocities. In the present work, which only focuses on twodimensional physical configurations, a D2Q9 velocity set is used. However, it should be mentioned that the following analysis does not rely on a particular velocity discretization. In the D2Q9 model, the velocity space is discretized on nine velocities, namely

c i =              (0, 0), i = 0, c cos( π(i-1) 2 ), sin( π(i-1) 2 ) , i ∈ [1, 4], √ 2c cos( π(2i-9)

4

), sin( π(2i-9)

4

)

, i ∈ [5, 8], ( 4
)
where c is the lattice speed. The particle densities at velocities {c i } are represented by the discrete-velocity distribution functions {f i (x, t)}, also called particle populations.

The macroscopic flow quantities are moments of the particle populations in the velocity space. In particular, the fluid momentum and density read

ρu = 8 i=0 f i c i , ρ = 8 i=0 f i , (5)
with ρ and u the fluid density and velocity.

Two-relaxation-time Boltzmann equation

In the discretized velocity space ξ = {c i }, the two-relaxation-time Boltzmann equation reads

∂f i ∂t + ∂ ∂x α (c i,α f i ) = - 1 τ + f + i -f eq+ i - 1 τ -f - i -f eq- i , (6)
where τ + and τ -are the two relaxation times and f + i and f - i are symmetric and antisymmetric parts of f i , expressed as

f + i = 1 2 (f i + f i ) ; f - i = 1 2 (f i -f i ), (7)
with the index i defined so that c i = -c i . When τ + = τ -, equation ( 6) reduces to the singlerelaxation-time (SRT) Boltzmann equation based on the Bhatnagar-Gross-Krook (BGK) collision operator [START_REF] Bhatnagar | A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems[END_REF].

Populations are decomposed into equilibrium and non-equilibrium parts, namely f i = f eq + f neq . The equilibrium functions are expressed as

f (eq) i = w i ρ 1 + u • c i c 2 s + (u • c i ) 2 2c 4 s - u • u 2c 2 s , (8)
where c s is the sound speed and {w i } are the lattice weights depending on the lattice model. In the case of the D2Q9 scheme, they read w 0 = 4/9, w 1 = w 2 = w 3 = w 4 = 1/9 and w 5 = w 6 = w 7 = w 8 = 1/36. The symmetric and anti-symmetric populations exhibit a series of properties that are described in §Appendix A. These identities will be particularly useful when projecting the Boltzmann equation onto Hermite polynomials in §3.1.

Space-time discretization

Space-time discretization of the LB equation ( 6) is performed through a method of characteristics. It can be shown that the second-order discretized equation can be expressed in an explicit form, called the lattice-Boltzmann equation, namely [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] fi (x + c i ∆t, t + ∆t) -fi (x, t) = -

1 τ + f neq+ i - 1 τ - f neq- i , (9)
through the change of variable

fi = f i + 1 2 1 τ + f neq+ i + 1 τ -f neq- i , (10)
and defining f + i = ( fi + fi )/2 and fi = ( fifi )/2, and the numerical relaxation times τ + = τ + + 1/2 and τ -= τ -+ 1/2, with τ ± > 0 and τ ± > 1/2. Equation ( 9) is defined on a two-dimensional uniform and Cartesian grid which consists of a series of nodes spaced by a distance ∆n in the x and y directions. This lattice satisfies the acoustic scaling condition c = ∆n/∆t, where ∆t is the numerical time step. Therefore, the streaming process, on the left-hand side of equation [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF], is characterized by the transport of the populations from one node to a neighboring one at the lattice speed c. It is balanced by the collision process, on the right-hand side of the equation. D'Humières and Ginzburg [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF] have shown that the discretization errors on steady-state bulk solutions and specific boundary rules are controlled by an effective relaxation time, also called magic collision parameter, Λ = (τ + -1/2)(τ --1/2). Small discretization errors are generally expected for values of Λ close to 0.1 [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. In particular, the value Λ = 1/4 ensures optimal stability conditions [START_REF] Ginzburg | Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium[END_REF][START_REF] Kuzmin | The role of the kinetic parameter in the stability of tworelaxation-time advection-diffusion lattice boltzmann schemes[END_REF].

From unsteady to steady Boltzmann equations

The TRT Boltzmann equation ( 6), or its discretized counterpart [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF], are not suitable to describe highly viscous unsteady flows. As shown in the following, the fluid viscosity is determined by the relaxation time τ + , namely ν = c 2 s τ + , where ν is the kinematic viscosity. A typical numerical time scale related to viscous diffusion can thus be expressed as t ν = ∆n 2 /ν = 3/τ + , using the lattice scaling ∆n = ∆t = 1 and c s = 1/ √ 3. In order to accurately describe unsteady diffusive phenomena, the diffusive time scale should be larger than numerical time step, namely t ν > 1. In theory, this condition is only satisfied when τ + < 3. Even though in practice this critical value may vary from one case to the other, this condition illustrates that the lattice-Boltzmann method is generally not reliable for the simulation of highly viscous unsteady flows involving large values of the relaxation time.

Other time integration strategies, as dual-time-stepping approaches [START_REF] Gsell | Multigrid dual-time-stepping lattice Boltzmann method[END_REF], can be employed to simulate highly viscous unsteady flows using the lattice-Boltzmann method. This method relies on the accuracy of the lattice-Boltzmann in simulating steady fluid flows. Therefore, in the following the focus is first placed on the analysis of the steady TRT-LB method. The simulation of unsteady flows will be addressed in a future work.

Here, the unsteady LB equation is thus only employed as a method to reach steady flow solutions. The transient part of the simulations is ignored, and a focus is placed on the reliability of the obtained steady solutions. Once a steady solution is reached, the system satisfies the steady TRT lattice-Boltzmann equation. Its continuous counterpart is the steady Boltzmann equation, which is obtained by ignoring the time derivative in equation [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF],

∂ ∂x α (c i,α f i ) = - 1 τ + f + i -f eq+ i - 1 τ -f - i -f eq- i . (11)

Analysis of the steady TRT Boltzmann equation

In the following, the steady TRT Boltzmann equation is analyzed in order to characterize the emerging macroscopic flow behavior, especially in the viscous flow regime. The analysis is based on the velocity-discretized Boltzmann equation [START_REF] Lauga | Continuous breakdown of Purcells scallop theorem with inertia[END_REF], in order to ignore numerical errors related to the space-time discretization, which have already been thoroughly examined in prior works [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF].

Hermite projection of the Boltzmann equation

To recover the macroscopic equations, equation ( 11) is projected onto Hermite polynomials of increasing order, through the discrete operator

< f i |H (n) i >= i f i H (n) i , (12)
where H (n) is the n-rank Hermite tensor. Hermite tensors from order n = 0 to order n = 2 can be expressed as [START_REF] Coreixas | Recursive regularization step for highorder lattice Boltzmann methods[END_REF] H

(0) i = 1, H (1)
i,α = c i,α , H (2)
i,αβ = c i,α c i,β -c 2 s δ α,β . (13)
It should be noted that the following analysis is based on the equilibrium distribution functions defined by equation [START_REF] Karniadakis | Microflows and nanoflows: fundamentals and simulation[END_REF], that are already truncated in the Hermite space. Indeed, this analysis only aims at deriving the macroscopic equations recovered by the present LB model, and the effect of Hermite truncation is not directly addressed. At order 0, the Boltzmann equation [START_REF] Lauga | Continuous breakdown of Purcells scallop theorem with inertia[END_REF] recovers the mass equation,

∂ ∂x α (ρu α ) = 0. ( 14
)
Then the first-order projection,

∂ ∂x α ( i f i c i,α c i,β ) = 0, (15)
involves the second-order moment Π αβ = i f i c i,α c i,β . Its equilibrium part Π eq αβ = i f eq i c i,α c i,β can be directly computed from the expression of f eq i ; it reads Π eq αβ = ρc 2 s δ αβ + ρu α u β . Expression (15) thus becomes

∂ ∂x α (ρc 2 s δ αβ + ρu α u β ) + ∂ ∂x α ( i f neq i c i,α c i,β ) = 0. ( 16
)
Expression [START_REF] Gsell | Explicit and viscosity-independent immersed-boundary scheme for the lattice Boltzmann method[END_REF] governs the flow momentum transport. The first term of the equation relates to inviscid fluxes. The flow pressure can be identified by introducing the isothermal ideal gaz state equation, p = ρc 2 s . However, the non-equilibrium second-order moment, Π neq αβ = i f neq i c i,α c i,β , remains to be determined. In order to recover a hydrodynamic behavior, Π neq αβ should be equal to the viscous stress tensor.

The tensor Π neq αβ can be explicitly computed by performing a second-order projection of the Boltzmann equation, as detailed in Appendix B. The resulting equation reads

2ρc 2 s S βγ + ∂ ∂x α (Q neq αβγ ) -u γ ∂ ∂x α (Π neq αβ ) -u β ∂ ∂x α (Π neq αγ ) = -Π neq βγ τ + , (17)
where Q αβγ is the third-order moment defined by

Q αβγ = i f i c α c β c γ and S αβ = 1 2 ∂uα ∂x β + ∂u β ∂xα
is the shear-rate tensor. The viscous stress tensor can be recognized in the first term of equation [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF]. However, other terms represent non-hydrodynamic contributions to the macroscopic flow behavior. These terms are quantified hereafter following the Chapman-Enskog procedure.

Chapman-Enskog expansion

The distribution functions are expanded around the equilibrium solution f i = f eq i using the expansion parameter 1, namely

f i = N n=0 n f (n) i , ( 18
)
where N is the order of the expansion. The definition of the expansion parameter will be discussed in the next section. At order 0, the distribution functions are equal to the equilibrium distributions, f

= f eq i . The multi-scale expansion consists in substituting this expansion in the Boltzmann equation and projecting the equation on different orders of .

The space derivative is considered to be first order, and expanded in the form

∂ ∂x = ∂ ∂x (1)
.

If is small enough, f i can be reasonably approximated by f i ≈ f eq i + f (1) . By expanding the shear-rate tensor as S βγ = S

βγ = -2τ + ρc 2 s S (1)
βγ -2 τ + ∂ (1) ∂x α (Q (1) αβγ ) + u γ τ + 2 ∂ (1) ∂x α (Π (1)
αβ ) + u β τ + 2 ∂ (1) ∂x α (Π (1) αγ ), (19)
which at first order reduces to Π

(1)

βγ = -2τ + ρc 2 s S (1) βγ , (20)
or, reversing the spatial derivative expansion,

Π neq βγ = -2τ + ρc 2 s S βγ . ( 21
)
Substituting Π neq βγ in the momentum equation ( 16) finally gives the expected steady momentum equation. By adding the mass equation ( 14), the Navier-Stokes equations are obtained,

∂ ∂x α (ρu α ) =0, (22a)
∂ ∂x β (ρu α u β ) = - ∂p ∂x α + ∂ ∂x β (2µS αβ ), (22b)
with p = ρc 2 s and µ = ρc 2 s τ + . The relaxation time τ + , which controls the fluid viscosity, is called the viscous relaxation time in the following. In contrast, τ -is not involved in the macroscopic equations; it is called the free relaxation time. Note that in the above analysis, no assumption has been made concerning the viscous relaxation time τ + ; in particular, τ + can be space dependent [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. The definition µ = ρc 2 s τ + thus applies to non-uniform viscosities, e.g. shear-dependent viscosities, as already done in prior works [START_REF] Ginzburg | A free-surface lattice Boltzmann method for modelling the filling of expanding cavities by Bingham fluids[END_REF][START_REF] Talon | On the determination of a generalized darcy equation for yield-stress fluid in porous media using a Lattice-Boltzmann TRT scheme[END_REF].

The above Chapman-Enskog analysis is equivalent for both the single-and two-relaxationtime collision models, since the multi-scale expansion only has to be performed on the second-order projection of the Boltzmann equation [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF], which only involves the viscous relaxation time. Therefore, both models may appear to be equivalent at the macroscopic level. In both cases, the derived hydrodynamic equations ( 22) rely on the small value of the expansion coefficient , as all terms ∼ O( 2) have been neglected. If 1 is not satisfied, additional macroscopic equations can emerge and lead to non-hydrodynamic behaviors, e.g. visco-elastic effects [START_REF] Giraud | A lattice Boltzmann model for Jeffreys viscoelastic fluid[END_REF]. The definition of is thus crucial to close the analysis. As discussed in the following, this aspect reveals important differences between the SRT and TRT models.

Discussion on the Knudsen number

The Chapman-Enskog expansion aims at expressing deviations from the equilibrium solution through small perturbations whose amplitude is quantified by the expansion parameter . Therefore, this parameter controls the magnitude of the non-equilibrium functions; in particular, f neq i = 0 must be satisfied when = 0. In the case of the SRT collision model, the solution tends to the equilibrium distribution when the relaxation time tends to 0. Therefore, τ + can be employed as a perturbative parameter, namely SRT ∼ τ + . Generally, the expansion parameter is rather expressed as a Knudsen number Kn, quantifying the ratio between the molecular free path length l and the typical hydrodynamic length scale D. The Knudsen number can be determined based on the Reynolds Re and Mach M a numbers through the von Kármán relation, namely Kn = M a/Re, which can be re-written as Kn = τ + c s /D. If the value of c s is fixed by the lattice as in the present LB model, it is seen that both definitions SRT ∼ τ + or SRT ∼ Kn are equivalent for a given value of D. The Chapman-Enskog expansion can only be performed if Kn 1; consequently, the SRT lattice-Boltzmann method only recovers the Navier-Stokes equations for low Knudsen number configurations. The relation ∼ τ + clearly indicates the limitation of the SRT collision model: the viscous relaxation time, i.e. the fluid viscosity, has an influence on the resulting macroscopic equations, since large relaxation times (or large Kn) will result in the growth of non-hydrodynamic modes in the momentum equation [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF].

The definition ∼ τ + can not be employed in the case of the TRT collision model, since both τ + and τ -are expected to impact the magnitude of the non-equilibrium functions. The relevant expansion parameter can be determined by expressing the non-equilibrium populations as functions of both relaxation times, as detailed in Appendix C. This analysis shows that, in the case of the TRT collision model, the equilibrium solution is achieved when the product of both relaxation times tends to 0, namely τ + τ -→ 0. Note that this result holds for continuous equations, as additional terms may emerge at the discrete level. The parameter T RT ∼ √ τ + τ -can thus be considered as a relevant parameter for the expansion. Interestingly, this number is equivalent to the magic collision number Λ = τ + τ -, identified by d'Humières and Ginzburg [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF], that determines the magnitude of the numerical errors of the TRT lattice-Boltzmann equation, i.e. the TRT Boltzmann equation ( 11) after spatial discretization. Here, it is proposed to express it as an effective Knudsen number,

Kn * = √ τ + τ -c s /D = √ Λc s /D
, corresponding to the relevant expansion parameter for the Chapman-Enskog analysis of both the SRT and TRT Boltzmann equations.

The Knudsen number controls the reliability of the Boltzmann equation in modeling fluid flows, as it determines the magnitude of high-order non-hydrodynamic terms in the expression of the shear-stress tensor (see equation [START_REF] Krüger | Shear stress in lattice Boltzmann simulations[END_REF]). Small values of Kn * should therefore be used in LB simulations. Using a SRT collision model, small Knudsen numbers can only be achieved by choosing small relaxation times, i.e. by performing weakly viscous simulations. In contrast, the TRT model allows one to freely increase the viscous relaxation time, provided that τ + τ -remains small. This property, already established based on the space-discretized equations [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF], is supported by the present continuous analysis, providing a consistent theory at both the discrete and continuous levels.

Non-dimensional macroscopic equations and incompressibility assumption

Generally, a hydrodynamic problem described by the Navier-Stokes equations ( 22) is expected to depend at least on five physical parameters, namely a reference density ρ 0 , a reference velocity u 0 , a reference kinematic viscosity (or relaxation time) ν 0 = c 2 s τ + 0 , the speed of sound c s and a reference length scale D. Following the Buckingham theorem, the non-dimensional flow solution should only depend on two non-dimensional parameters, that are often defined as the Reynolds number Re and Mach number M a,

Re = u 0 D ν 0 , M a = u 0 c s . ( 23
)
For inertial flows, a set of non-dimensional physical quantities can be defined as

u * = u/u 0 , ρ * = c 2 s (ρ -ρ 0 )/ρ 0 u 2 0 , ν * = ν/ν 0 and S * = SD/u 0 . The non-dimensional flow equations read ∇ • (u * ) + M a 2 ∇ • (ρ * u * ) = 0, (24a) ∇ • (u * u * ) + M a 2 ∇ • (ρ * u * u * ) = -∇(ρ * ) + 2 Re ∇ • (ν * S * ) + 2 Re M a 2 ∇ • (ρ * ν * S * ). (24b)
If the Mach number is small enough, these equations recover the steady incompressible Navier-Stokes equations,

∇ • (u * ) = 0, (25a) ∇ • (u * u * ) = -∇(ρ * ) + 2 Re ∇ • (ν * S * ). ( 25b
)
Note that ν * = 1 in the Newtonian case.

Alternatively, if the flow is expected to be dominated by viscous stresses, as in the present study, the non-dimensional density should be defined as ρ * = c 2 s (ρ -ρ 0 )D/µ 0 u 0 . A different set of non-dimensional equations is then obtained,

∇ • (u * ) + M a 2 Re ∇ • (ρ * u * ) = 0, (26a) Re∇ • (u * u * ) + M a 2 ∇ • (ρ * u * u * ) = -∇(ρ * ) + 2∇ • (ν * S * ) + 2 M a 2 Re ∇ • (ρ * ν * S * ). (26b)
In this case, the incompressible Navier-Stokes equations ( 25) are recovered if the Mach and Reynolds numbers satisfy M a 2 Re, or equivalently T = M a 2 /Re = τ + 0 u 0 /D 1, where T is the viscous incompressibility factor.

In summary, the above theoretical analysis indicates that highly viscous flow simulations can be performed using the proposed TRT-LB method, provided that the effective Knudsen number Kn * = τ + 0 τ - 0 c s /D = √ Λc s /D and the viscous incompressibility factor T = M a 2 /Re are small enough. These features are numerically investigated in the following on the basis of simulations of the flow past a square cylinder.

Physical and numerical setup

The theoretical analysis developed in §3 is applied in the following to the simulation of the flow past a square cylinder. Complex fluids are described using a Non-Newtonian Herschel-Bulkley model, described in §4.1. Details on the physical setup are presented in §4.2, and the numerical parameters are discussed in §4.3.

Non-Newtonian Herschel-Bulkley model

The Herschel-Bulkley model is employed to simulate non-Newtonian flows of viscoplastic fluids. The local dynamic fluid viscosity follows

µ = σ 0 γ + k γn-1 , (27)
where σ 0 is the yield stress, k is the flow consistency, n is the flow index, γ is the local shear-rate magnitude, expressed as [START_REF] Chai | Multiple-relaxation-time lattice Boltzmann model for generalized Newtonian fluid flows[END_REF][START_REF] Lashgari | First instability of the flow of shear-thinning and shear-thickening fluids past a circular cylinder[END_REF][START_REF] Mossaz | Non-recirculating and recirculating inertial flows of a viscoplastic fluid around a cylinder[END_REF][START_REF] Farnoush | An investigation on the body force modeling in a lattice Boltzmann BGK simulation of generalized Newtonian fluids[END_REF][START_REF] Wu | A modified lattice Boltzmann method for Herschel-Bulkley fluids[END_REF]

γ = √ 2S : S, ( 28
)
and S is the shear-rate tensor. When the fluid is Newtonian, σ 0 = 0 and n = 1, thus the fluid viscosity is equal to the fluid consistency. When n < 1, the fluid exhibits a shear-thinning behavior, i.e. its viscosity decreases as a function of the shear rate. The opposite trend, when n > 1, is called shear-thickening behavior. The yield stress introduces a rapid variation of the viscosity as the shear rate approaches σ 0 /µ 0 . This behavior aims at reproducing the shear-dependent solid /fluid transitions observed in many complex flows. The solid -like behavior is thus modeled by an infinitely large viscosity.

In practice, the local fluid viscosity is varied through the viscous relaxation time τ + ,

τ + = 1 2 + µ ρc 2 s . ( 29
)
The value of µ is determined by computing the local shear-rate tensor based on the nonequilibrium distribution functions, following equation ( 21),

S = - 1 2ρc 2 s τ + 8 i=1 f i -f (eq) i c i c i . ( 30
)
As τ + is involved in the expression of S, the value of τ + from the previous time step is used. This has no effect of the numerical results since only steady solutions are considered in this work.

In many numerical works, the Herschel-Bulkley law is truncated in order to avoid excessive viscosities in the computational domain that may deteriorate the numerical stability and accuracy. Regularized viscosity models have been proposed to achieve this truncation while avoiding sharp transitions between truncated and non-truncated regions [START_REF] Papanastasiou | Flows of materials with yield[END_REF]. The effect of this regularization has been analyzed in prior works [START_REF] Ginzburg | A free-surface lattice Boltzmann method for modelling the filling of expanding cavities by Bingham fluids[END_REF][START_REF] Mossaz | Criteria for the appearance of recirculating and non-stationary regimes behind a cylinder in a viscoplastic fluid[END_REF]. Here, an original Herschel-Bulkley law without truncation is employed, taking advantage of the viscosityindependence of the present numerical framework. However, the viscosity remains restricted by the floating point approximation and thus it can not reach infinite values. In particular, the shear-rate γ = 0 must be avoided in expression [START_REF] Kuzmin | The role of the kinetic parameter in the stability of tworelaxation-time advection-diffusion lattice boltzmann schemes[END_REF]. In practice, the maximum value of the viscosity is thus set by controlling the minimum shear-rate value in expression [START_REF] Kuzmin | The role of the kinetic parameter in the stability of tworelaxation-time advection-diffusion lattice boltzmann schemes[END_REF], by redefining γ as γ ← max( γ, γ ), with γ = 10 -14 .

Physical setup

The physical configuration is schematized in figure 1. A square cylinder of size D is immersed in an oncoming flow characterized by its velocity u 0 , density ρ 0 , consistency k, flow index n and yield stress σ 0 . A square computational domain is used; its size is denoted by L. The flow inlet is set by a velocity Dirichlet condition based on the bounce-back method. The density at the boundary is set to the reference velocity ρ 0 . At the outlet, a pressure Dirichlet condition is ensured through a non-equilibrium bounce-back method. The boundary density and tangential velocity are set to ρ 0 and 0, and the normal flow velocity is derived from density and momentum balance [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF]. Periodic boundary conditions are set on the upper and lower boundaries. Finally, the no-slip condition on the square cylinder is enforced by a halfway bounce-back method.

While the bounce-back method is expected to achieve viscosity-independent numerical errors at the boundary [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions[END_REF][START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF], ensuring an accurate solution at the cylinder surface, the viscosity independence is not theoretically ensured at the outlet where the non-equilibrium bounce-back approach is applied. However, this boundary has a smaller effect on the overall flow as it is placed far from the cylinder. In practice, a good viscosity-independence of the flow solution close to the cylinder is obtained in the present simulations, as illustrated in the following. Several non-dimensional numbers can be defined to characterize the physical setup. First, the Reynolds number reads

Re = ρ 0 u 0 D µ 0 , (31)
where µ 0 is the reference viscosity. It is associated with a reference viscous relaxation time τ + 0 . If the fluid is Newtonian or if it is a Bingham fluid (n = 1 and σ 0 ≥ 0), the reference viscosity is equal to the fluid consistency, µ 0 = k. When n = 1, the reference viscosity depends on the physical configuration. A reference shear rate can be defined as γ0 = u 0 /D, allowing one to derive the reference viscosity µ 0 = k γ0 n-1 . Overall, the general definition of the Reynolds number is [START_REF] Mossaz | Criteria for the appearance of recirculating and non-stationary regimes behind a cylinder in a viscoplastic fluid[END_REF]

Re = ρ 0 u 2-n 0 D n k . (32)
Then, the Bingham number,

Bn = σ 0 k D u 0 n , (33)
is used to quantify the ratio between the yield stress and the viscous stresses in the flow. The drag force exerted by the fluid on the body is used to examine the accuracy of the numerical simulations. The drag force is often normalized by the reference inertial stress to define the drag coefficient,

C x = 2F x ρ 0 u 2 0 D , (34)
where F x is the fluid force in the x direction. However, when the flow is governed by viscous stresses, it is more convenient to define another coefficient, referred to as viscous drag coefficient in the following,

C x,µ = C x Re = 2F x ku n 0 D 1-n . (35)
Finally, a third definition should be introduced when a creeping flow is dominated by the yield stress σ 0 at large Bingham numbers, namely

C x,σ 0 = C x Re Bn = 2F x σ 0 D . ( 36
)
In the present LB implementation, the fluid force F = (F x , F y ) exerted on the cylinder is determined through a momentum balance approach. Considering a surface Γ encompassing the cylinder, the fluid force is balanced, when the steady solution is reached, by the integrated stresses on Γ, namely

F + Γ σ • ndS = 0, (37)
where σ = -pI + 2µS is the stress tensor, I is the unit tensor and n is the unit normal vector on Γ. Note that the shear-rate tensor S is computed locally using expression [START_REF] Giraud | A lattice Boltzmann model for Jeffreys viscoelastic fluid[END_REF]. The force F is determined at each time step using equation [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF], after the update of the macroscopic quantities.

Numerical parameters

The theoretical analysis performed in §3 has shown that a number of numerical parameters has to be carefully set in order to ensure the reliability of the simulations. First, the two relaxation times τ + and τhave to be chosen so that the effective Knudsen number Kn * remains small. In practice, this is achieved by setting the relaxation time product Λ = (τ + -1/2)(τ --1/2) to a constant value, equal to 1/6 in the following. This value has been reported to be optimal for diffusive problems [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] and it keeps Kn * to a reasonably small value: Kn * = √ Λc s /D remains lower than 0.01 as long as the cylinder size D is larger than 10∆n. In non-Newtonian simulations, the free relaxation time τis thus varied together with τ + in order to maintain the constant value of Λ.

Then, the viscous incompressibility factor T = M a 2 /Re, with M a = u 0 /c s , is set to control the incompressibility of the flow. If the value of the Reynolds number is fixed, T determines the reference velocity u 0 which drastically impacts the computational cost: when the flow velocity is decreased, transient simulations get longer and more time steps have to be performed to reach a steady solution. The optimal value of T will be discussed in §5.

When the fluid is non-Newtonian, the local viscosity varies across the computational domain. This variation may be seen as a variation of the local Reynolds number, which in turn alters the incompressibility factor. The local Reynolds number is defined as the local ratio between inertial and viscous stresses, expressed as

Re l = ρ|u| 2 µ γ , ( 38
)
where |u| is the local velocity magnitude. A local Mach Number is also defined as M a l = |u|/c s . Consequently, the local incompressibility factor reads

T l = µ γ ρc 2 s . ( 39
)
According to the previous analysis, the fluid viscosity has no influence on the simulation accuracy as long as the incompressibility factor is small enough. Therefore, the local quantity T l is employed to control the variation of the viscosity and ensure the simulation reliability. Given a maximum value of T l that is acceptable to satisfy the incompressibility condition, called T max , the local maximum fluid viscosity is defined as µ max = T max ρc 2 s / γ. Following this procedure, it is seen that µ max tends to infinity when the shear rate tends to zero, which is expected to allow major viscosity variations in the simulations. In practice, when µ max is reached during a simulation, the Herschel-Bulkley law is artificially truncated, altering the accuracy of the non-Newtonian modelling. In this case, the reference incompressibility factor T can be further decreased until the whole flow solution satisfies the incompressibility condition, reproducing a realistic Herschel-Bulkley law without any cut-off and achieving major viscosity ratios across the computation domain.

Numerical results

The reliability of the present numerical method for the simulation of creeping non-Newtonian flows is examined in the following. First, a series of simulations is performed for intermediate values of the Reynolds number (Re ≥ 1) in order to set the grid resolution and illustrate the reliability of the method in this well-documented flow regime; these results are presented in §5.1. Then, the creeping flow regime (Re 1) is explored in §5.2, in particular to demonstrate the viscosity independence of the present numerical strategy. The method is then applied to viscoplastic fluids, namely yield-stress fluid in §5.3 and shear-thinning and shear-thickening fluids in §5.4.

Validation results at Re ≥ 1

Figure 2(a) shows the evolution of the drag force coefficient as a function of the mesh refinement, quantified by D/∆n. The Reynolds number is set to Re = 1 and the size of the computational domain is L = 50D. A clear converging trend is noted, although the overall variation of C x remains small over the considered range of D/∆n. Indeed, it should be noted that the relative amplitude of the C x axis in this plot is only about 5%. The influence of the domain size on the drag force coefficient at Re = 1 is examined in figure 2(b), for D/∆n = 5. The domain size significantly impacts C x , even when large domains are used, as expected due to the low value of the Reynolds number.

In the following, the domain size is limited in order to keep a reasonable computational cost and allow some parametric analyses. It should be mentioned that, as confinement effects are expected to increase as the Reynolds number decreases, the following simulations performed for Re 1 are necessarily expected to be altered by the domain size, unless a very large domain is used [START_REF] Pantokratoras | Steady flow of a power-law non-Newtonian fluid across an unconfined square cylinder[END_REF]. Convergence regarding this parameter is thus not required

This value of the drag force corresponds to simulations involving negligible blockage effects. It is seen in figure 2(c) that the present numerical setup accurately reproduces this drag evolution. In addition, as described in §5.3 and §5.4, the present numerical framework leads to accurate flow predictions in configurations involving viscoplastic fluids. The effect of the domain size, which varies as a function of fluid properties, is also further discussed in the following.

Creeping Newtonian flow past a square cylinder

The reliability of the present numerical framework in simulating highly viscous flows is examined in the following, in the case of a Newtonian fluid. Figure 3(a) shows the predicted drag force coefficient C x,µ = C x Re for two Reynolds numbers, Re = 10 -2 and Re = 10 -3 , over a range of the incompressibility factor T . In order to keep the Reynolds number constant, T is varied together with the viscous relaxation time τ + and Mach number M a. This procedure drastically impacts the computational efficiency, since more time steps have to be performed to achieve steady flow solutions when the flow velocity is reduced. According to the analysis performed in §3, T should be the only numerical parameter impacting the simulations, and numerical accuracy should be ensured as long as T 1. This is supported by the results in figure 3(a), where the drag coefficient exhibits a converging trend as the incompressibility factor is decreased. The constant value of C x,µ in the range T 1 indicates that other numerical parameters, in particular the viscous relaxation time, have no effect on the simulation result. Numerical simulations are considered to be reliable in this range, where numerical parameters approach their regular values, namely τ + ∼ 1, M a 1 and Kn 1 (see table 1 described hereafter). Some additional simulations have been performed using a domain length L = 100D to allow quantitative comparison with the numerical data reported by Pantokratoras [START_REF] Pantokratoras | Steady flow of a power-law non-Newtonian fluid across an unconfined square cylinder[END_REF] based on this particular geometry. In this configuration, for Re = 0.01, the drag predicted by the present method C x,µ = 6.74 is comparable to the value of 7.2 reported in that paper. Small differences between both simulations may however be due to the employed boundary conditions, since velocity Dirichlet conditions are set on the upper and lower boundaries of the domain in Pantokratoras' simulations, and thus the confinement conditions differ between both studies. Further comparison with the data issued from that study, including a discussion on the domain boundary effects, will be presented in §5. [START_REF] Higuera | Lattice gas dynamics with enhanced collisions[END_REF].

Values of C x,µ are very similar for both values of the Reynolds number, suggesting that the creeping flow regime Re ≈ 0 is reached. The relative variation of C x,µ when T is varied from 0.001 to 0.1 is lower than 1%. Therefore, T = 0.1 appears to be the optimal value in order to ensure numerical accuracy while minimizing the computational cost.

The effect of T on the flow in the vicinity of the cylinder is depicted in figures 3(b,c). Both the streamwise flow velocity and shear-rate magnitude remain almost unaltered for T ≤ 0.1. In this range, the flow patterns are consistent with the expected creeping flow features. In particular, the flow presents symmetries along both the streamwise and cross-flow directions. The highest shear regions are found close to the edges of the cylinder. In contrast, important flow variations occur for T > 0.1 due to the emergence of compressibility effects. A high shear region emerges downstream of the body and the flow symmetry is broken. Overall, the evolution of the flow pattern is consistent with the rapid variation of the drag coefficient observed in figure 3(a) for T > 0.1.

Details on the numerical setups associated with the data plotted in figure 3 are provided in table 1. It can be noted that the Mach number, Knudsen number and viscous relaxation time increase as functions of T . It should be recalled that Kn designates the von Kármán Knudsen number, namely Kn = M a/Re. The value of Kn varies over a wide range and it typically violates the condition Kn 1. Moreover, the viscous relaxation time greatly departs from unity. In particular, the case Re = 10 -3 and T = 10 -1 shows that accurate simulations can be performed up to Kn = 10 and τ + ≈ 360.

u x /u 0 -2 -1 0 1 2 -2 -1 0 1 2 y/D x/D -2 -1 0 1 2 -2 -1 0 1 2 1 x/D -2 -1 0 1 2 -2 -1 0 1 2 0 1 2 x/D -2 -1 0 1 2 -2 -1 0 1 2 3 -2 -1 0 1 2 -2 -1 0 1 2 (c) -2 -1 0 1 2 -2 -1 0 1 2 γD/u 0 -2 -1 0 1 2 -2 -1 0 1 2 y/D x/D -2 -1 0 1 2 -2 -1 0 1 2 1 x/D -2 -1 0 1 2 -2 -1 0 1 2 0 1 2 x/D -2 -1 0 1 2 -2 -1 0 1 2 3
To further illustrate that the simulation is not altered by the numerical parameters if the incompressibility factor is fixed, the flow past a square cylinder has been computed over a range of Reynolds numbers Re 1 for T = 0.1. The resulting evolution of the drag force coefficient is plotted in figure 4(a). As the Reynolds number decreases, the flow solution approaches the asymptotic creeping flow regime and the drag coefficient C x,µ rapidly converges to a constant value equal to 8.25, approximately. This limit value is close to the drag value of 8.4 obtained using a linear lattice-Boltzmann model, i.e. by keeping only the first two terms of the equilibrium functions in equation ( 8), to model the Stokes regime. The viscosity-independence of the drag force in the Stokes regime using the TRT model has been investigated in prior works [START_REF] Ginzburg | Multireflection boundary conditions for lattice Boltzmann models[END_REF][START_REF] Talon | Assessment of the two relaxation time Lattice-Boltzmann scheme to simulate Stokes flow in porous media[END_REF][START_REF] Khirevich | Coarse-and fine-grid numerical behavior of MRT/TRT lattice-Boltzmann schemes in regular and random sphere packings[END_REF]. Here, the constant evolution of the drag in the range Re 1 indicates that the non-linear simulations remain accurate even at very low Reynolds numbers. Other flow quantities also remain unaltered by variations of the Reynolds number in the region Re 1. This is depicted in figures 4(b,c), which show the evolutions of the streamwise flow velocity and shear rate along the cross-flow direction for Re = 10 -3 and Re = 10 -5 , in the vicinity of the body for x/D = 0. It can be noted that both simulations result in identical flow solutions. The numerical setups corresponding to the data plotted in figure 4 have also been reported in table 1. Knudsen numbers up to Kn ≈ 310 and relaxation times up to τ + ≈ 11000 can be achieved without altering the simulation accuracy, which further supports the theoretical analysis proposed in §3.

Creeping yield-stress flow past a square cylinder

According to the previous numerical analysis of the Newtonian flow past a square cylinder, the fluid viscosity can be freely varied in the present numerical framework as long as the condition T ≤ 0.1 is satisfied. In the following, the focus is placed on yield-stress flow simulations. In this case, large variations of the local incompressibility factor T l (39) are expected due to local variations of the flow viscosity. In particular, T l may locally exceed the nominal incompressibility factor T based on the reference viscosity. As described in §4.3, the numerical accuracy is thus ensured by setting the maximum value of T l in the simulations. Based on the previous Newtonian analysis, this value is chosen as T max = 0.1.

T max controls the local value of the maximum fluid viscosity, possibly resulting in the trun- cation of the Herschel-Bulkley law. The nominal incompressibility factor remains a free parameter, that should be chosen in the range T ∈]0, T max ]. Small values of T are expected to allow large viscosity variations. In contrast, viscosity truncation is expected to occur if T is too close to T max . In practice, the value of T l is controlled at each time step. In the present implementation, warning messages are sent to the simulation output when viscosity truncation occurs on one or more lattice nodes. Truncated simulations should be carefully considered: while small and localized viscosity truncations may have negligible effects on the overall flow solution, important truncations can significantly alter the flow accuracy, as illustrated hereafter. Therefore, it is generally recommended to avoid viscosity truncation (by decreasing T ) in order to ensure the numerical accuracy.

The effect of T is examined in the following in the case Bn = 1 and Re = 0.01. Figure 5(a) shows the evolution of the drag force coefficient over the range T ∈ [0.001, 0.1]. The drag force exhibits a constant evolution for T ≤ 0.01. In this region, no viscosity truncation occurs, i.e. the Herschel-Bulkley law is fully recovered in the limit of the floating point approximation (see §4.1). In addition, the predicted value of the drag is in agreement with the force computed by Nirmalkar et al. [START_REF] Nirmalkar | On creeping flow of a Bingham plastic fluid past a square cylinder[END_REF], which has been reported in the figure.

As expected, the case T = T max = 0.1 results in viscosity truncation. The effect of the truncation on the drag force is clearly noted in figure 5(a).

The viscosity patterns associated with the three simulations reported in figure 5(a) are depicted in figure 5(b) using iso-contours of the local viscous relaxation time. As expected from figure 5(a), similar patterns are observed for T = 0.001 and T = 0.01. The flow exhibits sharp interfaces between solid-and fluid -like regions, called yield surfaces [START_REF] Nirmalkar | On creeping flow of a Bingham plastic fluid past a square cylinder[END_REF]. These interfaces are associated with considerable viscosity ratios, extending over more than 10 orders of magnitude. In particular, substantially large viscosities appear in low-shear regions, as allowed by the present definition of µ max . The viscosity pattern is in agreement with previously reported simulation results The fluid region exhibits a figure-eight shape in the vicinity of the body. Three types of solid regions are observed: (i) the far-field region, where the fluid follows a solid motion, (ii) two isolated regions of solid rotation on the upper and lower sides of the body, and (iii) two triangular stagnation regions attached to the upstream and downstream sides of the cylinder. A major alteration of the flow is noted when the Herschel-Bulkley law is truncated, as illustrated by the case T = 0.1.

The described flow pattern has implications on the effect of the domain size on the flow solution. Indeed, the emergence of a solid -like region away from the body indicates that the flow is close to uniform in these regions. If the computational boundaries are far enough from the body surface to allow this fluid /solid transition, their effect on the flow is thus expected to be limited, even though the Reynolds number is small. This is confirmed by the domain size convergence study performed by Nirmalkar et al. [START_REF] Nirmalkar | On creeping flow of a Bingham plastic fluid past a square cylinder[END_REF].

The evolution of the drag force as a function of the Bingham number is examined in figure 6(a). The maximum incompressibility factor is maintained to T max = 0.1. As Bn increases, the global incompressibility factor has to be decreased in order to avoid viscosity truncation. The values T = 10 -2 , 10 -3 , 10 -4 , 10 -5 have been employed for the cases Bn = 1, 10, 100, 1000, respectively. Note that the relevant value of T can be estimated a priori, thus avoiding unnecessary numerical tests. Indeed, the typical value of T l can be expressed, using γ ≈ u 0 /D, as T l = T (1 + Bn). The condition T l < T max thus leads to T < T max /(1 + Bn), indicating that in the range Bn 1, T should be divided by 10 when Bn is multiplied by 10. The data issued from the numerical work of Nirmalkar et al. [START_REF] Nirmalkar | On creeping flow of a Bingham plastic fluid past a square cylinder[END_REF] have been reported in figure 6(a) for comparison purpose. An excellent agreement is noted between both studies.

The solid and fluid regions are determined by defining the yield surface as the iso-contour γ = , with 1. The value = 10 -3 is used in the following. The resulting patterns are plotted in figure 6(b) for three values of the Bingham number. The evolution of the yield surface is well-captured by the present simulations, and the predicted patterns are consistent with the numerical results of Nirmalkar et al. [START_REF] Nirmalkar | On creeping flow of a Bingham plastic fluid past a square cylinder[END_REF]. Overall, the size of the fluid region tends to decrease as a function of Bn. In contrast, solid regions are extended, leading to the formation of two equilateral triangular stagnation regions and two large solid -rotation regions at large Bingham numbers.

5.4. Creeping shear-thinning and shear-thickening flows past a square cylinder Simulations of shear-thinning and shear-thickening flows are considered in the following. The Bingham number is set to zero, and the simulations are performed over a range of the flow index n for Re = 0.01. Following the same procedure as that employed in the previous section, the maximum incompressibility parameter T max is set to 0.1. All simulations are performed with T = 0.01, a value that is small enough to avoid any viscosity cut-off over the range n ∈ [0.2, 1.8].

The evolution of the drag force coefficient is plotted in figure 7(a). The drag force significantly decreases as a function of the flow index. The numerical data of Pantokratoras [START_REF] Pantokratoras | Steady flow of a power-law non-Newtonian fluid across an unconfined square cylinder[END_REF] have also been reported in the figure. These data correspond to the unconfined simulations performed in that study, carried out with a domain size L = 50000D. A good agreement is observed between both studies in highly shear-thinning and highly shear-thickening regions. However, the drag force predicted by the present simulations is substantially higher when the flow index in close to one. These deviations are expected to relate to blockage effects, as Pantokratoras [START_REF] Pantokratoras | Steady flow of a power-law non-Newtonian fluid across an unconfined square cylinder[END_REF] pointed out that very large computational domains must be used in this case to achieve unconfined simulations. Additional computations with larger domains have been performed in order to examine this aspect. These simulations, also reported in figure 7(a), confirm that the present drag force tends to the one computed by Pantokratoras [START_REF] Pantokratoras | Steady flow of a power-law non-Newtonian fluid across an unconfined square cylinder[END_REF] as the domain size increases. The effect of the domain size is however contrasted over the range of flow indices, as important variations of the drag force are principally observed when n is close to 1. Figures 7(b,c) present visualizations of the flow in the vicinity of the cylinder for two values of the flow index, namely n = 0.2 and n = 1.8, corresponding to highly shear-thinning and highly shear-thickening cases. Figure 7(b) shows iso-contours of the viscous relaxation time. Overall, it can be noted that variations of the flow viscosity are much less pronounced than in the case of the yield-stress flow, visualized in figure 5. Shear-thinning and shearthickening simulations can thus be considered to be less critical regarding the effect of the fluid viscosity on the numerical accuracy. It should be noted that when n = 0.2, the flow tends to exhibit large viscosities in low-shear regions, i.e. the shear-thinning property does not systematically result in a viscosity decrease compared to the reference viscosity. These low-shear regions are encountered in the far-field region, but also close to the body near the upstream and downstream square surfaces. Similarly, low-viscosity regions emerge in the shear-thickening case. In both cases, viscosity variations are mostly found close to the body surface. The shear-thinning case however exhibits larger viscosity gradients; in particular, thin low-viscosity layers are noted on the upper and lower body surfaces. These smaller flow length scales can be clearly noted in figure 7(c), which shows iso-contours of the streamwise flow velocity. The important viscosity decrease close to the cylinder surface in the shear-thinning case results in remarkably thin boundary layers. Similar behaviors, especially characterized by a major increase of the wall vorticity when the flow index n decreases, have been observed in prior simulations of shear-thinning flows past bluff bodies [START_REF] Bharti | Steady flow of power law fluids across a circular cylinder[END_REF]. Consequently, the far-field flow velocity is rapidly recovered away from the body. This is consistent with the small blockage effects observed in the low-n region in figure 7(a). In contrast, the shear-thickening flow exhibits much smoother velocity variations, which might increase the effect of the domain size in this case. The important effect of n on the flow length scale is further illustrated in figure 8, which shows far-field iso-contours of the nondimensional streamwise flow velocity for n = 0.2 and n = 1.8. This evolution has important numerical implications, as either large computational domains or fine mesh resolutions may be required to achieve high numerical accuracy, depending on the value of the flow index.

Summary

This paper provides general guidelines for the simulation of highly viscous steady flows based on a general two-relaxation-time lattice-Boltzmann framework. The analysis of the space-continuous mesoscopic equations, performed through Hermite projection and multiscale expansion, shows that hydrodynamic behavior can be recovered at the macroscopic scale for any fluid viscosity, provided that the product of both relaxation times τ + τ -remains small enough. An alternative definition of the Knudsen number, controlling the magnitude of the non-equilibrium part of the mesoscopic solution, must however be considered, namely Kn * = √ τ + τ -c s /D. Its magnitude can be controlled by varying τ -which is a free numerical parameter. The present analysis supports previous works showing that bulk numerical solutions of the discrete TRT lattice-Boltzmann equation are viscosity-independent if Λ = τ + τ -, which controls the magnitude of numerical errors, is kept constant [START_REF] Humières | Viscosity independent numerical errors for lattice Boltzmann models: From recurrence equations to magic collision numbers[END_REF].

The main limitation concerning the fluid viscosity emerges when deriving the nondimensional macroscopic equations, which involves viscosity-dependent compressibility terms. In order to achieve incompressibility, the incompressibility parameter T = M a 2 /Re must satisfy T 1. A local definition of T is proposed to control the accuracy of the flow in simulations involving variable viscosities, i.e. non-Newtonian fluids. This analysis shows that considerably large fluid viscosities can be numerically allowed in low-shear regions, a feature that proved to be useful for yield-stress flow simulations.

These properties have been confirmed based on a series of numerical simulations of the creeping flow past a square cylinder. When the fluid is Newtonian, incompressibility of the flow is ensured for T = 0.1. This value of T should be used in future works to set the value of the Mach number in low-Reynolds-number simulations and to minimize the computational cost. Once T and Kn * are fixed, highly viscous simulations can be performed: relaxation times up to 10 4 are employed in the present work without any effect on the flow solution.

Non-Newtonian flows have been simulated using an Herschel-Bulkley fluid model. If the local value of T is properly controlled, major viscosity variations can be reproduced. Yieldstress flow solutions exhibit sharp yield-surfaces associated with viscosity ratios larger than 10 10 . Simulations have been performed over the range of Bingham numbers Bn ∈ [1, 1000], i.e. including regimes that had not been simulated using the LB method; the predicted flow solutions are in agreement with previous numerical results based on Navier-Stokes methods. More generally, this property may be expressed for any n-order moment P (n) = i f i c n i as This second-order projection is then simplified by linear combination with the first-order equation [START_REF] Gsell | Explicit and viscosity-independent immersed-boundary scheme for the lattice Boltzmann method[END_REF]. Multiplying [START_REF] Gsell | Explicit and viscosity-independent immersed-boundary scheme for the lattice Boltzmann method[END_REF]

P (n) = P (n)+ = i f + i c n i if n is even (A.
c iα f + i = 1 2 ∂ ∂x α (c iα f i ) - 1 2 ∂ ∂x α (c iα f i ) , (C.1)
which can be rewritten as a function of the collision operators using equation [START_REF] Lauga | Continuous breakdown of Purcells scallop theorem with inertia[END_REF],

∂ ∂x α c iα f + i = - 1 2τ + f neq+ i - 1 2τ -f neq- i + 1 2τ + f neq+ i + 1 2τ -f neq- i .
(C.2)

Noticing that f + i = f + i and f - i = -f - i , this equation becomes Developing f + i as a function of f i and f - i , namely

∂ ∂x α c iα f + i = - 1 τ -f neq-
f + i = f i -f - i , gives -f neq i = τ -∂ ∂x α (c iα f i ) + (τ --τ + ) ∂ ∂x α c iα f + i . (C.6)
A similar equation can be obtained by developing f - i as a function of f i and f + i ,

-f neq i = τ + ∂ ∂x α (c iα f i ) + (τ + -τ -) ∂ ∂x α c iα f - i . (C.7)
Finally, the non-equilibrium distribution functions can be expressed as functions of f

∂ ∂x α (c iα f i ) = τ + -τ - τ + + τ - ∂ ∂x α (c iα f i ) , (C.10a) ∂ ∂x α (c iα f i ) = τ + + τ - τ + -τ - ∂ ∂x α (c iα f i ) , (C.10b)
which is equivalent to τ + -τ - τ + + τ -= τ + + τ - τ + -τ -⇐⇒ τ + τ -= 0. (C.11)

Figure 1 :

 1 Figure 1: Schematic view of the physical configuration: a square cylinder of size D is immersed in a flow characterized by its velocity u 0 , density ρ 0 , consistency k, flow index n and yield stress σ 0 .

Figure 2 :

 2 Figure 2: Flow past a square cylinder for Re ≥ 1: evolution of the drag force coefficient C x as a function of (a) the mesh resolution D/∆n and (b) domain size L/D for Re = 1, and (c) as a function of the Reynolds number for D = 21∆n and L = 50D.

Figure 3 :Fig. 3 Re = 10 - 2

 33102 Figure 3: Effect of the incompressibility factor on Newtonian creeping flow simulations: (a) evolution of the drag force coefficient as a function of T , for two values of the Reynolds number Re 1, and iso-contours of the non-dimensional (b) streamwise flow velocity and (c) shear-rate magnitude close to the cylinder, for Re = 10 -2 and for three values of T indicated by the labels in (a).

Figure 4 :

 4 Figure 4: Effect of the Reynolds number Re 1 on Newtonian creeping flow simulations, for a fixed incompressibility factor T = 0.1: (a) evolution of the drag force coefficient C x,µ as a function of Re, and cross-flow evolution of the non-dimensional (b) streamwise velocity and (c) shear-rate magnitude close to the upper square surface at x/D = 0, for two values of the Reynolds number.

3 Figure 5 :

 35 Figure 5: Effect of the global incompressibility parameter T = M a 2 /Re on yield-stress flow simulations, for Re = 0.01 and Bn = 1: evolution of the yield drag coefficient C x,σ0 as a function of T , and (b) isocontours of the normalized viscous relaxation time for three values of T , indicated by the labels reported in (a). Note that the maximal local incompressibility factor T max is set to 0.1 (see §4.3). In (a), the horizontal dashed line indicates the drag value reported by Nirmalkar et al. [40].

22

 22

3 Figure 6 :

 36 Figure 6: Effect of the Bingham number Bn on yield-stress flow simulations, for Re = 0.01: (a) evolution of the yield drag coefficient C x,σ0 as a function of Bn, predicted by the present simulations and reported by Nirmalkar et al. [40], and (b) visualization of the yield surface for three values of Bn, indicated by the labels reported in (a).

Figure 7 :

 7 Figure 7: Effect of the flow index n of shear-thinning and shear-thickening simulations: (a) evolution of the drag coefficient as a function of n, for different computational domain sizes L and compared to the unconfined numerical results reported by Pantokratoras [39], and iso-contours of (b) the normalized relaxation time and (c) the non-dimensional streamwise flow velocity in the vicinity of the body in the case L = 100D, for n = 0.2 and n = 1.8.

Figure 8 :

 8 Figure 8: Effect of n on the flow length scale at Re = 0.01: far-field iso-contours of the non-dimensional streamwise flow velocity for (bottom) n = 0.2 and (top) n = 1.8. The size of the computational domain is L = 100D.

.

  Shear-thinning and shear-thickening flows are also accurately reproduced over a wide range of flow indices n ∈ [0.2, 1.8]. It should be recalled that the following identities are also satisfied by these functions: of f i are fully recovered by the symmetric populations at even orders and by the anti-symmetric ones at odd orders. At orders 0 and 1, this is expressed as

8 )

 8 Appendix B. Second-order Hermite projection of the TRT Boltzmann equationThe second-order projection of the TRT Boltzmann equation reads∂ ∂x α (Q αβγ -ρc 2 s u α δ βγ ) = -1 τ + Π neq+ βγ -1 τ -Π neq- βγ , (B.1)where Q αβγ is the third-order moment defined byQ αβγ = i f i c i,α c i,β c i,γ and Π neq+ βγ and Π neq- βγ are the symmetric and anti-symmetric parts of Π neq βγ . It can be shown that, due to symmetry (c i,β c i,γ = c i,β c i,γ ), Π neq- βγ = 0 and Π neq+ βγ β c i,γ = Π neq βγ . (B.2b)The equilibrium part of Q αβγ can be computed using the equilibrium distributions, namely Q eq αβγ = ρu α u β u γ + ρc 2 s (u α δ βγ + u β δ αγ + u γ δ αβ ). Equation (B.1) thus becomes ∂ ∂x α (ρu α u β u γ + ρc 2 s (u β δ αγ + u γ δ αβ ) + Q neq αβγ ) = -1 τ + Π neq βγ . (B.3)

4 )

 4 can be shown that the dynamics of the antisymmetric populations is governed by ∂ ∂x α c iα f - Summing equations (C.3) and (C.4), the non-equilibrium functions satisfy-f neq i = τ -∂ ∂x α c iα f + i + τ + ∂ ∂x α c iα f -

  i and f i by summing equations (C.6) and (C.7),-2f neq i = (τ + + τ -) ∂ ∂x α (c iα f i ) -(τ + -τ -) ∂ ∂x α (c iα f i ) . (C.8)On the other hand, non-equilibrium populations f neq i also satisfy-2f neq i = (τ + -τ -) ∂ ∂x α (c iα f i ) -(τ + + τ -) ∂ ∂x α (c iα f i ) . (C.9)When the non-equilibrium functions vanish, namely f neq i

Table 1 :

 1 Details on the numerical setups analyzed in figures 3 and 4.

			10 -3 3.16 × 10 -3 3.16 × 10 -1	12	8.27
			10 -2	10 -2	1	36.9	8.27
			10 -1 3.16 × 10 -2	3.16	115.4 8.25
			1	10 -1	10	364.2 8.13
			10 -3	10 -3	1	36.9	8.26
		Re = 10 -3	10 -2 3.16 × 10 -3 10 -1 10 -2	3.16 10	115.4 8.26 364.2 8.25
			1	3.16 × 10 -2	31.6	1149.9 8.13
	Fig. 4	Re = 10 -4 0.1 3.16 × 10 -3	31.6	1149 8.25
		Re = 10 -5 0.1	10 -3	100	3638 8.25
		Re = 10 -6 0.1 3.16 × 10 -4	316	11494 8.25

  by u γ gives∂ ∂x α (ρu α u β u γ ) -ρu α u βA similar expression can be obtained by expressing[START_REF] Gsell | Explicit and viscosity-independent immersed-boundary scheme for the lattice Boltzmann method[END_REF] on indices α and γ and multiplying it by u β , namely∂ ∂x α (ρu α u β u γ ) -ρu α u γ (ρu α u β u γ) due to the mass conservation equation (14),Appendix C. Equilibrium solutions of the TRT Boltzmann equationFirst, the dynamic equation for the symmetric populations f + i is derived using ∂ ∂x α

	Noticing that ρu α and defining the shear-rate tensor S αβ = ∂ ∂x α ∂ (u β u γ ) = ∂x α	1 2	(	∂u α ∂x β	+	∂u β ∂x α	), equation (B.6) simplifies to
	2ρc 2 s S βγ +	∂ ∂x α	(Q neq αβγ ) -u γ	∂ ∂x α	(Π neq αβ ) -u β	∂ ∂x α	(Π neq αγ ) =	-Π neq βγ τ + .	(B.7)
				∂ ∂x α	(u γ ) + u γ	∂ ∂x β	(ρc 2 s ) + u γ	∂ ∂x α	(Π neq αβ ) = 0.	(B.4)
				∂ ∂x α	(u β ) + u β	∂ ∂x γ	(ρc 2 s ) + u β	∂ ∂x α	(Π neq αγ ) = 0.	(B.5)
	Subtracting (B.4) and (B.5) to (B.3), a new second-order projection is obtained:
	-+	∂ ∂x α ∂ ∂x α	(ρu α u β u γ ) + ρu α (Q neq αβγ ) -u γ ∂ ∂x α (Π neq ∂ ∂x α (u β u γ ) + ρc 2 s αβ ) -u β ∂ ∂x α (Π neq ∂u γ ∂x β αγ ) = + -Π neq ∂u β ∂x γ βγ τ + .	(B.6)

Acknowledgements

We thank Song Zhao for useful discussions. The project MACBION leading to this publication has received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French Investissements d'Avenir programme and the SINUMER project (ANR-18-CE45-0009-01) of the French National Research Agency (ANR).

Appendix A. Symmetric and anti-symmetric population properties

Equilibrium and non-equilibrium populations are decomposed into symmetric and antisymmetric parts,

and