Modélisation thermodynamique des acides
Modélisation thermodynamique des acides
LIENS
THÈSE
présentée et soutenue publiquement le 26 avril 2022
pour l’obtention du
Composition du jury
Je tiens à remercier Messieurs Lhachmi Khamar et Saad Benjelloun pour leur co-encadrement et collabo-
ration, de m’avoir accueilli dans leur équipe de recherches au laboratoire MSDA. Je tiens à leur exprimer
ma profonde gratitude pour leur aide, conseils, et accompagnement. Merci beaucoup pour m’avoir ac-
compagné pendant ces années aussi bien professionnellement qu’à amicalement.
Mes remerciements vont aussi à Monsieur Brahim Benyahia et à Madame Isabel Pitault d’avoir accepté
de rapporter sur cette thèse. Je voudrais les remercier d’avoir pris le temps de la lecture et des correc-
tions du manuscrit. Faire partie du jury de ma soutenance m’honore beaucoup. Un remerciement spécial
à Professeur Brahim Benyahia pour tout le soutien qu’il a mis en œuvre pour mon après thèse.
Je remercie Madame Caroline gentric d’avoir accepté de présider le jury, Madame Lucie Jaubert et Ma-
dame Sara El Bouzidi d’avoir pris le temps pour examiner le manuscrit, et surtout, d’avoir accepté de
faire partie du jury. Votre présence me fait beaucoup d’honneur.
Un remerciement très spécial à mes parents Hassan et Rabia qui m’ont soutenu pendant tout ce cursus,
je voudrais leur exprimer ma profonde gratitude pour tout l’amour et l’affection qu’ils m’ont témoignés.
Je voudrais leur dédier ces travaux et j’espère que devenir docteur leur fera honneur et les comblera de
bonheur. Un grand merci à mon frère Zakaria et mes deux sœurs Kamilia et Saoussane, merci pour tout
l’amour et l’accompagnement moral que vous m’avez témoignés. Je suis très reconnaissant que nos liens
familiaux soient assez forts, et qu’on se retrouve l’un à coté de l’autre dès que c’est nécessaire.
Je voudrais dédier ce travail aussi à mon frère et cher ami d’enfance Yahya El Fariq, et à sa famille qui
m’ont également soutenu à plusieurs reprises. On a commencé ce cursus ensemble, et la vie nous a donnés
enfin deux chemins différents à suivre. Je voudrais te remercier pour tout les fous rires, l’accompagnement
et les conseils. Je suis toujours fier de toi, et j’espère que ces quelques mots t’expériment ma profonde
gratitude.
J’arrive maintenant à mes amis avec qui j’ai passé toute ces années de recherche. Je commence par Cris-
tian Cardenas « Les choses de France restent en France » et Alexis Courtais « J’ai peut-être une idée mon
gars » qui m’ont accueilli et intégré dans l’équipe assez facilement. Grâce à vous, j’ai pu découvrir plein
de choses et surtout de cultures, je n’oublierai jamais tous les fous rires et les fous moments partagés. Mes
remerciements vont à mes amis du LRGP aussi, Fatima Matamoros, Andres Pina, Christian Quintero,
Felipe Buendia, Thoma Delpietro, Nihad Kamar, Chakib Behloul, Jennifer Andrea. Je souhaite remercier
mes amis du labo au Maroc, Mohamed Hamou, Marwan El Karii, Abdelouahed Ouardghi, Moussa Zig-
gaf, Mustafa Bahari, Amine Hamadi, Mohamed El Guide, Souhail Maazioui, Sanae et Safae el Missaoui,
Amine Hanini. La liste est très longue.
i
ii
À mes parents Hassan et Rabia,
Cette thèse vous appartient.
Mes mots ne seraient jamais à la hauteur de l’amour
et l’affection que vous m’avez témoignés tout au
long de ce cursus. J’aimerais vous exprimer toute
ma gratitude et ma reconnaissance. Cette dédicace
est pour moi, la meilleure façon de vous honorer
et vous montrer à quel point vous avez été magnifique.
iii
iv
Sommaire
Nomenclature xvii
Introduction générale 1
Chapitre 1
Généralités et étude bibliographique
v
Sommaire
Chapitre 2
Modélisation thermodynamique du procédé de fabrication d’acide phos-
phorique dihydrate
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2 Revue de la littérature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.3 Modélisation thermodynamique . . . . . . . . . . . . . . . . . . . . . . . . 34
2.3.1 Procédé dihydrate . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.3.2 Bilans de matière et de charge et équilibres . . . . . . . . . . . . . . 35
2.4 Base de données expérimentales . . . . . . . . . . . . . . . . . . . . . . . . 39
2.4.1 Mesures effectuées en laboratoire dans ce travail . . . . . . . . . . . 39
2.4.2 Mesures de la littérature . . . . . . . . . . . . . . . . . . . . . . . . 39
2.5 Analyse d’estimabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.5.1 Approche locale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.5.2 Approche globale . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.6 Identification des paramètres . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.7 Résultats et discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
2.7.1 Analyse de sensibilité globale . . . . . . . . . . . . . . . . . . . . . 55
2.7.2 Analyse d’estimabilité . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.7.3 Identification des paramètres et seuil optimal d’estimabilité . . . . . 58
2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Chapitre 3
Exploitation du modèle thermodynamique pour l’optimisation multicri-
tère du procédé dihydrate
vi
3.4 Résultats de l’optimisation . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.4.1 Approche de résolution . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.4.2 Front de Pareto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.4.3 Méthode d’aide à la décision . . . . . . . . . . . . . . . . . . . . . . 75
3.4.4 Modèle de substitution . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
Chapitre 4
Modélisation thermodynamique des problèmes d’encrassement dans la
cuve d’attaque
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Encrassement des unités de production d’acide phosphorique . . . . . . . . 83
4.2.1 Impuretés liées aux procédés . . . . . . . . . . . . . . . . . . . . . . 84
4.2.2 Impuretés liés au minerai de phosphate . . . . . . . . . . . . . . . . 84
4.3 Modélisation thermodynamique . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3.1 Indice de saturation . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.3.2 Phases minérales considérées . . . . . . . . . . . . . . . . . . . . . . 86
4.4 Résultats et discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
Chapitre 5
Modélisation de la cristallisation du gypse dans un réacteur batch
vii
Sommaire
Chapitre 6
Validation expérimentale du modèle thermodynamique
Bibliographie
viii
Table des figures
ix
Table des figures
2.10 Valeurs de Jcc,L pour chaque cas mse. (A) : activité de l’eau, (B) : molalités
des espèces, (C) : activité de l’eau et molalités des espèces. . . . . . . . . . 59
4.1 Prédiction des indices de saturation des différentes phases minérales consi-
dérées. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.3 Distribution des tailles des particules normées à différents instants. . . . . 108
5.5 Distribution des tailles des particules (normée) dans le temps. . . . . . . . 110
x
5.6 Résultats des différences finies centrées (2D). (A) : Concentrations du phos-
phate tricalcique et d’acide sulfurique. (B) : Concentrations du gypse et
d’acide phosphorique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.7 Distributions des tailles des particules normées en fonction du temps. (A1-
B1) : à t=0 min. (A2-B2) : à t=10 min. (A3-B3) : à t=30 min. (A4-B4) :
à t=100 min. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.1 Montage expérimental. (1) : Bain thermostaté. (2) : Agitateur. (3) : Réac-
teur batch. (4) : Thermostat. (5) : Ampoule d’acide sulfurique. (6) : Pompe
à vide. (7) : Dispositif de filtration. (8) : Hotte de laboratoire. . . . . . . . 118
6.2 Composition chimique des phases liquides produites . . . . . . . . . . . . . 127
6.3 Comparaison des prédictions du modèle développé avec les mesures . . . . 128
6.4 Comparaison des prédictions du modèle développé avec les mesures . . . . 129
6.5 Comparaison des prédictions du modèle développé avec les mesures . . . . 130
6.6 Test de Kolmogorov-Smirnov pour vérifier la normalité des échantillons
utilisés pour la validation statistique. Les échantillons sont centrés et réduits.131
6.7 Résultats de diffraction aux rayons X des échantillons solides. (A) : expé-
rience 4. (B) : expérience 7. (C) : expérience 9. . . . . . . . . . . . . . . . . 135
6.8 Résultats de diffraction aux rayons X des échantillons solides. (A) : expé-
rience 10. (B) : expérience 13. (C) : expérience 17. . . . . . . . . . . . . . . 136
6.9 Résultats de diffraction aux rayons X des échantillons solides. (A) : expé-
rience 20. (B) : expérience 24. (C) : expérience 26. . . . . . . . . . . . . . . 137
6.10 Différentes zone de l’interface graphique développée . . . . . . . . . . . . . 140
6.11 Résultats de simulation de l’attaque du minerai de phosphate à l’aide de
l’interface graphique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
xi
Table des figures
xii
Liste des tableaux
xiii
Liste des tableaux
6.1 Description des réactifs utilisés pour la réalisation des expériences . . . . . 119
6.2 Description du premier plan d’expériences . . . . . . . . . . . . . . . . . . 120
6.3 Description du deuxième plan d’expériences . . . . . . . . . . . . . . . . . 120
6.4 Conditions opératoires pour le premier plan d’expériences . . . . . . . . . . 121
6.5 Conditions opératoires pour le deuxième plan d’expériences . . . . . . . . . 122
6.6 Analyses effectuées des deux phases liquide et solide . . . . . . . . . . . . . 123
6.7 Validation du modèle developpé par le test F de Fisher-Snedecor. . . . . . 133
6.8 Conditions opératoires et indices de saturation du gypse. . . . . . . . . . . 134
xiv
8 Paramètres d’interactions de Pitzer pour le système AlCl3 − H2 O. . . . . 149
9 Paramètres d’interactions de Pitzer pour le système CaCl2 − H2 O. . . . . 149
10 Paramètres d’interactions de Pitzer pour le système K2 SO4 − M gSO4 −
H2 O. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
11 Paramètres d’interactions de Pitzer pour le système F eSO4 −H2 SO4 −H2 O. 150
12 Paramètres d’interactions de Pitzer pour le système N aSO4 − N aCl − H2 O. 151
13 Paramètres d’interactions de Pitzer pour le système M gSO4 − M gCl2 −
H2 O. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
xv
Liste des tableaux
xvi
Nomenclature
Symboles latins
aw Activité de l’eau -
ai Activité de l’élément i -
Aφ Constante de Debye-Hückel kg · mol−1/2
b Paramètre universel de Debye-Hückel kg1/2 · mol−1/2
Cp Capacité thermique J.K −1
Cpid Capacité thermique idéale J.K −1
Cpex Capacité thermique d’excès J.K −1
Cp0 Capacité thermique de référence J.K −1
(0)
Ci,j Paramètre d’interaction binaire de Pitzer -
(1)
Ci,j Paramètre d’interaction binaire de Pitzer -
Cpas Capacité thermique d’acide sulfurique J.K −1
Cpps Capacité thermique d’une suspension de phosphate J.K −1
Cpout Capacité thermique d’un flux de sortie J.K −1
Cprph Capacité thermique de l’acide phosphorique de retour J.K −1
Cpw0 Capacité thermique de l’eau à l’état de référence J.K −1
d Nombre de paramètre du modèle de Pitzer -
Dn Statistique du test de Kolmogorov-Smirnov -
Dn∗ Statistique critique du test de Kolmogorov-Smirnov -
f Terme de Debye-Hückel -
Gi Paramètre d’interaction du modèle de Bromley -
G Énergie libre de Gibbs J
Gex Énergie de Gibbs d’excès J
Gid Énergie libre de Gibbs d’une solution idéale
Gx Terme de croissance dans la direction de la longueur m/s
Gy Terme de croissance dans la direction de la largeur m/s
H Enthalpie d’une solution J
H0 Hypothèse nulle de test statistique -
H id Enthalpie d’une solution idéale J
xvii
Nomenclature
xviii
Trph Température de l’acide phosphorique de retour K
Tout Température d’un flux de sortie K
UL Pourcentage des pertes inattaquées %
Vjk Paramètre d’interaction du modèle de Bromley -
V Volume molaire m3
Vps Volume d’une suspension de phosphate m3
Vrph Volume de l’acide phosphorique de retour m3
Vout Volume d’un flux de sortie m3
V id Volume idéal d’une solution m3
V ex Volume d’excès d’une solution m3
V0 Volume de référence d’une solution m3
Vm Volume molaire m3
w1 Pourcentage d’acide sulfurique en section 1 %
w2 Pourcentage d’acide sulfurique en section 2 %
w3 Pourcentage d’acide sulfurique en section 3 %
xj Fraction molaire du composant j n%
Z Matrice d’estimabilité -
Z̃ Statistique du test F de Fisher-Snedecor -
Z∗ Statistique critique du test F de Fisher-Snedecor -
zi Grandeur molaire partielle -
zi Charge de l’élément i -
Zm Charge de cation M -
Zx Charge d’anion X -
Z ex Grandeur molaire partielle d’excès -
Ziid Grandeur molaire partielle d’une solution idéale -
Symboles grecs
xix
Nomenclature
Abréviation
CP U Temps de calcul
CSD Distribution de taille des particules
DRX Diffraction aux rayons X
exp Expérience
mod Modèle
mse Erreur quadratique moyenne
NA Nombre d’anions
NC Nombre de cations
NM Nombre de molécule
P BE Modèle de bilan de population
P DH Pitze-Debye-Hückel
RM SE Racine de l’erreur quadratique moyenne
SI Indice de saturation
tot Total
xx
Introduction générale
1
Contexte de l’étude
De nos jours, la demande mondiale en phosphates et en dérivés phosphatés est en constante
augmentation. En effet, suite à un accroissement démographique très élevé qu’a connu le
monde ces dernières décennies, il devient crucial d’augmenter les surfaces agricoles et sur-
tout d’améliorer le rendement agricole afin d’assurer une suffisance alimentaire. Une des
solutions les plus prometteuses pour atteindre cet objectif est l’utilisation des engrais,
en particulier les engrais phosphatés. Toutefois, il est impératif de mieux exploiter ces
ressources naturelles, notamment en optimisant les procédés de leur valorisation.
Parmi les problèmes qui se posent, les pertes chimiques des phosphates lors de la réaction
de digestion du minerai méritent d’être analysées et étudiées afin d’améliorer davantage
les performances du procédé. En effet, trois types de pertes chimiques se distinguent. Il
s’agit (i) des pertes chimiques inattaquées qui sont dues à l’enrobage des particules de
phosphates par les sulfates de calcium produits les rendant ainsi inaccessibles par l’acide
d’attaque. (ii) Des pertes syncristallisés qui sont dues à l’insertion des espèces phosphatées
dans le réseau cristallin des sulfates de calcium qui sont considérés comme produit indési-
rable. (iii) Des pertes solubles en eau qui sont dues à l’évacuation de l’acide phosphorique
avec le gâteau de filtration des sulfates de calcium, lorsque les qualités de filtrabilité et de
lavabilité du gâteau sont dégradées.
Par ailleurs, la présence de plusieurs impuretés, qui sont issues du minerai de phosphates
ou de l’eau du procédé non traitée, sont à l’origine de plusieurs réactions secondaires qui
se produisent simultanément avec la réaction principale de l’attaque. La plupart de ces
réactions sont néfastes pour les performances du procédé. En effet, elle conduisent à la
production de plusieurs minéraux indésirables qui causent l’encrassement des équipements
du procédé, tels que les réacteurs, les échangeurs de chaleur, les conduites etc., impactant
ainsi directement leur durée de vie.
3
Introduction générale
Objectifs de la thèse
L’objectif principal de cette thèse est de parfaire les connaissances techniques et scien-
tifiques sur le fonctionnement des procédés de production d’acide phosphorique par voie
humide.
1. Un modèle thermodynamique fondé sur les phénomènes mis en jeu et décrit par des
équations de bilans de matière, de charges, et d’équilibre sera développé (connais-
sances).
2. Ce modèle mettra inévitablement en jeu des paramètres inconnus dont les valeurs
appropriées ne sont pas toujours disponibles. La question sera alors de savoir s’ils
(tous ou seulement une partie) peuvent être estimables à partir des mesures expé-
rimentales disponibles. D’où l’intérêt de l’analyse d’estimabilité des paramètres, et
de méthodes de conception d’expériences optimales pour éventuellement effectuer le
minimum d’expériences (souvent très coûteuses) qui contiennent le maximum d’in-
formations (informations).
3. Pour que le modèle soit validé, il doit être confronté à la réalité de l’expérience.
Ainsi, les paramètres inconnus et estimables du modèle seront d’abord identifiés à
partir de mesures disponibles, et ensuite les prédictions du modèle identifié seront
comparées aux mesures obtenues dans des conditions opératoires différentes de celles
qui ont servi pour son identification. Il est par conséquent clair que la qualité et la
quantité des mesures expérimentales sont capitales pour la compréhension et la mo-
délisation des phénomènes (données).
4
Plan du manuscrit
Le manuscrit de cette thèse est organisé autour des chapitres suivants :
Le troisième chapitre présente une première exploitation du modèle développé dans l’op-
timisation multicritère du procédé. Les fonctions objectifs sont ainsi prédites à l’aide du
modèle, et les résultats d’optimisation sont comparés à des expériences effectuées au ni-
veau industriel.
Le cinquième chapitre présente une troisième exploitation du modèle. Il traite plus spé-
cifiquement de la modélisation de la cristallisation des sulfates de calcium durant la pro-
duction de l’acide phosphorique. Le modèle de cristallisation fait intervenir des termes
sources qui sont prédits par le modèle.
Pour finir, le manuscrit se termine par une conclusion générale sur le travail réalisé ainsi
que quelques propositions de perspectives pour de futures études.
5
Introduction générale
6
Chapitre 1
Sommaire
1.1 Grandeurs thermodynamiques . . . . . . . . . . . . . . . . . . . 9
1.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.1.2 Grandeurs extensives . . . . . . . . . . . . . . . . . . . . . . . . 9
1.1.3 Énergie libre de Gibbs . . . . . . . . . . . . . . . . . . . . . . . 10
1.1.4 Potentiel chimique . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.5 État de référence et conditions standards . . . . . . . . . . . . 12
1.1.6 Solution idéale . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.1.7 Propriétés thermodynamiques d’excès . . . . . . . . . . . . . . 12
1.1.8 Activité de l’eau - coefficient osmotique . . . . . . . . . . . . . 13
1.1.9 Coefficient d’activité . . . . . . . . . . . . . . . . . . . . . . . . 13
1.1.10 Volume de solutions aqueuses . . . . . . . . . . . . . . . . . . . 14
1.1.11 Masse volumique des solutions aqueuses . . . . . . . . . . . . . 15
1.1.12 Capacité thermique - Enthalpie . . . . . . . . . . . . . . . . . . 15
1.1.13 Solubilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.2 Modèles thermodynamiques . . . . . . . . . . . . . . . . . . . . 18
1.2.1 Théorie de Debye–Hückel . . . . . . . . . . . . . . . . . . . . . 18
1.2.2 Modèle de Bromley - Zemaitis . . . . . . . . . . . . . . . . . . . 20
1.2.3 Modèle de Wilson . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.2.4 Modèle eUNIQUAC . . . . . . . . . . . . . . . . . . . . . . . . 23
1.2.5 Modèle LIQUAC . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.2.6 Modèle de Pitzer . . . . . . . . . . . . . . . . . . . . . . . . . . 25
1.2.7 Modèle eNRTL . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
7
Chapitre 1. Généralités et étude bibliographique
Résumé
8
1.1. Grandeurs thermodynamiques
Ce chapitre est organisé en deux parties. La première partie est consacrée à une étude
bibliographique sur les grandeurs thermodynamiques pertinentes pour le développement
de modèles de procédés. Quant à la deuxième, elle décrit les principaux modèles thermo-
dynamiques électrolytiques utilisés pour les prédire.
1.1.1 Introduction
Les modèles phénoménologiques mettent très souvent en jeu des grandeurs physiques dont
la qualité de prédiction est capitale pour l’optimisation de la conception, du dimension-
nement et du fonctionnement des procédés. Les grandeurs nécessaires pour les modèles
développés dans le ce travail sont définies ci-dessous.
m
X ni
Z(T, P, ~n) = ni zi (T, P, ~x), ~x = (x1 , ..., xm ), ~n = (n1 , ..., nm ), xi = Pm (1.1)
i=1 i ni
avec :
∂Z(T, P, ~n)
zi (T, P, ~x) = (1.2)
∂ni T,P,nj6=i
m
X m
X
dZ(T, P, ~n) = ni dzi + zi (T, P, ~x)dni (1.3)
i=1 i=1
9
Chapitre 1. Généralités et étude bibliographique
ou encore en fonction de la somme des dérivées partielles (Lach, 2015 ; Pitzer, 1973) :
∂Z (T, P, ~n) ∂Z (T, P, ni )
dZ (T, P, ~n) = dT + dP
∂T P,~
n ∂P T,~
n
(1.4)
m
X ∂Z (T, P, ~n)
+ dni
i=1
∂ni T,P,nj6=i
L’expression ci-dessus est générale et peut être utilisée pour décrire les variations de
plusieurs propriétés extensives. Nous l’exploiterons dans ce qui suit pour le calcul de
l’énergie libre de Gibbs (Baker et al., 1982) qui constitue la base de plusieurs modèles
thermodynamiques.
(1.6)
m
X ∂G (T, P, ~n)
+ dni
i=1
∂ni T,P,nj6=i
avec :
∂G (T, P, ~n)
= −S (1.8)
∂T P,~
n
10
1.1. Grandeurs thermodynamiques
∂G (T, P, ~n)
=V (1.9)
∂P T,~
n
m
X ∂G (T, P, ~n)
dG (T, P, ~n) = −SdT + V dP + dni (1.10)
i=1
∂ni T,P,nj6=i
G dG G
d = − dT (1.11)
RT RT RT 2
m
G S V 1 X ∂G (T, P, ~n) H S
d =− dT + dP + dni − 2
dT + dT (1.12)
RT RT RT RT i=1 ∂ni RT RT
soit :
m
G V H 1 X ∂G (T, P, ~n)
d = dP − dT + dni (1.13)
RT RT RT 2 RT i=1 ∂ni T,P,nj6=i
Il est à noter que l’énergie libre de Gibbs est une des propriétés thermodynamiques les plus
importantes, son calcul ou son estimation via des modèles thermodynamiques permet de
déterminer plusieurs propriétés de mélanges, ou d’espèces chimiques en solution, comme
l’entropie, le volume, le potentiel chimique, qui à leur tour permettent de déduire d’autres
propriétés micro et/ou macroscopiques (Lach, 2015 ; Zurich, 2018).
Le potentiel chimique peut être calculé soit à partir de l’énergie libre de Gibbs, ou à partir
de l’activité des espèces dans le milieu étudié (Kaufman, 2002 ; Mills et al., 1993) :
∂G (T, P, ~n)
µi (T, P, ~n) = (1.14)
∂ni T,P,nj6=i
11
Chapitre 1. Généralités et étude bibliographique
µi (T, P, ~n) = µ0aq,i (T, P ) + [Link] (ai (T, P, ~n)) , avecai = xi γi (1.15)
où µ0aq,i (T, P ) est le potentiel chimique de l’espèce i pure à l’état de référence à la tempé-
rature T et à la pression P . R est la constante universelle des gaz parfaits et ai et γi sont
respectivement l’activité et le coefficient d’activité de l’espèce i.
µid 0
aq,i (T, P, mi ) = µaq,i (T, P ) + RT ln (mi ) (1.16)
m
X
id
G (T, P, ~n) = ni µid
i (T, P, mi ) (1.17)
i=1
12
1.1. Grandeurs thermodynamiques
P
Mw i mi
lnaw (T, P, ~n) = φ(T, P, ~n) (1.20)
1000
∂Gex
− ∂ww
ni
φ−1= P (1.21)
RT i mi
ww est la masse d’eau dans la solution étudiée. Deux types de mesures permettent d’obte-
nir l’activité de l’eau ou le coefficient osmotique d’une solution. Le premier type concerne
les mesures indirectes qui nécessitent la connaissance de l’activité de l’eau d’une solution,
appelée solution de référence. Il s’agit par exemple des méthodes hygrométriques (El
Guendouzi et Errougui, 2005 ; El Guendouzi et Marouani, 2003) ou isopestiques
(Rard, 2018). Le second type de mesures concerne les mesures directes qui se basent
sur la mesure de la pression de vapeur (Lach, 2015 ; Shiah et Tseng, 1996) ou sur
l’abaissement du point de congélation (Ge et Wang, 2009).
13
Chapitre 1. Généralités et étude bibliographique
µex
aq,i (T, P, ~
n) = RT ln γi (T, P, ~n) (1.22)
Par conséquent :
Les coefficients d’activité sont alors obtenus en calculant la dérivée partielle de l’énergie
libre d’excès de Gibbs par rapport au nombre de moles. Ils sont très difficiles à mesurer
et il est plus simple de les estimer en utilisant des modèles thermodynamiques.
V ∂(G/RT )
= (1.25)
RT ∂P T,~
n
Soit :
∂Gid ∂Gex
V = + (1.26)
∂P T,~
n ∂P T,~
n
| {z } | {z }
V id V ex
m m
∂µid ∂µ0i
X X
V id
= ni i
= ni (1.27)
∂P T,~
n ∂P
i=1 i=1
| {z T,~n}
Vi0
14
1.1. Grandeurs thermodynamiques
Par conséquent, le volume d’une solution idéale est égal au volume résultant de la contri-
bution de tous ses constituants à l’état de référence. Vi0 est le volume de l’espèce i à l’état
de référence. Pour l’eau, on travaille à l’état de corps pur, et pour les solutés à l’état de
dilution infinie. On peut donc écrire :
X
V = nw .Vw0 + ni Vi0 + V ex (1.28)
n6=w
m m
ρ= = (1.29)
nw .Vw + n6=w ni Vi0 + V ex
0
P
V
avec :
X
m = nw Mw + ni Mi (1.30)
n6=w
Soit :
P
nw Mw + n6=w ni Mi
ρ= (1.31)
nw .Vw0 + n6=w ni Vi0 + V ex
P
∂H
Cp = (1.32)
∂T P,~
n
15
Chapitre 1. Généralités et étude bibliographique
avec :
H ∂(G/RT )
− = (1.33)
RT 2 ∂T P,~
n
Soit :
∂H id ∂H ex
Cp = + (1.35)
∂T P,~
n ∂T P,~
n
| {z } | {z }
Cpid Cpex
Gid m id
2 ∂ RT 2 ∂(µi /RT )
X
H id
= −RT = ni −RT (1.36)
∂T i=1
∂T
m 0
2 ∂(µi /RT )
X
H id
= ni −RT (1.37)
i=1
∂T
L’enthalpie d’une solution idéale est donc égale à la somme des enthalpies de tous ses
constituants dans l’état de référence. Ce constat permet de calculer la capacité thermique
d’une solution idéale comme suit :
m m
X ∂h0i X
Cpid = ni = 0
ni Cp,i (1.38)
i=1
∂T i=1
où Cp,i
0
est la capacité calorifique de l’espèce i en J/K/mol dans l’état de référence. Il
s’agit pour l’eau du corps pur, et pour les solutés de la dilution infinie. On peut donc
réécrire l’équation 1.35 :
X
0
Cp = nw Cpw + 0
ni Cp,i + Cpex (1.39)
n6=w
où Cpw
0
et Cp,i
0
correspondent à la capacité calorifique molaire à pression constante de l’eau
pure et à la capacité calorifique molaire de l’espèce i à dilution infinie.
16
1.1. Grandeurs thermodynamiques
1.1.13 Solubilité
La solubilité est la capacité d’une substance, appelée soluté, à se dissoudre dans une autre
substance, appelée solvant, pour former un mélange homogène appelé solution (Letcher,
2007 ; Pitzer, 2018a). La solubilité d’une substance dépend des propriétés physiques et
chimiques du soluté et du solvant, telles que la température, la pression et la présence
d’autres espèces chimiques (y compris les modifications du pH) dans la solution.
Le degré de solubilité d’une substance dans un solvant spécifique est mesuré en tant que
concentration de saturation, où l’ajout de plus de soluté n’augmente plus la concentration
de la solution et conduit à précipiter la quantité en excès du soluté. Le plus souvent, le
solvant est un liquide, qui peut être une substance pure ou un mélange.
Soit MνM XνX .ν0 H2 O un solide hydraté où νM est le nombre de cations M de charge ZM , νx
est le nombre d’anions X de charge Zx , ν0 est le nombre de molécules d’eau. En supposant
que ces coefficients sont obtenus à partir de la dissociation complète du solide, on a :
ln Ks = νM ln aM + νx ln ax + ν0 ln aw (1.42)
L’énergie libre molaire standard de dissolution de la réaction est déterminée par l’équa-
tion :
X
∆Gdis = ∆GSS Cr
f − ∆Gf (1.44)
17
Chapitre 1. Généralités et étude bibliographique
La majorité des modèles thermodynamiques pour les solutions électrolytiques sont basés
sur la théorie de Debye et Hückel (1923) . L’hypothèse principale de cette théorie
repose sur le fait que l’écart à l’idéalité est uniquement dû aux forces électrostatiques
entre les ions. Ceci est vrai pour les solutions très diluées, mais un peu moins pour les
solutions concentrées, surtout lorsque les forces d’interaction entre espèces deviennent
dominantes. Pour ces raisons, plusieurs modèles thermodynamiques ont été développés.
Les plus utilisés sont décrits ci-dessous.
L’équation de Debye–Hückel constitue une base de départ pour l’étude de solutions élec-
trolytiques non-idéales (Wright, 2007). Elle permet de prédire les coefficients d’activité
en fonction de la force ionique I et de la charge électrique zi comme suit :
log γi = −Azi2 I 1/2 (1.45)
avec :
1X
I= mi zi2 (1.46)
2 i
Debye et Hückel (1923) ont développé davantage cette équation pour décrire aussi les
solutions qui sont modérément concentrées. La nouvelle version est donnée par :
−Azi2 I 1/2
log γi = − (1.47)
1 + BaI 1/2
où a est le diamètre de l’ion (supposé sphérique) dont l’estimation de la valeur n’est pas
toujours aisée. Les paramètres A et B sont calculés comme suit :
3/2
e2
1 p
A= 2a ρw (1.48)
ln 10 4π0 kT
18
1.2. Modèles thermodynamiques
s
2Na ρw e2
B= (1.49)
0 kT
Cette équation présente un accord satisfaisant avec les mesures expérimentales pour les
faibles concentrations en espèces électrolytiques, généralement inférieures à 10−3 mol/L.
Des écarts par rapport à la théorie se produisent à des concentrations plus élevées et
avec des solutions qui contiennent des électrolytes qui ont des charges plus élevées, plus
particulièrement des électrolytes asymétriques.
Pour réduire ces écarts, Debye et Hückel (1923) ont encore amélioré leur équation en
ajoutant un paramètre ajustable C comme suit :
−Azi2 I 1/2
log γi = − + CI (1.50)
1 + BaI 1/2
D’autre part, Guggenheim (1935) propose une version pour le calcul du coefficient
d’activité moyen :
1/2
I
ln γ± = αz+ z− + 2νβm (1.51)
1 + I 1/2
avec :
2ν+ ν−
α = A ln 10 et ν = (1.52)
ν+ + ν−
où A correspond au paramètre défini par l’équation (1.48), ν est un paramètre qui s’ex-
prime en fonction des coefficients stœchiométriques et β est un paramètre d’interaction
binaire.
Dans (Davies, 1938), une nouvelle extension empirique de la loi de Debye-Hückel est
présentée. Elle s’exprime par :
I 1/2
log γ± = −2 − CI (1.53)
1 + I 1/2
Bromley (1973) propose finalement l’extension la plus flexible et qui peut être utilisée
pour de grandes valeurs de la force ionique. Elle a pour expression :
19
Chapitre 1. Généralités et étude bibliographique
B0 − B
= 0.06 + 0.6B (1.55)
|z+ z− |
Quant à la valeur de C elle est souvent prise égale à zéro. L’équation 1.54 devient donc :
log γi = fi + Li + Pi (1.58)
Avec :
−A | zi | I 1/2
fi = (1.59)
1 + I 1/2
NO 2
X | zi | + | zi |
Li = βij mj (1.60)
j=1
2
(0.06 + 0.6)Bij | zi zj |
βij = 2
+ Bij + Cij I + Dij I 2 (1.61)
1 + 1.5I | zi zj |)
20
1.2. Modèles thermodynamiques
NM
X zi2
Pi = Qij mj + Gj (1.62)
j
4I 2
NS
X
Gj = 0.86859mj Vjk mk (1.64)
k=1
(1)
Vjk = βjk [1 − (1 + 2I 1/2 + 2I)exp(−2I 1/2 )] (1.65)
Dans les équations 1.58-1.64, f est le terme de Debye-Hückel, Li sont les termes de
Bromley-Zemaitis pour prendre en considération les interactions ions-ions, N O est le
nombre d’ions de charges opposées à celle de l’électrolyte considéré, Bij , Cij et Dij sont
des paramètres dépendant de la température et qui permettent de prendre en compte
les interactions cation-anion. Pi sont les termes de Bromley-Zemaitis pour prendre en
compte les interactions ions-neutres. N M est le nombre de molécules en solution, N S
(0) (1)
est le nombre d’électrolytes en solution et βij et βij sont des paramètres d’interactions
qui dépendent de la température et de la pression. Les coefficients d’activité des espèces
neutres sont calculés comme suit :
NM
X
log γi = 2 Qij mi (1.66)
j=1
Bien que le modèle de Bromley-Zemaitis soit plus efficace que la théorie de Debye–Hückel,
et a été utilisé dans des travaux pour des conditions de température allant jusqu’à 100°C
et des molalités qui dépassent 20 mol/kgw (Li et Demopoulos, 2006), ses équations de
bases sont essentiellement empiriques et le rendent limité par rapport aux modèles à base
d’énergie libre de Gibbs. Son utilisation pour la prédiction des grandeurs thermodyna-
miques qui n’ont pas servi à son identification sont sujettes à beaucoup d’incertitudes. De
plus, la comparaison de ce modèle à d’autres montre qu’il présente toujours des limita-
tions en terme de précision (Foti et al., 1997). Les modèles les plus intéressants qui sont
basés sur l’énergie libre de Gibbs sont présentés dans les paragraphes suivants.
21
Chapitre 1. Généralités et étude bibliographique
Gex∗,pdh est l’expression de Pitzer Debye-Hückel et Gex∗,SR est exprimée comme suit :
∗
Gex ,SR Gex,SR X
xck ln γc∞k + xak ln γa∞k (1.69)
= −
RT RT k
avec :
Eek m Emek
ln γc∞k = Zck C − exp +1+ (1.70)
CRT CRT
Eek m Emek
ln γa∞k = Zak C − exp +1+ (1.71)
CRT CRT
GexSR
P
Xm + k (Xck + Xak ) Hek m X
= −CXm ln P −C Xck .
RT Xm + k (Xck + Xak ) k
(1.72)
" P !# " P !#
Xm Hek m + j Xaj Hej ek X Xm Hek m + j Xcj Hej ek
ln P −C Xak ln P
Xm + j Xaj k
Xm + j Xcj
X j = xj K j (1.73)
−Eij
Hij = exp (1.74)
CRT
xj est la fraction molaire de l’espèce j, Kj est égal à 1 pour les espèces moléculaires et égal
22
1.2. Modèles thermodynamiques
à zj pour les espèces chargées. Hij sont des paramètres ajustables exprimés en fonction
de Eij .
Eck m = Eak m = Eek m ; Emck ,ak ck = Emak ,ck ak = Em,ek ; Eal ck = Ecl ak = Eel ek (1.75)
Le modèle de Wilson donne une bonne représentation de l’enthalpie libre d’excès pour un
grand nombre de mélanges (Zhao et al., 2000). Elle se montre particulièrement efficace
pour des températures et des concentrations élevées (jusqu’à 120 °C et 25mol/kg d’eau
(Mazloumi et Shojaeian, 2019)). Toutefois, très peu de contributions dans la littérature
se sont intéressées à l’utilisation de ce modèle pour la modélisation de systèmes similaires
aux nôtres. Il serait alors intéressant de le comparer à d’autres modèles qui se basent sur
l’énergie libre de Gibbs. Parmi ces derniers, le modèle eUNIQUAC, décrit ci-dessous, est
à considérer.
En effet, le terme combinatoire prend en compte la différence de taille entre les particules.
Il est exprimé selon Sander et al. (1986) :
Gex
X φi zX φi
Combinatorial
= xi ln − qi xi ln (1.76)
RT i
xi 2 i θi
xi q i
θi = P (1.78)
l xl q l
23
Chapitre 1. Généralités et étude bibliographique
avec :
uki − uii
ψki = exp − (1.80)
T
Gex b2 I
Debye-Hückel 4Aφ 1/2 1/2
(1.81)
= −xw Mw 3 ln 1 + bI − bI +
RT b 2
ion
X
−1
Gex
LR = − (3DS ) si zi2 e2 τ (κa) (1.82)
i=1
24
1.2. Modèles thermodynamiques
Les constantes a1 et a2 ont été déterminées de manière empirique. Les meilleurs résultats
ont été obtenus avec a1 = −1 pour les interactions ions-ions, et a1 = −1.2 pour les inter-
actions ions-solvant. Dans les deux cas, a2 = 0.13.
Un dernier terme courte-portée, est pris en compte. Il exprime la contribution des interac-
tions des espèces non chargées ; il correspond aux interactions de courte portée dans une
solution non électrolytique et est représenté par le modèle UNIQUAC. Son expression est
donnée par :
Gex ex ex
SR = GCombinatorial + GResidual (1.86)
Ces deux termes sont ceux décrits par les équations (1.76) et (1.79), excepté que Li et al.
(1994) fixent les paramètres de volume et de surface du composé i, respectivement ri et
qi , à 1.
Récemment, Mohs et Gmehling (2013) ont présenté une revue bibliographique sur
le modèle LIQUAC. Leurs travaux ont montré l’efficacité de ce modèle pour plusieurs
systèmes électrolytiques et aussi pour des conditions de température et de concentration
élevées. La seule limite des deux modèle eUNIQUAC et LIQUAC est que très peu de
contributions dans la littérature se sont intéressées à leur utilisation pour la modélisa-
tion des systèmes électrolytiques qui nous intéressent. Ceci rend l’accès aux paramètres
inconnus de ces modèles très limité. Ce qui n’est pas le cas du modèle de Pitzer pour le-
quel les paramètres sont disponibles dans plusieurs bases de données thermodynamiques
pour un grand nombre de système électrolytiques. Le modèle de Pitzer est décrit dans le
paragraphe suivant.
Gex X X
= nw f + mi mj λi,j + mi mj ψi,j,k + ... (1.87)
Wn RT i,j i,j,k
25
Chapitre 1. Généralités et étude bibliographique
où λi,j représente les interactions binaires à courte distance entre deux espèces i et j, ψi,j,k
représente les paramètres d’interactions ternaires, T est la température, Wn est l’inverse
de la masse molaire du solvant, et f est la fonction de Debye et Hückel (1923) qui
dépend de la force ionique selon la relation suivante :
4Aφ I
ln 1 + I 0.5 (1.88)
f =−
b
avec :
1 1 3
Aφ = (2πN0 dw /1000) 2 (z 2 DkT ) 2 (1.89)
3
Les coefficients d’activité γi sont exprimés en fonction de la fonction des paramètres des
interactions binaires et ternaires, et de la fonction de Debye-Hückel :
zi2 X X 0
X
ln(γi ) = f +2 mj λi,j (I) + zi2 mj mk λj,k (I) + 3 mj mk ψi,j,k + ... (1.90)
2 j j j,k
où :
1X dλi,j ni
I= mi zi2 , λ0i,j = , mi = (1.91)
2 i dI nw
!
X
log(aw ) = Mw φ zi mi /1000 (1.93)
i
Pitzer donne des formes plus appropriées de sa théorie en utilisant des paramètres βi,j
0
pour exprimer les paramètres λi,j et λi,j en fonction de la force ionique :
(0) (1) 2
√ √
λi,j (I) = βi,j + βi,j 1 − (1 + α I)e−α I (1.94)
b
(1)
0 βi,j √ √
λi,j (I) = 2 −1 + (1 + α I + αI)e−α I
(1.95)
α I
26
1.2. Modèles thermodynamiques
De plus, ces paramètres dépendent de la température (Simoes et al., 2017) et sont géné-
ralement exprimés sous la forme suivante où T0 est égale à 298 K.
(0)
(0) (0) (0) ci,j
βi,j (T ) = ai,j + bi,j .(T − T0 ) + (1.96)
T − T0
(1)
(1) (1) (1) ci,j
βi,j (T ) = ai,j + bi,j .(T − T0 ) + (1.97)
T − T0
(2)
(2) (2) ci,j,k
ψi,j,k (T ) = ai,j,k + bi,j,k .(T − T0 ) + (1.98)
T − T0
Le modèle de Pitzer met en jeu des paramètres d’interactions binaires et ternaires, qui
dépendent du système thermodynamique étudié, et qui doivent être identifiés à partir de
mesures expérimentales, telles que la spéciation, l’activité de l’eau, le coefficient osmotique
etc.
où Gex
m est l’énergie de Gibbs d’excès du système électrolytique étudié, Gm
ex,lc
est la contri-
bution des interactions locales et Gm
ex,pdh
est la contribution qui découle des interactions
27
Chapitre 1. Généralités et étude bibliographique
ion-ion de longue distance. L’équation 1.99 permet d’exprimer le coefficient d’activité par :
Le modèle eNRTL permet de décrire l’énergie libre d’excès de Gibbs des interactions
locales. En effet, pour un système électrolytique, cette énergie peut s’écrire comme suit
(Chen et al., 2001) :
P P !
Gex,lc j Xj Gjm τjm j Xj Gjc,ac τjc,ac
X X X
m
= nm P + zc nc za P
RT m k X k Gkm c a k Xk Gkc,ac
P !
X X j X j Gja,ca τ ja,ca
+ za na zc P (1.101)
a c k Xk Gka,ca
avec :
X
nj = nI rj,I (1.102)
I
X j = C j xj (1.103)
nj
xj = P (1.104)
i ni
j = m, c, a (1.105)
Dans l’équation 1.101, le premier terme est la contribution lorsqu’une espèce moléculaire
se trouve au centre d’autres espèces (chargées ou neutres), le second est la contribution
lorsqu’une espèce cationique est au centre d’autres espèces (chargées ou neutres), et le troi-
sième terme est la contribution lorsqu’une espèce anionique se trouve au centre d’autres
espèces (chargées ou neutres) comme dans (Chen et al., 2001). Dans les équations (1.102-
1.105), i, j, k représentent différentes espèces, nI est le nombre de moles du constituant i
différent de j, Cj = zj , et Cj = 1 pour les espèces moléculaires. xj et Xj sont les fractions
molaires de l’espèce j.
28
1.3. Conclusion
Pour tenir compte des interactions ion-ion de longue distance, le modèle utilise la formule
asymétrique suivante (Pitzer, 2018a) :
1/2
Gex,P DH
1000 4AφI
= −N ln (1 + ρI 1/2 ) (1.107)
RT Ms ρ
avec :
1/2 3/2
Q2e
1 2πN0 ds
Aφ = (1.108)
3 1000 s kT
1X
I= xi zi2 (1.109)
2 i
Le modèle eNRTL est utilisé dans la littérature dans des conditions de température et
de concentration très élevées , allant jusqu’à 400 K pour la température (Honarparvar
et al., 2018) et jusqu’à 40 mol/kg d’eau pour la concentration (Wang et al., 2017). Il
serait par conséquent intéressant pour nos cas d’étude. Cependant, l’accès aux paramètres
inconnus de ce modèle le rend moins flexible par rapport au modèle de Pitzer.
1.3 Conclusion
L’analyse bibliographique a permis de définir l’ensemble des grandeurs thermodynamiques
nécessaires pour la modélisation des procédés, en se focalisant sur les méthodes qui per-
mettent de les estimer. Il s’est avéré que l’énergie libre de Gibbs est une grandeur fon-
damentale qui constitue la base pour la prédiction de ces grandeurs nécessaires. Elle est
souvent calculée à l’aide de modèles thermodynamiques.
Une revue bibliographique de ces modèles a ensuite été effectuée afin d’en déduire ceux qui
peuvent potentiellement être utilisés dans notre cas. Plusieurs modèles électrolytiques ont
ainsi été analysés et le modèle de Pitzer s’est révélé comme un des modèles les plus utilisés,
flexible et surtout fiable. Il est choisi comme modèle thermodynamique électrolytique dans
ce travail de thèse.
29
Chapitre 1. Généralités et étude bibliographique
30
Chapitre 2
Modélisation thermodynamique du
procédé de fabrication d’acide
phosphorique dihydrate
Sommaire
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2 Revue de la littérature . . . . . . . . . . . . . . . . . . . . . . . 33
2.3 Modélisation thermodynamique . . . . . . . . . . . . . . . . . . 34
2.3.1 Procédé dihydrate . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.3.2 Bilans de matière et de charge et équilibres . . . . . . . . . . . 35
2.4 Base de données expérimentales . . . . . . . . . . . . . . . . . . 39
2.4.1 Mesures effectuées en laboratoire dans ce travail . . . . . . . . . 39
2.4.2 Mesures de la littérature . . . . . . . . . . . . . . . . . . . . . . 39
2.5 Analyse d’estimabilité . . . . . . . . . . . . . . . . . . . . . . . . 44
2.5.1 Approche locale . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.5.2 Approche globale . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.6 Identification des paramètres . . . . . . . . . . . . . . . . . . . 47
2.7 Résultats et discussions . . . . . . . . . . . . . . . . . . . . . . . 55
2.7.1 Analyse de sensibilité globale . . . . . . . . . . . . . . . . . . . 55
2.7.2 Analyse d’estimabilité . . . . . . . . . . . . . . . . . . . . . . . 56
2.7.3 Identification des paramètres et seuil optimal d’estimabilité . . 58
2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
31
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
Résumé
Dans ce chapitre, un modèle est développé en se basant sur des bilans de matière et de
charge, des équations d’équilibres chimiques et du modèle thermodynamique de Pitzer.
Ce dernier met en jeu plusieurs paramètres d’interactions, binaires et ternaires, qu’il faut
identifier à partir de mesures expérimentales. Elle consiste en des données de spéciation
de sous-systèmes à base de soufre et de phosphate, des données de solubilité de dix mi-
néraux, et des données de l’activité de l’eau de huit systèmes binaires. Une partie des
mesures utilisées est réalisée en laboratoire, alors que le reste est collecté de la littéra-
ture. Par ailleurs, une des questions qui se posent souvent est de savoir si les mesures
disponibles permettent d’identifier tous les paramètres inconnus du modèle ou seulement
certains d’entre eux. Pour répondre à cette question, une approche d’estimabilité basée
sur les sensibilités globales des sorties du modèle par rapport aux paramètres inconnus est
développée. Cette approche est utilisée dans ce travail afin de déterminer les paramètres
les plus estimables pour chaque sous-système thermodynamique considéré. Les paramètres
les plus estimables sont alors identifiés, et leur précision est évaluée à l’aide de leur in-
tervalle de confiance. Les valeurs optimales des paramètres sont ensuite utilisées pour
comparer les prédictions du modèle aux mesures expérimentales. Les résultats montrent
un bon accord pour l’ensemble des sous-systèmes étudiés. Cet accord est confirmé par la
grande valeur du coefficient de corrélation de Pearson.
32
2.1. Introduction
2.1 Introduction
Les modèles thermodynamiques sont importants pour le développement des modèles de
procédés industriels. Ils permettent de prédire des grandeurs mises en jeu dans les équa-
tions des bilans de matière, de chaleur, d’équilibre de phase etc. Ces modèles peuvent
être simples ou très compliqués en fonction des espèces chimiques présentes dans les dif-
férentes phases (nombre, nature, et concentrations), et aussi des conditions opératoires
comme la température et la pression. Ces modèles thermodynamiques mettent souvent
en jeu de nombreux paramètres inconnus qui sont généralement identifiés à partir de
mesures expérimentales disponibles. La question qui se pose souvent est de savoir si les
mesures disponibles contiennent suffisamment d’informations pour permettre d’identifier
l’ensemble des paramètres inconnus ou seulement un sous-ensemble d’entre eux. Si l’on ne
répond pas correctement à cette question, l’identification peut conduire à des prédictions
imprécises des grandeurs thermodynamiques, et par conséquent à une dégradation de la
qualité des modèles des procédés.
33
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
systèmes à base de fluor dans les unités de production d’acide phosphorique, ainsi que la
possibilité de récupération du fluor à partir de ces systèmes. (Abutayeh et Campbell,
2009 ; Messnaoui et al., 2012) ont utilisé respectivement le modèle de Pitzer et le mo-
dèle NRTL-électrolyte pour l’étude et l’analyse des pertes chimiques du procédé d’acide
phosphorique dihydrate.
Au cours de cette réaction, l’acide sulfurique et l’acide phosphorique recyclé sont utili-
sés pour la dissolution et l’attaque du minerai de phosphate pour produire une bouillie
constituée d’un mélange d’acide phosphorique liquide et de sulfates de calcium solides.
34
2.3. Modélisation thermodynamique
Ces produits sont séparés dans une unité de filtration sous vide située en aval. L’acide
phosphorique est récupéré pour être valorisé et commercialisé, tandis que les sulfates de
calcium sont éliminés comme produit indésirable (Figure 2.1).
Le tableau 2.1 résume les réactions principales et secondaires qui sont dues à la spéciation
des réactifs utilisés lors de la digestion du minerai de phosphate par les acides sulfurique
et phosphorique. En effet, ces équilibres doivent être pris en compte par le modèle ther-
modynamique développé notamment dans les équations de bilans de matière, de charge
et d’équilibre. Ces équations sont décrites dans le paragraphe suivant.
NC
X
(M ) tot
= δM,i mi ; M ∈ {S, P, Ca, F, Si, M g, Al, F e, N a, Cl, K} (2.2)
i=1
35
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
NC
X
zi .mi = 0 (2.3)
i=1
NC NC
α
Y Y
Kj = ai ij = (mi .γi )αij , j = 1, .., N R (2.4)
i=1 i=1
∆Hj 1 1
log(Kj ) = log(Kj0 ) + − (2.5)
R T T0
À partir du système ci-dessus (Eqs. 2.2-2.5), les équations du modèle développé peuvent
être détaillées en se basant sur la spéciation des réactifs, sur les charges des différentes
espèces considérées en phase liquide, et sur les équilibres (tableau 2.1).
(P )tot = mH3 P O4 + mH2 P O4− + mHP O42− + mP O43− + 2.mH5 P2 O8− (2.7)
36
2.3. Modélisation thermodynamique
mHSO4− + 2.mSO42− + mH2 P O4− + 2mHP O42− + 3mP O43− + mH5 P2 O8− + mF − + mSiF6−
mH + mSO42− γH + .γSO42−
k1 = . (2.18)
mHSO4− γHSO4−
mH + mP O43− γH + γP O43−
k4 = . (2.21)
mHP O42− γHP O42−
mH + mF − γH + γF −
k6 = . (2.23)
mHF γHF
mHF2− γHF2−
k7 = . (2.24)
mHF mF − γHF γF −
37
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
m2H + mSiF62− γH
2
+ γSiF 2−
k8 = . 6
(2.25)
mH2 SiF6 γH2 SiF6
mCaSO4 γCaSO4
k9 = . (2.26)
mCa2+ mSO42− γCa2+ γSO42−
mN a+ mSO42− γN a+ γSO42−
k14 = . (2.31)
mN aSO4− γN aSO4−
m2K + mSO42− γK
2
+ γSO 2−
k15 = . 4
(2.32)
mKSO4− γKSO4−
mN a+ mCl− γN a+ γCl−
k16 = . (2.33)
mN aCl γN aCl
mK + mCl− γK + γCl−
k17 = . (2.34)
mKCl γKCl
38
2.4. Base de données expérimentales
Des mesures concernant l’activité de l’eau des systèmes à base de fluor et de silice ont
été récemment effectués. Pour réaliser ces mesures, des solutions d’acide fluorhydrique et
d’acide hexafluorosilicique ont été préparées. En effet, pour ces deux systèmes, l’activité
de l’eau a été mesurée par la méthode hygrométrique développée par (El Guendouzis
et Dinane, 2000). Elle est basée sur la mesure de l’activité de l’eau dans des solutions
aqueuses contenant des électrolytes non volatils. Pour toutes les mesures, la température
et la concentration des échantillons varient respectivement de 298 K à 353 K et de 0 à
6 moles/kgw pour le cas de l’acide fluorhydrique, et de 0 à 3 moles/kgw pour le cas de
l’acide hexafluorosilicique.
39
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
Tableau 2.1 – Réactions d’équilibre en phase aqueuse avec leur constante d’équilibre et
enthalpie à 298 K.
j Ki0 ∆H (J/mol) Références
1 K
1 HSO4 – )
−−
−*− SO4 2 – + H+ 0.0103 -1.274.104 Sippola et Taskinen (2014)
2 K
2 −−
H3 PO4 )−*− H2 PO4 – + H+ 0.0071 -4.271.103 Cherif et al. (2000)
3 K
3 H2 PO4 – )
−−
−*− HPO4 2 – + H+ 0.2550 -9.266.103 Cherif et al. (2000)
4 K
4 HPO4 2 – )
−−
−*− PO4 3 – + H+ 6.3.10−8 -1.800.104 Cherif et al. (2000)
5 K
5 H2 PO4 – + H3 PO4 )
−−
−*− H5 P2 O8 – 4.2.10−12 -1.500.104 Cherif et al. (2000)
K6
6 HF −
)−
−− F – + H+
* 7.2.10−4 -1.329.104 Elguendouzi et al. (2019)
7 K
7 HF + F – )
−−
−*− HF2 – 5.5000 -1.740.104 Elguendouzi et al. (2019)
8K
8 H2 SiF6 −
)−
−*− SiF6 2 – + 2 H+ 0.3.102 -6.799.104 Elguendouzi et al. (2015)
9 K
9 Ca2+ + SO4 2 – −
)−
−*− CaSO4 0.6394 -7.200.103 Shenn et al. (2020)
10 K
10 MgSO4 −
−
)−−− Mg2+ + SO4 2 –
*
− 5.623.10−3 -1.179.104 Gustafsson (2011)
K
11
11 MgCl2 −
−
)−−− Mg2+ + 2 Cl –
*
− 0.350 -1.729103 Gustafsson (2011)
K
12
12 AlCl3 + −
−
)−−− Al3+ + 3 Cl –
*
− 9.549.10−4 -8.995.103 J. Wang et al. (2016)
K13
13 FeSO4 −
−
)−−
*
−− Fe2+ + SO4 2 – 5.623.10−3 1.351.104 Kobylin et al. (2012)
K14
– −
14 NaSO4 )−−− Na+ + SO4 2 –
−−* 0.1995 -4.686.103 Christov (2000)
K15
– −
15 KSO4 )−−− 2 K+ + SO4 2 –
−−* 0.1412 -9.414.103 Christov (2000)
K16
16 −−
NaCl )−−
−− Na+ + Cl –
* 0.381.102 6.568.103 Christov (2000)
K17
17 −−
KCl )−−
−− K+ + Cl –
* 0.079.102 3.840.103 Christov (2000)
K18
18 −−
CaCl2 )−−
−− Ca2+ + 2 Cl –
* 0.201 2.020.103 Christov (2000)
40
2.4. Base de données expérimentales
41
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
42
2.4. Base de données expérimentales
43
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
Parmi les méthodes d’estimabilité paramétrique développées dans la littérature, les mé-
thodes basées sur l’orthogonalisation se sont avérées particulièrement pertinentes pour
classer les paramètres du plus estimable au moins estimable. La méthode utilisée dans
ce travail est basée sur l’utilisation d’une matrice de coefficients de sensibilité Z qui est
calculée à partir des coefficients de sensibilité individuels Sij (Lei et al., 2013 ; Ngo et al.,
2014 ; Yao et al., 2003) :
S1,1 |m1 ,T1 ··· S1,d |m1 ,T1
.. .. ..
Z= . . . (2.36)
Sn,1 |mn ,Tn · · · Sn,d |mn ,Tn
∆Yi
Si,j ≈ , i ∈ [1; n], j ∈ [1; d] (2.37)
∆Pj
où Yi est la sortie i du modèle et Pj est le paramètre inconnu j. Les éléments Si,j sont
ensuite multipliés par des coefficient de pondération, à savoir Pj et Yi , afin de garantir le
même ordre de grandeur des éléments Si,j ∗
:
Pj
∗
Si,j = Si,j (2.38)
Yi
La matrice Z est ensuite implémentée dans l’algorithme d’estimabilité développé dans
(Yao et al., 2003) afin de classer les paramètres du plus estimable au moins estimable. De
plus, l’implémentation de cette méthode nécessite la définition d’un seuil d’estimabilité
dont la valeur est choisie presque arbitrairement. Le nombre de paramètres les plus esti-
mables dépend donc de la valeur choisie du seuil et rend ainsi la méthode moins robuste.
44
2.5. Analyse d’estimabilité
Pour éviter ce choix arbitraire, Wu et al. (2011) ont proposé une méthode basée sur
le calcul des erreurs quadratiques moyennes. L’algorithme de cette méthode est présenté
ci-dessous.
— (i) Classer les paramètres du plus estimable au moins estimable en utilisant l’algo-
rithme d’orthogonalisation (Yao et al., 2003).
— (ii) Identifier le premier paramètre estimable en minimisant les erreurs quadratiques
J mse entre les mesures expérimentales et les prédictions du modèle. Ensuite, iden-
tifier les deux premiers paramètres estimables, puis les trois premiers et ainsi de
suite jusqu’à ce que tous les paramètres estimables soient identifiés. Les valeurs des
paramètres les moins estimables sont fixées à des valeurs nominales. Les erreurs
quadratiques moyennes J mse sont calculées comme suit :
n X k exp
!2
mod
X y i,j − yi,j
min J mse =
i=1 j=1
ωi
où yi,j
exp
et yi,j
mod
sont respectivement les mesures expérimentales et les sorties du
modèle, et ωi sont des facteurs de pondération.
— (iii) Noter la valeur de la fonction objectif pour l’identification des L paramètres les
plus estimables, tout en fixant les d − L paramètres restants constants.
— Calculer le rapport critique rc,L pour L = 1, ..., d comme suit :
JLmse − Jdmse
rc,L =
d−L
où JLmse définit la fonction objectif pour l’identification des L paramètres les plus
estimables.
— (iv) Calculer le rapport critique corrigé rcc,L selon :
d−L
rcc,L = (rcKub,L − 1)
nk
avec :
2
rcKub,L = max rc,L − 1, rc,L
d−L+2
rcc,L est utilisé pour comparer plusieurs modèles simplifiés avec différents nombres
de paramètres.
— (v) En déduire la valeur de L qui minimise la valeur de rcc,L . L correspond au nombre
de paramètres estimables à identifier à partir des mesures disponibles.
Bien que l’algorithme permette d’éviter le choix arbitraire de la valeur du seuil d’estima-
bilité, le classement des paramètres dépend toujours de la méthode de calcul de la matrice
Z. Cette dernière change avec les valeurs initiales des paramètres utilisés pour calculer
ses éléments, c’est-à-dire les sensibilités locales. Pour éviter ce problème, nous proposons
de remplir la matrice Z avec des sensibilités globales calculées selon l’approche de Sobol
(Sobol, 1990). C’est l’approche globale ; elle est détaillée ci-dessous.
45
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
E (V (Y |X∼i ))
STi = 1 − (2.39)
V(Y )
où V et E sont respectivement la variance et l’espérance mathématique de la sortie Y .
X∼i désigne toutes les entrées du modèle sauf Xi . Cependant, presque tous les praticiens
utilisent l’algorithme suivant pour estimer les indices de sensibilité totale de Sobol :
ˆ2
ˆ i = 1 − U∼ i − f0
ST (2.40)
V̂
où :
N
1 X
U∼ i = YA,k .YCi ,k (2.41)
N k=1
N
ˆ 1 X
f0 = YA,k (2.42)
N k=1
N
1 X 2 2
V̂ = YA,k − fˆ0 (2.43)
N k=1
avec :
(A) (A) (A)
YA,k = f (Xk,1 , ..., Xk,i , ..., Xk,d ) (2.44)
(A) (B)
Les vecteurs d’entrée Xi,j et Xi,j sont générés au moyen d’une méthode d’échantillonnage
de Monte-Carlo et utilisés comme éléments des matrices suivantes A et B comme suit :
46
2.6. Identification des paramètres
Ces matrices sont utilisées à leur tour pour générer les matrices Ci qui sont spécifiques à
chaque paramètre inconnu i comme suit :
k correspond aux grandeurs mesurées, c’est-à-dire l’activité de l’eau, les molalités des
composants et la solubilité, ωk sont des facteurs de pondération, T et m font référence à
la température et aux concentrations totales. La précision de l’identification est évaluée au
moyen des intervalles de confiance (Langman, 1986). En effet, on suppose que la fonction
objectif peut être exprimée sous la forme suivante :
47
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
où θ(P ) est le vecteur résiduel, c’est-à-dire la différence entre les prédictions du modèle
et les valeurs mesurées des sorties, W est une matrice de pondération (égale à la matrice
identité, puisque nous supposons que toutes les mesures ont le même poids dans la fonction
objectif), et P est le vecteur des paramètres inconnus. En supposant également que les
erreurs de mesure sont indépendantes et normalement distribuées, la matrice de covariance
C du problème des moindres carrés est approximée comme suit (Langman, 1986) :
F (P ∗ ) T −1
C≈ (J J) (2.51)
d−n
P ∗ est le vecteur des paramètres qui minimise la fonction objectif F (P ), n et d sont res-
pectivement le nombre d’expériences et le nombre de paramètres inconnus, et J est la
matrice jacobienne du vecteur θ(P ) (Eq. 2.52). Cette approximation est acceptable tant
que les non-linéarités ne sont pas fortes.
La matrice jacobienne est calculée en utilisant une méthode locale "One-At-Time" (OAT).
Elle consiste à perturber la valeur de chaque paramètre Pj∗ de 10% à gauche et à droite,
et la méthode des différences finies centrées est utilisée pour approximer les éléments de
J. Il convient de noter qu’une perturbation de 10% est choisie d’après les nombreux tra-
vaux antérieurs qui ont abordé les mêmes types de problèmes (Benyahia et al., 2009 ;
Benyahia, 2009 ; Benyahia et al., 2010 ; Benyahia et al., 2011 ; Benyahia et al., 2013 ;
Benyahia, 2018 ; Benyahia et al., 2019 ; Benyahia et al., 2021b ; Bouchkira, 2022 ;
Cardenas et al., 2020 ; Cardenas et al., 2021 ; Cristian et al., 2022 ; Fysikopoulos
et al., 2019 ; Lei et al., 2013 ; Ngo et al., 2014 ; Onyemelukwe et al., 2018 ; Renom
Carrasco et al., 2020).
Par ailleurs, J peut être utilisée pour apprécier si la réalisation d’expériences supplé-
mentaires peut augmenter le nombre de paramètres estimables (analyse d’estimabilité),
tandis que la matrice d’information de Fisher (F IM ) (exprimée par F IM = J T V −1 J,
où V est la matrice de variance-covariance des résultats) peut être utilisée pour définir
ces expériences (plan d’expériences optimal). En outre, la F IM est souvent utilisée pour
déterminer les intervalles/régions de confiance des paramètres.
∂θ1 ∂θ1
∂P1
··· ∂Pd
J = ... ... ..
. (2.52)
∂θn ∂θn
∂P1
··· ∂Pd
L’incertitude sur un paramètre j est déterminée au moyen de l’équation 2.53, où cjj est
l’élément diagonal j de la matrice de covariance C. Cette dernière est généralement es-
timée à partir de mesures expérimentales, puis utilisée pour estimer les intervalles de
48
2.6. Identification des paramètres
confiance.
Cependant, puisque dans ce travail l’estimation de C n’a pas été facile à réaliser, elle a été
approximée en utilisant F (P ∗ ) qui représente les résidus à l’optimum entre les prédictions
du modèle et les expériences utilisées. De plus, la division par P ∗ nous permet de quantifier
la précision de l’identification d’un paramètre P :
√
cjj t1−α/2,ν
Pj = ± .100% (2.53)
Pj∗
√ √
Pj ∈ Pj∗ − cjj t1−α/2,ν ; Pj∗ + cjj t1−α/2,ν (2.54)
D’autre part, le coefficient de corrélation de Pearson est calculé pour évaluer la précision et
la fiabilité des prédictions du modèle thermodynamique développé selon (Keith, 2014) :
Sijm − Sijm kije − kije
P
i,j
r = r 2 P 2 (2.55)
m m e e
P
i,j (Sij − Sij ) (k
i,j ij − k ij )
où Sijm est la prédiction du modèle, kije est la valeur expérimentale correspondante, Sijm et
kije sont respectivement leur moyenne sur le nombre total de mesures disponibles effectuées
à différentes molalités et températures. Les indices i et j font référence aux molalités et
aux températures respectivement.
49
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
50
2.6. Identification des paramètres
51
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
Figure 2.4 – Prédictions du modèle par rapport aux mesures de la littérature pour les
données de solubilité de huit minéraux différents.
52
2.6. Identification des paramètres
Figure 2.5 – Comparaison des prédictions et des mesures de l’activité de l’eau dans
quatre systèmes aqueux binaires.
53
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
Figure 2.6 – Comparaison des prédictions et des mesures de l’activité de l’eau dans
quatre systèmes aqueux binaires.
54
2.7. Résultats et discussions
Les résultats de l’analyse de sensibilité sont utilisés pour deux objectifs principaux. Le pre-
mier objectif est d’évaluer la sensibilité des sorties du modèle par rapport aux paramètres
considérés. Pour ce faire, les sensibilités globales des 21 paramètres d’interactions sont cal-
culées et comparées à un critère proposé par Zhang et al. (2015) qui supposent qu’un
paramètre dont la sensibilité globale est inférieure à 0.05 n’a pas une grande influence sur
les sorties du modèle. Le deuxième objectif est d’utiliser ces indices de sensibilité globale
calculés pour remplir la matrice Z.
Les figures 2.7A-D montrent les indices de sensibilité globale des sorties mesurées du mo-
dèle (c’est-à-dire l’activité de l’eau et les molalités des ions) par rapport aux 21 paramètres
inconnus. Les sensibilités sont calculées en utilisant un échantillonnage de Monte-Carlo de
10 000 vecteurs de taille 21 (nombre de paramètres). En tenant compte de l’algorithme de
calcul de sensibilité décrit précédemment, 230 000 simulations sont effectuées. Le temps
CP U nécessaire pour effectuer toutes les simulations a été estimé à 64 heures en utilisant
une station de calcul Dell Precision T 7810 Bi-Xeon 12 x Core 64GB. On peut remarquer
que les sorties du modèle sont sensibles à la variation de l’ensemble des paramètres du
modèle (Figs. 2.7A-D). L’élimination arbitraire d’un ou de plusieurs de ces paramètres
causerait sans doute de l’incertitude sur les sorties du modèle. L’analyse d’estimabilité
nous permettra donc de savoir quels paramètres sont estimables à partir des données dis-
ponibles et ceux dont les valeurs doivent être fixées à partir de la littérature ou de travaux
antérieurs.
55
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
La figure 2.8B, présente les résultats de l’estimabilité des 21 paramètres inconnus sur la
base de 168 valeurs mesurées de molalités des différents ions(c’est-à-dire H + , HSO4− , et
SO42− ) à quatre températures différentes et 14 concentrations totales d’acide sulfurique.
Dans ce cas, la matrice de sensibilité Z est composée de 21 colonnes et de 168 mesures.
L’analyse de la figure 2.8B montre que les magnitudes des paramètres inconnus sont plus
élevées par rapport aux résultats du premier cas, ceci montre que ces mesures contiennent
plus d’informations pour identifier plus de paramètres inconnus.
La figure 2.8C présente le cas où les mesures de l’activité de l’eau et les molalités des
espèces sont simultanément prises en compte. Dans ce cas, la matrice de sensibilité Z
comprend 21 colonnes et 224 points de mesures. Les magnitudes de tous les paramètres
sont plus élevées par rapport aux cas précédents. Ceci est cohérent avec le fait qu’en
ajoutant plus de mesures, le nombre de paramètres estimables augmente en raison de
l’augmentation d’information que ces mesures contiennent.
56
2.7. Résultats et discussions
Figure 2.7 – Indices de sensibilité globale de Sobol des 21 paramètres du modèle calculés
à partir d’un ensemble de 10 000 échantillons. (A), (B), (C) et (D) sont les sensibilités
globales par rapport à la molalité HSO4− , à la molalité SO42− , à l’activité de l’eau et à la
molalité H + respectivement.
57
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
(0)
Paramètres βi,j
(0) (0) (0)
1 aHSO− ,SO2− 8 bHSO− ,SO2− 15 cHSO− ,SO2−
4 4 4 4 4 4
(0) (0) (0)
2 aHSO− ,H + 9 bHSO− ,H + 16 cHSO− ,H +
4 4 4
(0) (0) (0)
3 aSO2− ,H + 10 bSO2− ,H + 17 cSO2− ,H +
4 4 4
(1)
Paramètres βi,j
(1) (1) (1)
4 aHSO− ,SO2− 11 bHSO− ,SO2− 18 cHSO− ,SO2−
4 4 4 4 4 4
(1) (1) (1)
5 aHSO− ,H + 12 bHSO− ,H + 19 cHSO− ,H +
4 4 4
(1) (1) (1)
6 aSO2− ,H + 13 bSO2− ,H + 20 cSO2− ,H +
4 4 4
Paramètres ψi,j,k
(2) (2) (2)
7 aHSO− ,SO2− ,H + 14 bHSO− ,SO2− ,H + 21 cHSO− ,SO2− ,H +
4 4 4 4 4 4
58
2.7. Résultats et discussions
Figure 2.9 – Valeurs de J mse pour chaque cas mse. (A) : activité de l’eau, (B) : molalités
des espèces, (C) : activité de l’eau et molalités des espèces.
Figure 2.10 – Valeurs de Jcc,L pour chaque cas mse. (A) : activité de l’eau, (B) : molalités
des espèces, (C) : activité de l’eau et molalités des espèces.
La figure 2.10A montre que 7 paramètres sont estimables lorsque seule la mesure de
l’activité de l’eau est utilisée. Le seuil d’estimabilité correspondant est d’environ 0.03
(l’ordonnée du paramètre 7 sur la figure 2.9A), ce qui signifie que pour qu’un paramètre
soit estimable, une variation de 10 % de sa valeur doit entraîner une variation d’au moins
1.75 % des sorties du modèle. Lorsque seules les mesures des molalités des espèces sont
utilisées, le nombre de paramètres les plus estimables ainsi que la valeur correspondante
du seuil d’estimabilité sont plus élevés que dans le cas précédent, et sont respectivement
égaux à 11 et 0.07 (Fig. 2.10B). Cela signifie que les molalités des espèces contiennent plus
d’informations que l’activité de l’eau. La variation minimale des sorties du modèle causée
par 10 % d’un paramètre estimable est d’environ 2.65 %. Enfin, comme prévu, lorsque
les deux ensembles de mesures, c’est-à-dire l’activité de l’eau et les molalités des espèces
sont utilisées simultanément, le nombre de paramètres les plus estimables est encore plus
élevé et est égal à 17 (Fig. 2.10C). La valeur du seuil d’estimabilité est d’environ 0.05,
59
Chapitre 2. Modélisation thermodynamique du procédé de fabrication d’acide phosphorique dihydrate
c’est-à-dire que la variation minimale des sorties du modèle causée par une variation de
10 % doit alors être supérieure ou égale à 2.24 %.
La même méthodologie est appliquée aux autres systèmes et les résultats sont résumés
dans le tableau 2.6. Il présente l’ensemble des systèmes étudiés, les références aux équi-
libres considérés pour chaque cas dans le tableau 2.1, le nombre de paramètres binaires et
ternaires, le seuil d’estimabilité déterminé pour chaque système, le nombre de paramètres
estimables à partir des mesures disponibles, le nombre de paramètres fixés de travaux
antérieurs ou de la littérature, et finalement les références des figures qui comparent les
prédictions du modèle et les mesures expérimentales. L’activité de l’eau et la spéciation ont
été utilisées comme mesures disponibles pour les systèmes H2 SO4 −H2 O et H3 P O4 −H2 O,
l’activité de l’eau a servi comme mesure disponible pour l’étude des systèmes binaires
HF −H2 O, H2 SiF6 −H2 O, M gSO4 −H2 O, M gCl2 −H2 O, AlCl3 −H2 O et CaCl2 −H2 O,
alors que la solubilité de plusieurs minéraux a été utilisée pour l’étude des systèmes
ternaires CaSO4 − H2 SO4 − H2 O, M gSO4 − M gCl2 − H2 O, F eSO4 − H2 SO4 − H2 O,
N aSO4 − N aCl − H2 O et K2 SO4 − M gSO4 − H2 O.
Les résultats montrent que des mesures contiennent plus d’informations nécessaires pour
l’identification des paramètres que d’autres. Notamment, la combinaison de l’activité de
l’eau et de la spéciation permet d’identifier plus de paramètres, ce qui est l’exemple des cas
des systèmes binaires d’acide sulfurique et d’acide phosphorique. Par ailleurs, l’activité de
l’eau seule comme mesure ne s’est pas avérée très efficace pour l’identification de certains
paramètres dans des systèmes binaires comme le cas des sulfates de magnésium (seuls 2
paramètres estimables), cas des chlorures de magnésium (1 seul paramètre estimable), cas
des chlorures d’aluminium (seuls 2 paramètres estimables). Finalement, la solubilité de cer-
tains minéraux a permis d’identifier un nombre modeste de paramètres. En effet, certains
paramètres des systèmes ternaires ont été identifiés avec les mesures des systèmes binaires
et ont été fixés. Par exemple, pour le cas du système ternaire CaSO4 − H2 SO4 − H2 O, les
βH + ,SO42− ont été fixés à partir de l’identification des paramètres du système H2 SO4 −H2 O.
Les valeurs de ces paramètres auraient pu être identifiées pour les systèmes ternaires, mais
leur identification avec les mesures des systèmes binaires s’est avérée plus satisfaisante,
vu que pour les systèmes binaires, le rapport du nombre de mesures et du nombre de
paramètres inconnus est plus important. Ceci permet d’une part une identification plus
précise, et d’autre part, minimise le temps de calcule pour l’identification dans les systèmes
ternaires.
2.8 Conclusion
Ce chapitre a pour but le développement d’un modèle thermodynamique nécessaire pour la
prédiction des propriétés à utiliser pour la simulation et l’optimisation du procédé d’acide
phosphorique dihydrate. Les équations des bilans de matière, de charge, des équilibres
et du modèle de Pitzer forment le modèle développé. L’utilisation du modèle de Pitzer
nécessite l’identification de ses paramètres à partir des mesures disponibles. A cette fin,
une revue concernant les données disponibles est effectuée et a abouti à une base de don-
60
2.8. Conclusion
nées contenant les mesures de solubilité de dix minéraux et l’activité de l’eau dans huit
systèmes aqueux binaires mis en jeu dans le procédé étudié. En outre, des mesures de
spéciation ont été effectuées pour l’acide sulfurique et l’acide phosphorique (réactifs du
procédé de production).
Pour chaque sous-système thermodynamique considéré, une analyse d’estimabilité est réa-
lisée afin de déterminer les paramètres les plus estimables. De plus, à chaque paramètre
identifié est associé un intervalle de confiance, et un coefficient de corrélation est calculé
pour évaluer la fiabilité du modèle développé dans la prédiction de plusieurs propriétés
thermodynamiques, telles que l’activité de l’eau, la spéciation et les solubilités.
Bien que les résultats obtenus soient très prometteurs, le modèle mérite d’être davantage
amélioré. En effet, plusieurs paramètres inconnus ne sont pas estimables à partir des don-
nées disponibles et leur valeur sont fixées à partir de travaux antérieurs. Pour répondre
à cette préoccupation, plus de mesures d’équilibre pouvant être utilisées sont en cours de
réalisation et seront prises en compte dans nos futurs travaux.
61
62
Tableau 2.6 – Résultats de l’analyse d’estimabilité et d’identification pour l’ensemble des systèmes étudiés
Système Equilibres βi,j ψi,j,k mesures Seuil d’estimabilité Pj estimables Pj fixés Comparaisons
Exploitation du modèle
thermodynamique pour l’optimisation
multicritère du procédé dihydrate
Sommaire
3.1 Contexte de l’étude . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.2 Revue de la littérature . . . . . . . . . . . . . . . . . . . . . . . 67
3.3 Problème d’optimisation . . . . . . . . . . . . . . . . . . . . . . 67
3.3.1 Fonctions objectifs . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.3.2 Variables de décision . . . . . . . . . . . . . . . . . . . . . . . . 69
3.3.3 Contraintes d’égalité et d’inégalité . . . . . . . . . . . . . . . . 70
3.4 Résultats de l’optimisation . . . . . . . . . . . . . . . . . . . . . 73
3.4.1 Approche de résolution . . . . . . . . . . . . . . . . . . . . . . . 73
3.4.2 Front de Pareto . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.4.3 Méthode d’aide à la décision . . . . . . . . . . . . . . . . . . . . 75
3.4.4 Modèle de substitution . . . . . . . . . . . . . . . . . . . . . . . 77
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
63
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
Résumé
64
3.1. Contexte de l’étude
La digestion s’effectue dans une cuve cylindrique constituée d’une série de neuf réacteurs
continus, de même volume et uniformément répartis à l’intérieur de la cuve. Cette dernière
est divisée en trois sections de trois réacteurs continus chacune. Le débit d’alimentation en
acide sulfurique est donc divisé en trois flux, chacun alimentant une des trois sections. Le
minerai de phosphate est introduit dans le premier réacteur (R1) de la première section et
le produit final quitte la cuve au niveau du dernier réacteur (R9) de la troisième section.
Une partie de la pulpe est recyclée depuis le dernier réacteur (R9) vers le premier réac-
teur (R1) en passant par un Flash permettant de contrôler la température de la réaction.
(Figure 3.1).
65
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
Au cours de la réaction de digestion, comme le montre la réaction globale (Eq. 3.1), l’acide
sulfurique est dissocié en ions H + et SO42− . Les ions H + sont mis en jeu dans l’attaque
des éléments phosphates et leur dissolution dans la phase liquide (solution d’acide phos-
phorique), tandis que les ions SO42− cristallisent avec les ions calcium Ca2+ pour former
un sulfate de calcium solide. Ce dernier peut être soit de l’anhydrite CaSO4 , du gypse
CaSO4 : 2H2 O ou de la bassanite CaSO4 :0.5H2 O, dépendant des conditions de fonction-
nement, notamment de la température et de la concentration en sulfates (Becker, 1989 ;
Dorozhkin, 1996).
L’acide phosphorique et les sulfates de calcium sont ensuite séparés dans une unité de
filtration sous vide située en aval de la cuve. L’acide phosphorique est récupéré pour être
valorisé et commercialisé, tandis que les sulfates de calcium solides sont éliminés comme
produit indésirable.
Par ailleurs, la quantité d’acide sulfurique requise doit être utilisée de manière optimale
afin d’optimiser les performances des unités de production. En effet, un déficit d’acide
sulfurique lors de la digestion entraîne une diminution de la concentration en ions sulfate
dans la cuve. Par conséquent, les ions calcium libres (provenant du minerai) commencent
à capturer les ions phosphatés HP O42− et à les cristalliser sous forme de brushite solide
CaHP O4 :2H2 O. Ceci diminue le rendement chimique du procédé puisque la brushite for-
mée est évacuée avec la phase solide comme produit indésirable, et entraîne des pertes de
phosphate par syncristallisation (Becker, 1989 ; Dorozhkin, 1996). Par contre, l’uti-
lisation d’un grand excès d’acide sulfurique provoque non seulement une augmentation
locale de la température en raison de la chaleur de dilution, mais aussi à la formation
d’une couche de gypse autour des particules de phosphate (enrobage des particules), les
protégeant ainsi de toute attaque ultérieure par l’acide sulfurique et entraînant également
des pertes de phosphate inattaqué. La quantité appropriée d’acide sulfurique ainsi que
sa répartition optimale sur les trois sections de la cuve d’attaque sont donc importantes
pour le fonctionnement optimal du procédé.
66
3.2. Revue de la littérature
En effet, Gioia et al. (1977) ont abordé dans leur travail la modélisation du procédé
hémihydrate (attaque et filtration). Ils se sont intéressés au dimensionnement de différents
procédés pouvant être plus économiques en terme de consommation de réactifs et ainsi
plus rentables en terme de production. Papadopoulos et Seferlis (2009) ont égale-
ment traité l’optimisation du même procédé (hémihydrate) en déterminant les conditions
opératoires optimales qui maximisent la quantité d’acide produit. Plus tard, Joao et al.
(2016) se sont intéressés à la modélisation des unités de concentration d’acide phospho-
rique, et ont utilisé des outils de simulation dynamique pour l’analyse et l’optimisation de
ces unités. Récemment, Ma et al. (2018) ont développé et optimisé des diagrammes de
flux pour la production d’acide phosphorique dihydrate. Toutefois, leur objectif principal
était de minimiser la quantité d’eau douce introduite dans le procédé. L’optimisation mul-
ticritère a été utilisée par Hemalatha et Rani (2017) pour déterminer les paramètres
de marche optimale dans un réacteur de refroidissement exploité lors de la production
d’acide phosphorique par les procédés dihydrate. Très récemment, (Bouchkira et al.,
2021b ; Bouchkira et al., 2021c ; Bouchkira et al., 2022) ont développé des approches
de modélisation et optimisation multicritère de la cuve de digestion d’un procédé industriel
de fabrication d’acide phosphorique dihydrate.
Les deux premiers objectifs sont pris en compte à l’aide de deux critères dans le problème
d’optimisation, tandis que les deux autres objectifs sont considérés comme des contraintes.
67
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
En effet, les deux premiers objectifs peuvent être prédits directement par la modélisation
thermodynamique, ce qui n’est pas le cas pour les deux derniers. Par exemple, le modèle
thermodynamique développé ne peut pas être utilisé pour prédire la quantité de parti-
cules qui subiront le phénomène d’enrobage qui les protégeront de l’attaque ultérieure
de l’acide sulfurique (objectif 3). Le modèle ne peut non plus être utilisé pour prédire
la filtrabilité du gâteau du gypse produit (objectif 4). En effet, des travaux expérimen-
taux antérieurs dans la littérature nous ont suggérés de considérer ces deux objectifs sous
forme de contraintes dans le problème. Il s’agit notamment des travaux expérimentaux
de Becker (1989) qui ont montré que l’excès optimal d’acide sulfurique devait se si-
tuer entre 1% et 3%. Il a notamment montré qu’en opérant au-dessus de 3% d’excès
d’acide sulfurique, la nucléation du gypse à la surface des particules de phosphate était
très favorable et conduisait à la formation d’une couche de gypse autour d’elles, ce qui
limitait leur contact avec l’acide sulfurique. De plus, de nombreux travaux de la littéra-
ture ont montré qu’un excès de 1% d’acide sulfurique pouvait être considéré comme la
quantité minimale pour provoquer une sursaturation du milieu réactionnel et ensuite pré-
cipiter le gypse (Becker, 1989 ; Dorozhkin, 1996 ; Peng et al., 2015 ; Zhu et al., 2016).
Enfin, une dernière inégalité est utilisée pour tenir compte du dernier objectif (4), et sug-
gère d’opérer en dessous de 355 K pour éviter la cristallisation de la bassanite, qui est
caractérisée par une viscosité élevée et des particules de petite taille, diminuant ainsi la
filtrabilité du gypse.
Les deux critères du problème d’optimisation sont exprimés en terme d’indices de satura-
68
3.3. Problème d’optimisation
N
!
1 X Y c
f1 = log akk .Ks−1 (3.2)
N i=1 k
1
i
N
!
1 X Y c
f2 = log akk .Ks−1 (3.3)
N i=1 k
2
i
ak est l’ activité de l’ion k mis en jeu. ks1 et ks2 sont les produits de solubilité de la
brushite et du gypse respectivement. ck est le coefficient stœchiométrique (négatif pour
les réactifs et positif pour les produits). i désigne une section de la cuve de digestion et
N est le nombre de sections.
Comme rapporté par Becker (1983) , dans les trois sections, la brushite est formée
par l’association des ions Ca2+ provenant de la roche phosphatée et des ions phosphates
HP O42− provenant de l’acide phosphorique produit (Eq. 3.4). De la même manière, le gypse
est formé par l’association des ions calcium Ca2+ avec les ions sulfate SO42− provenant de
la dissociation de l’acide sulfurique en milieu aqueux (Eq. 3.5).
2H O
(3.4)
2
Ca2+ + HPO4 2− )
−−
−−
−*− CaHPO4 · 2 H2 O
2H O
(3.5)
2
Ca2+ + SO4 2− )
−−
−−
−*− CaSO4 · 2 H2 O
Les fonctions objectifs peuvent finalement être exprimées en utilisant les activités des
espèces mises en jeu dans les réactions comme suit :
1X
f1 = log aCa2+ .aHP O42− .aH2 O Ks−1
2
(3.6)
3 i i
1
1X
f2 = log aCa2+ .aSO42− .a2H2 O Ks−1 (3.7)
3 i i
2
avec i=1,2,3. Ces deux relations expriment les indices de saturation moyennés sur les trois
sections de la cuve d’attaque.
69
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
1 % ≤ Sc ≤ 3 % (3.8)
T ≤ 355 K (3.9)
Les contraintes d’égalité sont fournies par les équations du modèle thermodynamique
(chapitre 2). Elles consistent en des bilans de matière (Eq 2.2), de bilan de charge (Eq
2.3) des équations des constantes d’équilibre (Eqs 2.4-2.5) et de l’équation de Pitzer
(Eq. 1.90). En plus, un bilan thermique est utilisé pour prendre en compte l’effet de la
température (Eqs 3.10-3.15) ; il est exprimé par :
où QH , QAG , QR+D , Qout , Qloss et Qcool se réfèrent respectivement à la chaleur des réactifs
d’entrée, à la chaleur d’agitation, à la chaleur de dilution de l’acide sulfurique, à la chaleur
de la suspension quittant le réacteur, aux pertes thermiques de la cuve d’attaque et à la
chaleur à évacuer pour contrôler la température. Les chaleurs des flux d’entrée et de sortie
sont exprimées par :
où msa , Cpsa , Tsa sont respectivement le débit massique, la capacité thermique et la tem-
pérature de l’acide sulfurique. Vps , Cpps , ρps , Tp sont respectivement le débit volumique, la
70
3.3. Problème d’optimisation
Il convient de mentionner que les corrélations développées par Becker (1989) sont
utilisées pour estimer la chaleur d’agitation et la chaleur dégagée par la dilution de l’acide
sulfurique. Il faut noter que Qloss est plus importante pour les petites unités que pour les
grandes en raison de leur grand rapport surface/volume. Sa valeur est souvent trop faible
pour être prise en compte et est donc négligée.
71
72
Tableau 3.1 – Réactions impliquées dans la digestion et leur constante d’équilibre et enthalpie correspondantes à 298 K.
(-) signifie que la réaction est totale et n’implique pas d’équilibre.
1 K 2– +
2 )−
HSO4 – −−*− SO4 + H 0.0103 -1.2.104 Gustafsson (2011) ; Pitzer (2018)
2 K – +
3 −−
H3 PO4 )−*− H2 PO4 + H 0.0071 -4.2.103 Gustafsson (2011)
3 K 2– +
4 )−
H2 PO4 – −−*− HPO4 + H 4.2.10−12 -1.5.104 Gustafsson (2011)
1
4 K –
5 )−
H2 PO4 – + H3 PO4 −−*− H5 P2 O8 0.2550 -9.2.103 Gustafsson (2011)
K5– +
6 −−
HF )−*
−F +H 7.2.10−4 -1.3.104 Ball et Nordstrom (1991)
6 K –
7 )−
HF + F – −−*− HF2 5.5000 -1.7.104 Ball et Nordstrom (1991)
K
7 2– +
8 )−
H2 SiF6 −−*− SiF + 2 H 0.3.102 -6.7.104 Gustafsson (2011)
8 K
9 −−
Ca2+ + SO4 2 – )−*− CaSO4 0.6394 -7.2.103 Gustafsson (2011) ; Pitzer (2018)
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
3.4. Résultats de l’optimisation
gj (x) = 0 j ∈ [1; J]
hk (x) ≥ 0 k ∈ [1; K]
(L) (U )
xi ≤ xi ≤ xi i ∈ [1; N ]
73
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
problèmes d’optimisation. Le temps CPU total est fixé à environ dix-sept heures puisque
le nombre de solutions optimales qui constituent le front de Pareto est égal à 50, ce qui
signifie que le problème d’optimisation multicritère est transformé en 50 problèmes d’op-
timisation mono-critère avec 20 minutes de temps CPU chacun.
Figure 3.2 – Front de Pareto du problème d’optimisation avec quatre points numérotés
de 1 à 4 et détaillés dans le tableau 3.2.
Par ailleurs, quatre solutions sont extraites du front de Pareto et sont détaillées dans le
tableau 3.2. Elles montrent que la distribution optimale d’acide sulfurique dans les trois
sections varie entre 63% et 73% (w1 ) pour la première section, entre 27% et 37% (w2 )
pour la deuxième et 0% (w3 ) pour la dernière section. Ces valeurs sont cohérentes avec les
74
3.4. Résultats de l’optimisation
distributions utilisées à l’échelle industrielle, notamment au sein d’une usine d’acide phos-
phorique située à Jorf Lasfar au Maroc (Bouchkira et Fariq, 2017 ; El Bouzidi, 2016).
Les températures optimales sont également présentées dans le tableau 3.2, leurs valeurs
permettent de déterminer la quantité de chaleur globale à évacuer via l’échangeur de
chaleur pour assurer de bonnes conditions de cristallisation du gypse et aussi pour limiter
la formation de la bassanite qui dégrade les performances de l’unité de filtration située en
aval de la cuve de digestion.
α1
f1max − f1 (x)
u1 (x) = (3.16)
f1max − f1min
α2
f2 (x) − f2min
u2 (x) = (3.17)
f2max − f2min
f1max et f2max sont les valeurs maximales des fonctions objectifs sur le front de Pareto et
f1min et f2min sont leurs valeurs minimales. α1 et α2 sont liés à la tolérance relative sur
les fonctions d’utilité individuelles. La fonction d’utilité multi-attribut est définie comme
suit :
75
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
Les figures 3.3A-B montrent l’effet des tolérances relatives α1 et α2 sur les fonctions d’uti-
lité individuelles. Elles représentent la sensibilité de u1 et u2 par rapport à f1 et f2 . Leurs
valeurs sont fixées dans notre cas à α1 = 0.5 et α2 = 0.5. Les valeurs des paramètres
de pondération y1 et y2 sont prises égales à 0.5 ce qui veut dire que les deux fonctions
objectifs ont la même importance dans la décision.
La fonction d’utilité individuelle U (x) est calculée pour chaque point du front de Pareto et
utilisée pour classer les solutions optimales. La meilleure solution optimale est présentée
dans le tableau 3.3 ainsi que les huit essais effectués à l’échelle industrielle (Bouchkira
et Fariq, 2017 ; El Bouzidi, 2016). Par ailleurs, la figure 3.4 présente toutes les solu-
tions optimales de Pareto avec la meilleure solution (en jaune) et tous les essais industriels.
Il est important de souligner que les pertes totales semblent faibles et que l’on pourrait
être tenté de les ignorer. Cependant, à titre d’exemple, à l’échelle industrielle, les unités de
production d’acide phosphorique située à Jorf Lasfar au Maroc (étudiées dans ce travail)
comprennent 16 unités traitant chacune un débit de 1400 tonnes de minerai de phosphate
par jour. Par conséquent, 1% de perte correspond à 224 tonnes de perte de phosphate par
jour, ce qui est très significatif et ne peut être négligé.
Figure 3.3 – Effet des facteurs de tolérance relative αi . (A) Indice de saturation du gypse,
(B) Indice de saturation de la bassanite.
76
3.4. Résultats de l’optimisation
Pour cette raison, un plan d’expériences (Park, 2007) est utilisé pour développer un
modèle de substitution qui permet de prédire les pertes de phosphate syncristallisés et
inattaquées. Compte tenu de la complexité du procédé et du manque de mesures in-
dustrielles, deux facteurs, c’est-à-dire la température et l’excès d’acide sulfurique sont
supposés connus. Par conséquent, seuls trois facteurs sur cinq, c’est-à-dire w1 , w2 et w3
sont considérés avec deux niveaux.
Les facteurs du plan d’expériences (c’est-à-dire w1 , w2 et w3 ) ainsi que leurs niveaux (c’est-
à-dire L(−1) et L(+1) ) sont indiqués sur le tableau 3.4. Le modèle de substitution adopté
est le suivant :
X XX
Y mod = a0 + bi x i + ci,j xi xj (3.19)
i i j
où Y mod désigne les prédictions du modèle de substitution. Y exp est le vecteur des valeurs
mesurées. X est la matrice du plan d’expérience. xi et xj sont les niveaux des trois
variables considérées. θ est le vecteur des coefficients inconnus du modèle de substitution
donné par :
La figure 3.4 montre que les solutions calculées sont très proches des mesures industrielles
qui sont obtenues en effectuant plusieurs expériences à l’échelle de l’usine, et en chan-
geant les valeurs des variables de décision d’une expérience à l’autre pour améliorer la
productivité et minimiser les pertes de phosphate. De plus, les valeurs des variables de
décision utilisées pour réaliser les expériences industrielles, sont sous-optimales et doivent
être ajustées en fonction de la composition du minerai de phosphate ce qui augmente
encore le nombre des expériences.
Il est important de noter dans ces conditions que l’obtention de ces mesures est très
coûteuse puisqu’elle implique un coût de main d’œuvre important, consomme du temps et
d’énormes quantités de réactifs, diminue la productivité de l’usine et pose des problèmes de
sécurité. La modélisation et l’optimisation multicritère se sont avérées efficaces puisqu’elles
77
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
fournissent des solutions optimales qui minimisent le nombre, la durée et le coût des
expériences et les risques de sécurité.
Figure 3.4 – Comparaison entre les solutions optimales, les pratiques industrielles et la
meilleure solution optimale.
3.5 Conclusion
Les conditions optimales de fonctionnement de la cuve d’attaque d’un procédé industriel
de production d’acide phosphorique sont déterminées à l’aide d’une méthode d’optimisa-
tion multicritère. Les objectifs utilisés (c’est-à-dire les indices de saturation) se sont avérés
efficaces pour minimiser les pertes de phosphate et augmenter la productivité de l’acide.
En outre, les deux contraintes d’inégalité sur la température et sur l’excès d’acide sulfu-
rique au-dessus de la quantité stœchiométrique requise par les réactions sont considérées
et se sont avérées aussi pertinentes pour assurer un fonctionnement optimal de la cuve
d’attaque.
Les résultats de l’optimisation sont cohérents avec les mesures effectuées sur une instal-
lation industrielle, montrant ainsi que l’optimisation est un outil puissant pour réduire le
coût et la durée des expériences industrielles dans des conditions de sécurité maximale.
Cependant, il est à noter que certaines réactions secondaires peuvent se produire en raison
des impuretés présentes dans le minerai de phosphate qui ne sont pas considérées dans
ce travail. Leur prise en compte permettrait sans aucun doute d’améliorer davantage les
performances de la cuve de digestion.
78
3.5. Conclusion
1 0.50 0.20 0.30 0.94 0.06 Essai industriel (El Bouzidi, 2016)
2 0.40 0.20 0.40 0.73 0.12 Essai industriel (El Bouzidi, 2016)
3 0.40 0.30 0.30 0.87 0.12 Essai industriel (El Bouzidi, 2016)
4 0.70 0.30 0 0.62 0.14 Essai industriel (Bouchkira et Fariq, 2017)
5 0.64 0.36 0 0.70 0.15 Essai industriel (Bouchkira et Fariq, 2017)
6 0.75 0.25 0 0.64 0.12 Essai industriel (Bouchkira et Fariq, 2017)
7 0.75 0.17 0.08 0.68 0.14 Essai industriel (Bouchkira et Fariq, 2017)
8 0.60 0.30 0.10 0.63 0.12 Essai industriel (Bouchkira et Fariq, 2017)
9 0.65 0.35 0 0.63 0.14 Meilleure solution optimale
Tableau 3.4 – Niveaux du plan d’expériences : L(+1) and L(−1) . SL and U L se réfèrent aux
pertes syncristallisées et inattaquées.
w1 w2 w3 SL(%) U L(%)
1 +1 +1 -1 0.62 0.14
2 +1 -1 -1 0.64 0.12
3 -1 -1 +1 0.94 0.06
4 -1 +1 +1 0.87 0.12
L(+1) [0.5 ;0.7] [0.2 ;0.3] [0.2 ;0.4] - -
L(−1) [0.4 ;0.5] [0.1 ;0.2] [0 ; 0.2] - -
79
Chapitre 3. Exploitation du modèle thermodynamique pour l’optimisation multicritère du procédé dihydrate
80
Chapitre 4
Sommaire
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Encrassement des unités de production d’acide phosphorique 83
4.2.1 Impuretés liées aux procédés . . . . . . . . . . . . . . . . . . . 84
4.2.2 Impuretés liés au minerai de phosphate . . . . . . . . . . . . . 84
4.3 Modélisation thermodynamique . . . . . . . . . . . . . . . . . . 85
4.3.1 Indice de saturation . . . . . . . . . . . . . . . . . . . . . . . . 86
4.3.2 Phases minérales considérées . . . . . . . . . . . . . . . . . . . 86
4.4 Résultats et discussions . . . . . . . . . . . . . . . . . . . . . . . 87
4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
81
Chapitre 4. Modélisation thermodynamique des problèmes d’encrassement dans la cuve d’attaque
Résumé
Dans le procédé de fabrication d’acide phosphorique par voie humide, plusieurs réactions
secondaires indésirables se produisent en raison des impuretés provenant des réactifs uti-
lisés. Ces réactions entraînent la cristallisation de sels, d’oxydes et d’hydroxydes solides,
ce qui provoque des problèmes d’encrassement au niveau des unités de production. Ces
problèmes nécessitent l’arrêt de la production pour maintenance, ce qui est très coûteux
en termes de rendement des unités, de leur productivité et surtout de leur sécurité. Ce
chapitre présente une exploitation du modèle thermodynamique développé, afin de pré-
dire les indices de saturation de plusieurs produits solides et d’évaluer leur potentiel de
précipitation, pour différentes conditions opératoires. Une approche de sensibilité globale
est par ailleurs utilisée pour définir et identifier les paramètres opératoires qui peuvent
limiter, voire éviter les problèmes d’encrassement au niveau des unités de production. Plus
spécifiquement, trois paramètres sont étudiés : la température de la réaction, le taux de
sulfates libres dans le milieu réactionnel et la concentration en phosphate dans le milieu
réactionnel. Les prédictions du modèle sont cohérentes avec des mesures expérimentales
issues de la littérature.
82
4.1. Introduction
4.1 Introduction
L’encrassement peut être défini comme étant l’accumulation d’éléments solides indési-
rables sur une interface (Guo et al., 2012 ; Hamedi et al., 2019 ; Watkinson et Wilson,
1997). C’est un phénomène indésirable qui affecte un grand nombre d’opérations unitaires
industrielles. Selon la classification établie par Epstein (1997) , il existe cinq grands types
d’encrassement différents, à savoir :
Quel que soit le phénomène d’encrassement considéré, il est plus ou moins critique selon la
nature des fluides employés et/ou la conception des unités industrielles. Dans les situations
réelles d’encrassement, les différents mécanismes ci-dessus interagissent ou se superposent
de manière plus ou moins complexe.
Les impuretés qui sont à l’origine des encrassements et de dépôts peuvent être classées en
deux grandes catégories (Azaroual et al., 2012 ; Khamar, 2012 ; Khamar et al., 2014 ;
Zmemla et al., 2020) :
83
Chapitre 4. Modélisation thermodynamique des problèmes d’encrassement dans la cuve d’attaque
Fluor et silice
Dans les roches phosphatées d’origine sédimentaire, la teneur en fluor varie entre 3% et
4%. Durant le processus de fabrication d’acide phosphorique, environ 20% du fluor sont
évacués avec les gaz, 65% restent sous forme solide et évacués avec le gypse, et 15% sont
solubles dans la solution d’acide phosphorique produite. Par ailleurs, la teneur des roches
phosphatées en silice valorisable varie entre 1% et 6%. Lors du processus d’attaque des
phosphates, la silice se combine avec les autres impuretés contenues dans la phase liquide
pour produire d’autres espèces chimiques aqueuses et solides. En effet, une grande partie
du fluor dissout s’associe à la silice et conduit à la précipitation de composés fluorosilicates
comme les fluorosilicates de sodium N a2 SiF6 , et de potassium K2 SiF6 , et d’autres solides
comme le quartz sous différentes formes cristallines, et la chukhrovite qui est un cristal
complexe contenant principalement des ions Ca2+ , des ions SiF62− , des ions Al3+ et des
ions SO42− .
84
4.3. Modélisation thermodynamique
Les roches phosphatées contiennent aussi des teneurs en fer comprises entre 0.1% et 0.2%.
Ces teneurs restent toujours inférieures à celles de l’aluminium dans la plupart des roches
phosphatées qui peuvent atteindre 1%. La présence de ces deux éléments pendant la
réaction de digestion cause la précipitation de solides comme les (F e ou Al)3 KH14 (P O4 )8
qui provoquent de l’encrassement, mais aussi impactent le rendement chimique en captant
des ions phosphatés.
Potassium et sodium
Les analyses périodiques des phosphates valorisables au niveau des unités de production
montrent que ces derniers contiennent de 0.04 à 0.16% en potassium. Lors de la réaction
de digestion, il se combine avec d’autres impuretés contenues dans la phase aqueuse pour
produire des phases solides notamment à base de fluor et de silice comme décrit précédem-
ment. Par ailleurs, la teneur en sodium varie généralement entre 0.3 et 1.2%. Cet élément
conduit, dans des conditions de température et de concentration dans la phase aqueuse,
à la formation de composés fluorosilicates.
Chlorures
Les chlorures présentent des risques de corrosion lorsque leur teneur dépasse 1000 ppm.
Pour les phosphates valorisés, le suivi des analyses chimiques montre que cette teneur
varie entre 133 et 337 ppm. En général, lors de l’attaque des phosphates, la totalité des
chlorures passe en phase liquide. Vu la difficulté de leur séparation de cette phase et
leur effet fortement corrosif, ils peuvent endommager les différentes unités du procédé de
valorisation.
Calcium
L’élément calcium fait partie des impuretés majeures présentes dans la roche phosphatée,
mais n’est pas critique pour la production. Sa séparation du minerai de phosphate attaqué
est assuré par l’ajout d’acide sulfurique concentré, ce qui produit des sulfates de calcium
sous différentes formes telles que le gypse CaSO4 .2H2 O, la bassanite CaSO4 : 0.5H2 0,
ou l’anhydrite CaSO4 . L’effet de ces produits solide sur l’encrassement des unités dépend
des lieux de leur précipitation.
85
Chapitre 4. Modélisation thermodynamique des problèmes d’encrassement dans la cuve d’attaque
Une revue bibliographique est détaillée dans la partie suivante afin de définir la liste des
phases minérales à prendre en compte dans le cas de la digestion de la roche phosphatée
impure.
Qn
(mi γi )ai
SI = log i=1
(4.1)
Ks
En outre, des analyses quantitatives de diffraction des rayons X concernant le quartz ont
été réalisées et ont prouvé son existence sous plusieurs formes cristallines dont la cristo-
balite et la calcédoine. Elles sont produites principalement en raison de la présence de la
86
4.4. Résultats et discussions
silice, et leur forme cristalline dépend des conditions opératoires, notamment de la tem-
pérature et de la pression (Khamar, 2012 ; Waerstad, 1985 ; Zmemla et al., 2020).
L’ensemble des phases solides qui sont susceptibles de causer des problèmes d’encrasse-
ment sont listés dans le tableau 4.1. Les équilibres sur lesquels le modèle thermodynamique
se base pour calculer les indices de saturation sont listés dans le tableau 4.2.
Les indices de saturation des minéraux à base de calcium (c’est-à-dire le gypse, la bas-
sanite et l’anhydrite) sont positifs, ce qui signifie que ces minéraux sont susceptibles de
cristalliser durant la réaction de digestion du minerai. En effet, au cours du processus de
digestion, les ions calcium libérés par le minerai de phosphate ont tendance à capturer les
ions sulfates libres de même polarité et de charges opposées conduisant à la cristallisation
des sulfates de calcium. Les résultats de l’analyse de sensibilité montrent que les indices
de saturation de ces minéraux sont très sensibles à la température de fonctionnement, ce
qui est très cohérent avec des études précédentes (Peng et Robinson, 1976 ; Shen et al.,
2020 ; Zhu et al., 2016). En général, ces minéraux sont considérés comme des produits
principaux de la réaction de digestion et leur implication dans l’encrassement des unités
de production peut être négligée.
87
Chapitre 4. Modélisation thermodynamique des problèmes d’encrassement dans la cuve d’attaque
(Frazier et al., 1966 ; Zmemla et al., 2020), où les auteurs ont particulièrement étudié
l’existence de ces minéraux dans le gâteau de gypse issu des unités de filtration. Il a été
également démontré dans (Khamar, 2012) que la cristallisation de ces phases minérales
se produit lors du refroidissement de l’acide phosphorique industriel dans les réservoirs
de stockage. L’analyse de sensibilité montre que la production de ces minéraux dépend
principalement de la concentration en soufre, ce qui est cohérent avec les résultats de
Becker (1989) et de Zmemla et al. (2020) .
Les résultats montrent également que le quartz a tendance à précipiter sous différentes
formes cristallines comme la calcédoine et la cristobalite. L’analyse de sensibilité montre
leur grande sensibilité à la température de fonctionnement du procédé. Pour limiter leur
précipitation, les industriels conseillent de doper le mélange réactionnel avec de l’alumi-
nium, ce qui favorise la fixation de la silice et la production d’argile et de la chukhrovite
qui sont favorables à l’amélioration de la filtrabilité du gâteau de gypse compte tenu de
leur forme cristalline (Khamar, 2012 ; Khamar et al., 2014 ; Zmemla et al., 2020).
Enfin, les indices de saturation des minéraux à base de magnésium sont positifs et dé-
pendent principalement de la concentration en soufre. En effet, le magnésium est une des
impuretés critiques impliquées dans plusieurs réactions secondaires qui provoquent des
phénomènes d’encrassement et diminuent la filtrabilité du gâteau de gypse. Il est à noter
que ces minéraux contenant du magnésium sont généralement caractérisés par de très
petites tailles, ce qui favorise les phénomènes d’encrassement (Bouchkira et al., 2021d ;
Khamar, 2012) et qui minimisent la porosité du gâteau de gypse et ainsi dégradent sa
filtrabilité (Becker, 1983).
4.5 Conclusion
Dans ce chapitre, le modèle thermodynamique développé est utilisé pour prédire la suscep-
tibilité des minéraux à précipiter pendant la digestion du minerai de phosphate en raison
de réactions secondaires et à conduire à des phénomènes d’encrassement. Il a été montré
que dans les conditions actuelles de production, la malladrite, la hiératite, la magadiite,
le quartz et l’hexafluorosilicate de magnésium sous leurs différentes formes cristallines
sont les minéraux les plus critiques susceptibles de précipiter et de provoquer des phé-
88
4.5. Conclusion
Figure 4.1 – Prédiction des indices de saturation des différentes phases minérales consi-
dérées.
89
Chapitre 4. Modélisation thermodynamique des problèmes d’encrassement dans la cuve d’attaque
Figure 4.2 – Indice de sensibilité globale de Sobol. (A) : Sensibilité par rapport à la
concentration d’acide sulfurique. (B) : Sensibilité par rapport à la concentration d’acide
phosphorique. (C) : Sensibilité par rapport à la température.
90
Tableau 4.1 – Phases minérales et leur réaction de formation-dissolution.
2+
Fluorapatite ) HF + 3PO3−
Ca5 (PO4 )3 F + H+ * 4 + 5Ca -0.910 (Ball et Nordstrom, 1991 ; Gustafsson, 2011)
3+
Chukhrovite ) 4Ca2+ + SiF2−
Ca4 SO4 AlSiF6 * 6 + Al + SO4 2− -19.55 (Zmemla et al., 2020)
+
Malladrite ) SiF2−
Na2 SiF6 * 6 + 2Na -6.475 (Skafi et al., 2019)
+
Hiératite ) SiF2−
K2 SiF6 * 6 +2 K -7.039 (Skafi et al., 2019)
Fluorite ) Ca+2 + 2 F−
CaF2 * -10.96 (Ball et Nordstrom, 1991 ; Gustafsson, 2011)
Quartz ) H4 SiO4
SiO2 + 2H2 O * -3.018 (Ball et Nordstrom, 1991 ; Gustafsson, 2011)
2+
MgSiF6 ) SiF−2
MgSiF6 * 6 + Mg -9.960 (Gustafsson, 2011 ; Khamar, 2012)
2+
MgSiF6 · 6H2 O ) SiF−2
MgSiF6 : 6H2 O * 6 + Mg + 6H2 0 -12.96 (Gustafsson, 2011 ; Khamar, 2012)
2+
MgSO4 .7H2 O ) SiF−2
MgSiF6 : 7H2 O * 6 + Mg + 7H2 0 -14.96 (Gustafsson, 2011 ; Khamar, 2012)
91
4.5. Conclusion
92
Tableau 4.2 – Phases aqueuses utilisées par le modèle développé
* 2− +
1 HSO−
4 ) SO4 + H 0.0103 -1.274 .104 (Gustafsson, 2011 ; Pitzer, 2018a)
+
2 H3 PO4 *
) H2 PO−
4 +H 0.0071 -4.271 .103 (Gustafsson, 2011)
* 2− +
3 H2 PO−
4 ) HPO4 + H 6.3 .10−8 -1.800 .104 (Gustafsson, 2011)
* 3− +
4 HPO2−
4 ) PO4 + H 4.2 .10−12 -1.500 .104 (Gustafsson, 2011)
* −
5 H3 PO4 + H2 PO−
4 ) H 5 P2 O 8 0.2550 -9.266 .103 (Khamar, 2012)
6 HF *
) F− + H + 7.2 .10−4 -1.329 .104 (Gustafsson, 2011)
+
7 ) SiF2−
H2 SiF6 * 6 + 2H 0.3 .102 -4.799 .104 (Gustafsson, 2011)
8 H4 SiO4 *
) SiO2 + 2H2 O 3.52 .10−7 -1.121 .105 (Gustafsson, 2011)
9 CaF+ *
) Ca2+ + F− 1.148 .10−7 -1.255 .104 (Ball et Nordstrom, 1991)
10 ) Ca2+ + SO2−
CaSO4 * 4 0.6394 -7.200 .103 (Pitzer, 2018a ; Shen et al., 2020)
11 AISO+ * 3+ + SO2−
4 ) Al 4 9.549 .10−4 -8.368 .103 (Gustafsson, 2011)
12 AlF− * 3+ + 4 F−
4 ) Al 9.772 .10−8 -8.368 .103 (Gustafsson, 2011)
13 ) Mg2+ + SO2−
MgSO4 * 4 5.623 .10−3 -4.184 .103 (Ball et Nordstrom, 1991)
* + 2−
14 NaSO−
4 ) Na + SO4 0.1995 -4.170 .103 (Ball et Nordstrom, 1991)
* + 2−
15 KSO−
4 ) K + SO4 1.1412 -8.368 .103 (Ball et Nordstrom, 1991)
Chapitre 4. Modélisation thermodynamique des problèmes d’encrassement dans la cuve d’attaque
Chapitre 5
Modélisation de la cristallisation du
gypse dans un réacteur batch
Sommaire
5.1 Contexte de l’étude . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.2 Travaux de la littérature . . . . . . . . . . . . . . . . . . . . . . 95
5.3 Modèle de bilan de population . . . . . . . . . . . . . . . . . . 97
5.3.1 Nucléation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5.3.2 Croissance cristalline . . . . . . . . . . . . . . . . . . . . . . . . 99
5.3.3 Bilan de matière . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5.4 Méthodes de résolution . . . . . . . . . . . . . . . . . . . . . . . 100
5.4.1 Méthode des moments (1D) . . . . . . . . . . . . . . . . . . . . 100
5.4.2 Méthode des différence finies (1D) . . . . . . . . . . . . . . . . 101
5.4.3 Méthode des différences finies (2D) . . . . . . . . . . . . . . . . 102
5.5 Résultats et discussions . . . . . . . . . . . . . . . . . . . . . . . 104
5.5.1 Estimbilité et identification . . . . . . . . . . . . . . . . . . . . 104
5.5.2 Méthode des moments (1D) . . . . . . . . . . . . . . . . . . . . 107
5.5.3 Méthode des différences finies (1D) . . . . . . . . . . . . . . . . 109
5.5.4 Méthode des différences finies (2D) . . . . . . . . . . . . . . . . 110
5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
93
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
Résumé
Le présent chapitre traite du développement d’un modèle de cristallisation du gypse lors
de la production d’acide phosphorique. Le modèle de bilan de population impliquant les
taux de nucléation et de croissance des cristaux est utilisé. Le modèle thermodynamique
développé est exploité pour prédire la sursaturation du milieu réactionnel, ce qui permet
le calcul des termes sources du modèle de bilan de population. Une base de données qui
contient des valeurs du rapport de sursaturation du gypse est utilisée pour identifier les
paramètres inconnus du modèle. Ces derniers sont utilisés pour comparer les prédictions
du modèle avec les mesures expérimentales. Les résultats sont prometteurs, et le modèle
identifié peut être utilisé dans la conception et l’optimisation de la cristallisation du gypse.
94
5.1. Contexte de l’étude
En lien avec la cristallisation des sulfates de calcium, El Moussaouiti et al. (1997) ont
couplé la PBE monodimensionnelle avec une équation empirique qui permet de calculer
la sursaturation. Leur modèle permet de prédire l’évolution de la distribution des tailles
95
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
Il est à noter que plusieurs travaux de la littérature se focalisent plutôt sur les approches
de résolution des systèmes d’équations résultants du couplage de la PBE avec les modèles
thermodynamiques. Parmi ces approches l’on peut citer, la méthode des moments, la mé-
thode des volumes finis, la méthode des éléments finis, et la méthode des caractéristiques.
En effet, (Andalibi et al., 2018 ; Andalibii et al., 2020 ; Li et al., 2019) ont proposé des
approches basées sur la méthode des moments comme une classe de méthodes qui per-
mettent de résoudre la PBE en la transformant en un système d’équations différentielles
ordinaires. Liu et Rigopoulos (2019) ont utilisé une méthode des volumes finis pour la
discrétisation des équations du modèle de population. La comparaison de leurs résultats
avec des solutions analytiques a montré la qualité de leur approche. Dans (Mantzaris
et al., 2001 ; Tsang et Rao, 1990) la méthode des éléments finis a été utilisée pour la
résolution de la PBE. Finalement, plusieurs travaux ont montré l’efficacité de la méthode
des caractéristiques pour la résolutions de plusieurs types de PBE (Peng et al., 2015 ;
Ramkrishna, 2000).
96
5.3. Modèle de bilan de population
Pour déterminer l’expression de ζ, il faut faire des hypothèses sur le processus de cristal-
lisation. Dans ce travail, il est supposé que la cristallisation est mise en œuvre dans un
réacteur batch isotherme et parfaitement mélangé, et que la vitesse de croissance est indé-
pendante de la taille des cristaux. Dans ces conditions, l’équation du bilan de population
s’écrit (Hanhoun et al., 2013) :
— Cas monodimensionnel (1D) :
∂n(t, x) ∂n(t, x)
= −G(t) (5.2)
∂t ∂x
RN (t)
avec : n(t, 0) = (5.3)
G(t)
RN (t)
avec : n(t, 0, y) = (5.5)
Gx (t)
RN (t)
et : n(t, x, 0) = (5.6)
Gy (t)
97
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
5.3.1 Nucléation
La nucléation est le phénomène de naissance des cristaux à partir d’une solution sursa-
turée. Ce phénomène impose le nombre final de cristaux, il a donc un effet direct sur la
distribution de leur taille. Deux types de nucléation peuvent être distingués : la nucléa-
tion primaire et la nucléation secondaire. En effet, dans une solution exempte de cristaux,
les molécules de soluté diffusent, et certaines finissent par se rencontrer et forment des
multimères (complexes macromoléculaires). À partir d’une certaine taille, ces agrégats
sont considérés comme des germes cristallins. On peut alors définir la nucléation primaire
comme étant l’apparition des germes dans le milieu où il n’existe aucun cristal (Chow
et al., 2005 ; El Bazi, 2011).
RN = RN 1 + RN 2 (5.7)
Ne tenant pas en compte l’effet de la température (réacteur isotherme), les termes sources
de nucléation peuvent s’exprimer par les modèles suivants (Barbier et al., 2009 ; Hanhoun
et al., 2013) :
−B
RN1 (t) = A × exp , nucléation primaire (5.8)
ln S(t)2
98
5.3. Modèle de bilan de population
!
aCa2+ aSO42− a2H2 O
où : S= ; avec : ai = mi .γi , pour la dissociation du gypse (Eq.3.5)
Ksp
(5.10)
Z ∞
et : µ3 (t) = x3 n(t, x)dx, pour le cas (1D) (5.11)
0
Z ∞ Z ∞
et : µ3 (t) = xy 2 n(t, x, y)dxdy, pour le cas (2D) (5.12)
0 0
Dans les équations ci-dessus, A, B, ks et s sont des paramètres inconnus qui sont géné-
ralement identifiés à partir de mesures expérimentales, S est le rapport de sursaturation
qui dépend notamment de la spéciation du milieu réactionnel, c’est-à-dire des mi et des
coefficients d’activité γi et de Ksp qui désigne le produit de solubilité du gypse. µ3 est le
moment d’ordre trois de la fonction de densité (n).
Ces termes source de nucléation et de croissance cristallines sont pris en compte dans
l’étude de la cristallisation du gypse lors de la production d’acide phosphorique. Leur
prédiction s’appuie sur les équations de bilans de matière et le modèle thermodynamique
précédemment développé. Ces équations de bilan sont établies ci-dessous.
99
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
1
0
CCa3 (PO4 )2 (t) = CCa3 (PO4 )2
− C∆gyp (t) (5.15)
3
2
CH3 PO4 (t) = CH0 3 PO4 + C∆gyp (t) (5.18)
3
où Ci0 sont les concentrations initiales des réactifs. Les équations (5.15-5.18) doivent être
résolues à chaque pas de temps afin de prédire le rapport de sursaturation à l’aide du
modèle thermodynamique développé. Le moment d’ordre trois de la distribution de la taille
des cristaux permet de calculer la quantité de gypse produite dans le temps (Hanhoun
et al., 2013 ; Peng et al., 2015 ; Zhou et al., 2013). Elle est donnée par :
kv ρs
C∆gyp (t) = (µ3 (t) − µ3 (0)) (5.19)
Mc
où kv est un facteur de forme volumique, ρs est la densité des cristaux et Mc est la masse
molaire du gypse.
Z ∞
µk (t) = xk n(t, x)dx (5.20)
0
100
5.4. Méthodes de résolution
dµ0 (t)
= RN (t) (5.21)
dt
dµk (t)
= kG(t)µk−1 (t), k ≥ 1 (5.22)
dt
avec les conditions initiales suivantes :
Z ∞
µk (0) = xk n0 (x)dx (5.23)
0
X = {x0 , x1 , . . . , xN } (5.24)
avec :
xi ∈ X, xi − xi−1 = ∆x (5.25)
101
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
RN (t)
n(t, 0) = (5.28)
G(t)
RR
— Les conditions initiales sont supposées suivre une distribution normale et sont ex-
primées par :
" 2 #
1 1 xi − µ x
n(0, xi ) = √ exp − (5.29)
2π 2 σx
Enfin, pour le calcul du moment d’ordre 3, la méthode des trapèzes est utilisée. Ainsi, le
moment s’exprime comme suit :
N −1
∆x X 3
xi n (xi , t) + x3i+1 n (xi+1 , t) (5.30)
µ3 (t) =
2 i=0
X = {x0 , x1 , . . . , xN } (5.31)
Y = {y0 , y1 , . . . , yN } (5.32)
102
5.4. Méthodes de résolution
Avec :
xi ∈ X, xi − xi−1 = ∆x (5.33)
yi ∈ Y, yj − yj−1 = ∆y (5.34)
Comme pour le cas (1D), l’équation du bilan de population discrétisée aux points (xi , yj )
s’écrit comme suit :
— Pour i et j = 1 à N-1 :
RR
— Pour i=N et j=N :
RN (t)
n(t, 0, yj ) = , ∀j (5.37)
Gx (t)
— Pour j=0 :
RN (t)
n(t, xi , 0) = , ∀i (5.38)
Gy (t)
RR
— Les conditions initiales sont supposées suivre la distribution normale suivante :
" 2 2 #
1 1 xi − µ x 1 y j − µy
n(0, xi , yj ) = √ exp − − (5.39)
2π 2 σx 2 σy
Enfin, comme pour le cas (1D), le moment d’ordre trois est approximé à l’aide de la
méthode des trapèzes détaillée dans (Peng et al., 2015).
103
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
Afin d’éviter le choix d’un seuil d’estimabilité arbitraire, l’approche d’estimabilité globale
développée en chapitre 2 est utilisée. Un seuil optimal d’estimabilité est déterminé pour
chaque cas étudié. Pour le premier cas où la PBE monodimensionnelle est résolue à l’aide
de la méthode des moments, l’approche d’estimabilité globale a permis de définir un seuil
d’estimabilité optimal d’environ 0.013. Cette valeur signifie que pour qu’un paramètre
soit estimable, une variation de 10% de sa valeur doit provoquer au moins 1.14% de va-
riation au niveau de la sortie du modèle. La comparaison des seuils d’estimabilité des
paramètres inconnus (tableau 5.1) montre que tous les paramètres sont estimables. Pour
le deuxième cas où la PBE monodimensionnelle est résolue par la méthode des différences
finies, le seuil optimal d’estimabilité est estimé à environ 0.015. Cette valeur signifie que
pour qu’un paramètre soit estimable, une variation de sa valeur de 10% doit provoquer au
moins 1.34% de variation de la sortie du modèle. La comparaison des seuils d’estimabilité
des paramètres avec cette valeur optimale (Tableau 5.2) montre aussi que l’ensemble des
paramètres sont estimables. Concernant le cas de la PBE bidimensionnelle, l’approche
d’estimabilité globale a permis de définir un seuil d’estimabilité globale d’environ 0.011.
Cette valeur signifie que pour qu’un paramètre soit estimable, une variation de sa valeur
de 10% doit provoquer au moins 1.04% de variation de la sortie du modèle. La compa-
raison des seuils d’estimabilité des paramètres avec cette valeur optimale (tableau 5.3),
montre également que la majorité des paramètres sont estimables sauf gy qui a un seuil
d’estimabilité inférieur au seuil optimal, ce qui signifie qu’il n’est pas estimable à partir
des mesures disponibles. Sa valeur est par conséquent fixée à 2 à partir des travaux de la
littérature (Barbier et al., 2009 ; Peng et al., 2015 ; Zhu et al., 2016).
Les paramètres estimables sont ensuite identifiés et leur précision est caractérisée à l’aide
des intervalles de confiance. Les valeurs des paramètres A et B montrent que la nucléation
primaire est très faible, et peut même être négligée, ce qui est en accord avec les études
de (Becker, 1983 ; Becker, 1989 ; Peng et al., 2015 ; Zhu et al., 2016) qui ont montré
que le phénomène de nucléation primaire est négligeable dans le cas de la cristallisation
des sulfates de calcium en réacteur batch.
104
Tableau 5.1 – Estimabilité et identification. : Cas de PBE monodimensionnelle (méthode des moments).
Paramètres A B ks kg g s
Rang d’estimabilité 6 4 2 1 3 5
Valeur identifiée 5.32.10−6 5.001 1.17 10−4 0.746 10−4 2.01 2.04
Intervalle de confiance [0.83 ;9.81].10−6 [4.35 ;5.65] [0.13 ;2.21].10−4 [0.55 ;0.94].10−4 [1.87 ;2.15] [1.66 ;2.36]
Tableau 5.2 – Estimabilité et identification. Cas de PBE monodimensionnelle (méthode des différences finies).
Paramètres A B ks kg g s
Rang d’estimabilité 6 4 2 1 3 5
Valeur identifiée 4.93.10−6 4.05 3.76 10−4 1.76 10−6 1.23 1.52
Intervalle de confiance [4.53 ;5.33].10−6 [1.36 ;6.74] [2.28 ;5.24].10−4 [1.03 ;2.49].10−6 [0.44 ;2.02] [0.94 ;2.10]
105
5.5. Résultats et discussions
106
Tableau 5.3 – Résultats d’estimabilité et d’identification. Cas de PBE bidimensionnelle (méthode des différences finies).
Paramètres A B ks kg x gx s
Rang d’estimabilité 7 5 2 1 4 6
Valeur identifiée 3.84.10−8 4.36 7.65 10−4 8.48 10−5 2.13 2.01
Intervalle de confiance [3.36 ;4.32].10−8 [1.98 ;6.74] [7.06 ;8.24].10−4 [7.47 ;9.49].10−6 [2.07 ;2.19] [0.92 ;3.10]
Paramètres kg y gy
Rang d’estimabilité 3 8
Unité m3 s−1 -
−6 −4
Intervalles de confiance AA [0.16 ;9.81].10 A [4.34 ;5.65] A [0.30 ;6.21].10 A [0.27 ;3.34].10−4 [1.32 ;4.15] A [1.92 ;5.36]
5.5. Résultats et discussions
Figure 5.1 – Résultats de modélisation. (A) : Diagramme de parité pour les mesures
d’équilibre. (B) : Diagramme de parité pour les données de sursaturation.
La figure 5.2 montre la concentration des réactifs (acide sulfurique et phosphate trical-
cique), et des produits (acide phosphorique et gypse). Les résultats sont cohérents puisque
les quantités d’acide phosphorique et de gypse augmentent avec le temps, tandis que celles
d’acide sulfurique et de phosphate tricalcique diminuent. De plus, comme attendu, le ta-
bleau 5.4 montre que la taille moyenne des cristaux de gypse augmente avec le temps, ce
qui signifie que la croissance des particules a bien lieu dans le réacteur. D’autre part, les
paramètres des distributions des tailles des particules sont calculés en utilisant les valeurs
des moments centrés d’ordre k (µ∗k ) à différents
√ ∗ temps de cristallisation. Ainsi, l’écart-
type est déduit du moment d’ordre 2, σ = µ2 , le coefficient de dissymétrie est calculé en
∗3/2
utilisant le moment d’ordre 3 sous la forme α = µ∗3 /µ2 , et le coefficient d’aplatissement
est calculé à partir du moment d’ordre 4 comme suit : β = µ∗4 /µ∗2 2 .
107
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
Figure 5.3 – Distribution des tailles des particules normées à différents instants.
Il est intéressant de noter que les valeurs de ces paramètres ne changent pas au cours du
temps et que les deux dernières valeurs sont typiques d’une distribution normale, c’est-
à-dire α = 0 et β = 3. De plus, étant donné que la distribution initiale des tailles des
particules est supposée normale et que seules les cinétiques de nucléation et de croissance
sont considérées dans le réacteur, les paramètres σ, α et β restent inchangés comme illus-
tré dans la figure 5.3.
Le modèle développé montre qu’il est capable de prédire les variables pertinentes pour la
conception, l’optimisation et le contrôle de la cristallisation du gypse. De plus, il peut être
utilisé pour déterminer les conditions opératoires optimales qui permettent d’extraire la
quantité maximale de phosphate du minerai, et de mieux valoriser le gypse. Bien que les
résultats obtenus soient intéressants, le modèle de cristallisation développé à ce stade est
un modèle monodimensionnel et ne permet pas de prédire les distributions des tailles des
particules. Ces dernières sont déterminées par des approches de reconstruction, qui sont
sujettes à des incertitudes et par conséquent influencent la fiabilité du modèle développé
(John et al., 2007).
108
5.5. Résultats et discussions
Tableau 5.4 – Variations de la taille moyenne des cristaux de gypse avec le temps - Méthode
des moments (1D).
Temps (min) AAAAA 0 AAA 20AAA 25AAA 50AAA 100 AAA
x(µm)AAAAA 10AAA 34.46AAA 39.87AAA 59.31AAA 72.84AAA
Les valeurs des paramètres inconnus du modèle sont ainsi identifiées et présentées dans
le tableau 5.2. Elles sont ensuite utilisées dans la modélisation de la cristallisation du
gypse dans les mêmes conditions opératoires que précédemment. Il est également supposé
que le réacteur contient la même quantité de gypse, dont la distribution de taille des
particules est normale avec x = 30 µm et σ = 7 µm. La figure 5.4 présente les prédictions
de consommation d’acide sulfurique et de phosphate tricalcique, et la production d’acide
phosphorique et de gypse. Les résultats sont également cohérents et en bon accord avec
les résultats obtenus à l’aide de la méthode des moments.
Figure 5.4 – Résultats des différences finies centrées. (A) : Concentrations du phosphate
tricalcique et d’acide sulfurique. (B) : Concentrations du gypse et d’acide phosphorique
Tableau 5.5 – Variations de la taille moyenne des cristaux de gypse avec le temps - méthode
des différences finies (1D).
Time (min) AAAAA 0 AAA 16AAA 33AAA 83AAA 100 AAA
x(µm)AAAAA 30.00 AAA 42.83AAA 54.90AAA 73.97AAA 76.29 AAA
Par ailleurs, le tableau 5.5 montre que la taille moyenne des cristaux de gypse augmente
avec le temps, ce qui signifie que la croissance des particules a bien lieu dans le réacteur.
109
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
D’autre part, les paramètres des distributions des tailles des particules sont calculés à
différents temps de cristallisation. En particulier, l’écart-type des distributions des tailles
des particules à différents instants est calculé et vaut σ = 7 µm, les facteurs de dissymétrie
et d’aplatissement sont également calculés et ont respectivement comme valeurs 1.10−5 et
3. Il est intéressant de noter que les valeurs de ces paramètres ne changent pas au cours du
temps et que les deux dernières valeurs sont typiques d’une distribution normale. Ce qui
est davantage illustré sur la figure 5.5, qui présente la distribution des tailles des particules
dans le temps.
Figure 5.5 – Distribution des tailles des particules (normée) dans le temps.
Le modèle développé montre qu’il est capable de prédire les variables pertinentes pour
la conception, l’optimisation et le contrôle des processus de cristallisation du gypse. Ce-
pendant, il ne permet pas de prédire la non-homogénéité de la cristallisation, et plus
globalement la forme des cristaux et l’influence des impuretés sur leur croissance. Un
modèle bidimensionnel permettrait de prédire la forme plane des cristaux. Ce modèle est
développé ci-dessous.
110
5.5. Résultats et discussions
Les valeurs des paramètres inconnus du modèle pour ce cas sont identifiées et présentées
dans le tableau 5.3. Elles sont ensuite utilisées dans la modélisation de la cristallisation
du gypse dans les mêmes conditions opératoires des cas précédents. Il est supposé égale-
ment que le réacteur contient la même quantité de gypse, dont la distribution de taille
des particules est normale avec x = 30 µm et σ = 7 µm.
La figure 5.6 montre les profils des concentrations des réactifs et des produits. Ces profils
sont cohérents avec les résultats obtenus dans les deux cas précédents. Cependant, les
réactifs se consomment d’une manière plus rapide en les comparant aux cas précédents,
ce qui fait aussi que l’acide phosphorique et les sulfates de calcium se produisent d’une
manière aussi rapide. En effet, pour les deux cas précédents monodimensionnels, les ré-
actifs se consomment pour produire du gypse qui croît dans une seule direction, alors
que dans le cas bidimensionnel, les réactifs se consomment pour produire du gypse qui
croît dans les deux directions. Ceci explique la différence de consommation des réactifs
dans ce cas, et montre que les deux modèles développés en (1D) sous-estiment les taux
de consommation et de production.
Figure 5.6 – Résultats des différences finies centrées (2D). (A) : Concentrations du phos-
phate tricalcique et d’acide sulfurique. (B) : Concentrations du gypse et d’acide phospho-
rique
111
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
Figure 5.7 – Distributions des tailles des particules normées en fonction du temps. (A1-
B1) : à t=0 min. (A2-B2) : à t=10 min. (A3-B3) : à t=30 min. (A4-B4) : à t=100 min.
112
5.6. Conclusion
Tableau 5.6 – Variations de la taille moyenne des cristaux de gypse avec le temps - méthode
des différences finis (2D).
Time (min) AAAAA 0 AAA 10 AAA 30AAA 83AAA 100 AAA
x(µm)AAAAA 30.00 AAA 44.98AAA 51.22AAA 51.81AAA 53.01 AAA
y(µm)AAAAA 30.00 AAA 32.43AAA 33.52AAA 40.18AAA 43.64 AAA
Les figures 5.7(A-B, 1-4) montrent l’évolution de la distribution des tailles des particules
dans le temps selon deux dimensions. En effet, les résultats montrent la croissance des cris-
taux dans les deux directions, avec un taux de croissance de la longueur plus important que
celui de la largeur. Ceci est très cohérent avec les résultats expérimentaux de (Becker,
1983 ; Khamar et al., 2014 ; Peng et Robinson, 1976 ; Zhang et Li, 2015) qui ont
prouvé que la forme finale des cristaux de gypse produits lors de la fabrication d’acide
phosphorique est parallélépipédique. En plus, il est intéressant de caractériser ces distribu-
tions dans le temps à l’aide des paramètres (x, y, σx , σy , αx , αy , βx , βy ) qui correspondent
respectivement aux moyennes de la longueur et de la largeur, aux écarts-types, et aux fac-
teurs de dissymétrie et d’aplatissement (tableau 5.3). L’analyse des valeurs moyennes de
x et y montrent effectivement la croissance des cristaux dans les deux directions (tableau
5.6). Les écarts-types gardent des valeurs constantes (σx = σy = 7 µm) ce qui signifie
que l’ensemble des cristaux croissent dans les deux directions, les facteurs de dissymétrie
sont très faibles (de l’ordre de 10−5 ) et les facteurs d’aplatissements, βx = βy = 3, sont
typiques pour une distribution normale.
5.6 Conclusion
Le modèle thermodynamique développé dans le chapitre 2 est couplé à un modèle de
bilan de population afin d’étudier la cristallisation du gypse pendant la production de
l’acide phosphorique. Dans une première étape, l’étude de la cristallisation en une seule
dimension a été considérée, et le modèle résultant a été résolu à l’aide de la méthode des
moments. Les résultats obtenus sont prometteurs, mais la méthode de résolution présente
quelques limites. En effet, elle ne permet pas de prédire directement la distribution des
tailles des particules, et nécessite une approche de reconstruction basée sur les moments
calculés, ce qui introduit davantage d’incertitudes, et rend la méthode moins fiable. Afin
de surmonter ce problème, la méthode des différences finies monodimensionnelle a été
utilisée comme alternative. Les distributions des tailles des particules sont directement
déterminées à partir de la résolution du modèle de cristallisation. Finalement, afin d’amé-
liorer les résultats et les prédictions du modèle de cristallisation développé, une dernière
partie est considérée où on prend en compte deux dimensions des cristaux de gypse, c’est-
à-dire la longueur et la largeur. Les résultats sont cohérents avec la littérature concernant
la forme finale des cristaux, et des différences ont pu être remarquées, en particulier par
rapport au taux de consommation des réactifs et de synthèse des produits.
113
Chapitre 5. Modélisation de la cristallisation du gypse dans un réacteur batch
Il convient cependant de noter que bien que le présent modèle permette de prédire des
informations importantes et fiables pour la cristallisation du gypse, il ne prend pas en
compte l’influence des impuretés du minerai de phosphate et de l’eau de procédé (non
traitée). Ces impuretés ont beaucoup d’effets sur la croissance des cristaux, et leur prise
en considération permettra sans doute d’améliorer davantage le modèle de cristallisation
développé.
114
Chapitre 6
Sommaire
6.1 Protocole expérimental . . . . . . . . . . . . . . . . . . . . . . . 117
6.2 Plans d’expériences . . . . . . . . . . . . . . . . . . . . . . . . . 120
6.3 Procédure expérimentale et analyses . . . . . . . . . . . . . . . 123
6.4 Analyse de la composition chimique . . . . . . . . . . . . . . . 124
6.4.1 Composition de la phase liquide . . . . . . . . . . . . . . . . . . 124
6.4.2 Composition de la phase solide . . . . . . . . . . . . . . . . . . 124
6.5 Validation du modèle développé . . . . . . . . . . . . . . . . . . 125
6.5.1 Comparaison entre les prédictions et les expériences . . . . . . 125
6.5.2 Test de Kolmogorov-Smirnov . . . . . . . . . . . . . . . . . . . 130
6.5.3 Test F de Fisher-Snedecor . . . . . . . . . . . . . . . . . . . . . 132
6.5.4 Analyses DRX et indices de saturation . . . . . . . . . . . . . . 134
6.6 Développement d’une interface graphique . . . . . . . . . . . . 138
6.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
115
Chapitre 6. Validation expérimentale du modèle thermodynamique
Résumé
Ce chapitre présente la partie expérimentale mise en œuvre pour valider le modèle ther-
modynamique développé dans le chapitre 2. Le montage expérimental et les méthodes
d’analyse sont décrits avant de comparer les mesures aux prédictions du modèle. Le test
statistique F de Fisher-Snedecor est ensuite utilisé pour valider les prédictions du modèle,
après avoir vérifié la normalité des distributions à l’aide du test de Kolmogorof-Smirnov.
116
6.1. Protocole expérimental
Ce chapitre présente les mesures expérimentales réalisées pour notamment valider le mo-
dèle thermodynamique développé dans le chapitre 2. Ces mesures auraient dû être effec-
tuées plutôt dans le travail, mais malheureusement la crise sanitaire nous a contraints de
modifier le planning initialement établi. De plus, seule une partie des mesures initialement
prévues a pu être mise en œuvre dès que nous avions accès aux locaux du laboratoire, et
surtout aux moyens d’analyse.
Le montage expérimental est constitué d’un réacteur batch immergé dans un bain thermo-
staté à circulation (Figure 6.1). La température du bain est réglée par un thermostat, qui
est également muni d’un souffleur d’air pour assurer l’homogénéisation du liquide chauf-
fant. Le réacteur, dont le volume est de 2 litres, est en polyéthylène rigide pour résister
aux acides utilisés. L’agitation est assurée par un agitateur à hélice digitale. À la fin de la
réaction, un dispositif de filtration est mis à disposition ; il est constitué d’une pompe à
vide régulée, un erlenmeyer et un entonnoir filtrant. La filtration permet la séparation des
produits de la réaction, c’est-à-dire, l’acide phosphorique et les sulfates de calcium. Ces
produits subissent plusieurs analyses physiques et chimiques par la suite. Il est finalement
à noter que la réaction de synthèse de l’acide phosphorique est effectuée sous une hotte
qui permet d’aspirer les gaz de la réaction, notamment ceux produits par l’utilisation de
l’acide sulfurique concentré. Le tableau 6.1 contient l’ensemble des réactifs utilisés pour
la fabrication de l’acide phosphorique à l’échelle du laboratoire.
117
Chapitre 6. Validation expérimentale du modèle thermodynamique
Figure 6.1 – Montage expérimental. (1) : Bain thermostaté. (2) : Agitateur. (3) : Réacteur
batch. (4) : Thermostat. (5) : Ampoule d’acide sulfurique. (6) : Pompe à vide. (7) :
Dispositif de filtration. (8) : Hotte de laboratoire.
118
Tableau 6.1 – Description des réactifs utilisés pour la réalisation des expériences
Phosphate tricalcique pur Solide 99% massique Loba Chemie 100 g par expérience
Acide sulfurique Liquide 98% massique Loba Chemie 1%, 2% et 3% massique (Excès)
Acide phosphorique Liquide 85% massique Loba Chemie 16%, 18% et 20% massique
Oxide de magnesium Solide 99% massique Sigma Aldrich 0%, 0.5% et 1% massique
Oxide d’aluminium Solide 99% massique Sigma Aldrich 0%, 0.5% et 1% massique
119
6.1. Protocole expérimental
Chapitre 6. Validation expérimentale du modèle thermodynamique
Un deuxième plan d’expériences est mis au point afin d’évaluer les prédictions du modèle
développé en présence de trois impuretés, à savoir l’oxyde de magnésium, l’oxyde d’alu-
minium et la silice. Ces impuretés sont majoritaires dans le minerai de phosphate. De
plus, trois niveaux pour chaque facteur sont considérés. Le tableau 6.3 décrit les facteurs
et les niveaux du deuxième plan. Les conditions opératoires des expériences à effectuer
sont détaillées dans le tableau 6.5.
120
Tableau 6.4 – Conditions opératoires pour le premier plan d’expériences
Expérience Niveau H2 SO4 Niveau H3 P O4 Niveau T (°C) M gO (%) Al2 O3 (%) SiO2 (%) Temps (min) Agitation (rpm)
1 -1 -1 -1 - - - 120 120
2 0 -1 -1 - - - 120 120
3 +1 -1 -1 - - - 120 120
4 -1 0 -1 - - - 120 120
5 0 0 -1 - - - 120 120
6 +1 0 -1 - - - 120 120
7 -1 +1 -1 - - - 120 120
8 0 +1 -1 - - - 120 120
9 +1 +1 -1 - - - 120 120
10 -1 -1 0 - - - 120 120
11 0 -1 0 - - - 120 120
12 +1 -1 0 - - - 120 120
13 -1 0 0 - - - 120 120
14 0 0 0 - - - 120 120
15 +1 0 0 - - - 120 120
16 -1 +1 0 - - - 120 120
17 0 +1 0 - - - 120 120
18 +1 +1 0 - - - 120 120
19 -1 -1 +1 - - - 120 120
20 0 -1 +1 - - - 120 120
21 +1 -1 +1 - - - 120 120
22 -1 0 +1 - - - 120 120
23 0 0 +1 - - - 120 120
24 +1 0 +1 - - - 120 120
25 -1 +1 +1 - - - 120 120
26 0 +1 +1 - - - 120 120
27 +1 +1 +1 - - - 120 120
121
6.2. Plans d’expériences
122
Tableau 6.5 – Conditions opératoires pour le deuxième plan d’expériences
Expérience H2 SO4 (%) H3 P O4 (%) T (°C) Niveau M gO Niveau Al2 O3 Niveau SiO2 Temps (min) Agitation (rpm)
1 2.0 18.0 80.0 -1 -1 -1 120 120
2 2.0 18.0 80.0 0 -1 -1 120 120
3 2.0 18.0 80.0 +1 -1 -1 120 120
4 2.0 18.0 80.0 -1 0 -1 120 120
5 2.0 18.0 80.0 0 0 -1 120 120
6 2.0 18.0 80.0 +1 0 -1 120 120
7 2.0 18.0 80.0 -1 +1 -1 120 120
8 2.0 18.0 80.0 0 +1 -1 120 120
9 2.0 18.0 80.0 +1 +1 -1 120 120
10 2.0 18.0 80.0 -1 -1 0 120 120
11 2.0 18.0 80.0 0 -1 0 120 120
12 2.0 18.0 80.0 +1 -1 0 120 120
13 2.0 18.0 80.0 -1 0 0 120 120
14 2.0 18.0 80.0 0 0 0 120 120
15 2.0 18.0 80.0 +1 0 0 120 120
16 2.0 18.0 80.0 -1 +1 0 120 120
17 2.0 18.0 80.0 0 +1 0 120 120
18 2.0 18.0 80.0 +1 +1 0 120 120
19 2.0 18.0 80.0 -1 -1 +1 120 120
20 2.0 18.0 80.0 0 -1 +1 120 120
21 2.0 18.0 80.0 +1 -1 +1 120 120
Chapitre 6. Validation expérimentale du modèle thermodynamique
123
Chapitre 6. Validation expérimentale du modèle thermodynamique
L’attaque du phosphate tricalcique par les acides sulfurique et phosphorique produit une
phase liquide (acide phosphorique) et une phase solide (sulfates de calcium). Ces phases
subissent plusieurs analyses chimiques et physiques. L’ensemble des analyses effectuées
sont résumées dans le tableau 6.6 ci-dessous.
124
6.5. Validation du modèle développé
(Bichri et al., 2020 ; Bichri et al., 2022 ; Khamar, 2012 ; Khamar et al., 2014). Toute-
fois, ces valeurs sont légèrement faibles vue l’utilisation de produits purs qui augmentent
le rendement de la réaction et minimisent les pertes en phosphate. Quant à la concentra-
tion en calcium, la figure 6.2E montre qu’elle varie entre 31% et 35.5% avec une moyenne
d’environ 32.5%. Ces concentrations sont aussi en bon accord avec les analyses effectuées
par (Bichri et al., 2020 ; Khamar, 2012 ; Zmemla et al., 2020). Enfin, la concentration
des sulfates libres en phase solide est présentée sur la figure 6.2F, et montre que la quan-
tité de sulfates est élevée. La comparaison des valeurs de la concentration de cet élément
avec les travaux de (Bichri et al., 2020 ; Khamar et al., 2014) montre que les résultats
des analyses sont cohérents.
Sijm − Sijm kije − kije
P
i,j
r = r (6.1)
2 P 2
Sijm − Sijm e e
P
i,j i,j kij − kij
v
uX Sijm − kije 2
u
RM SE = t (6.2)
i,j
N
où Sijm est la prédiction du modèle, kije est la valeur expérimentale correspondante, Sijm et
kije les moyennes des prédictions et des mesures. Les indices i et j se réfèrent aux concen-
trations et aux températures respectivement.
Les figures 6.3A-C montrent les prédictions du modèle comparées à des mesures de pH à
différentes concentrations et températures. Comme attendu, les résultats montrent que les
valeurs du pH diminuent avec l’augmentation de la concentration des acides sulfurique et
phosphorique. Par ailleurs, les résultats montrent également que le pH augmente avec la
température. Ceci est cohérent avec les résultats de Khamar (2012) où il a été montré
que le pH des systèmes très similaires varient entre 0.3 et 0.85, pour des températures
125
Chapitre 6. Validation expérimentale du modèle thermodynamique
allant de 60°C à 80°C. Le coefficient de Pearson r pour ces cas varie entre 0.64 et 0.94,
alors que les RM SE varient entre 0.03 et 0.09. Ces valeurs montrent que l’accord entre
les prédictions et les mesures est satisfaisant.
Les figures 6.3D-F présentent une comparaison entre les prédictions du modèle et des
mesures de densités à différentes concentrations et températures. Les résultats montrent
que la densité augmente avec la concentration en acide phosphorique puisque la quantité
d’eau diminue et que la masse volumique de l’acide est supérieure à celle de l’eau. Par
ailleurs, les concentrations d’acide sulfurique utilisées n’ont pas une influence significative
sur cette grandeur vue que les quantités d’excès utilisées ne sont pas très grandes. Par
contre, la densité augmente en diminuant la température du fait que les molécules de la
phase liquide se rapprochent (Lide et al., 2005). Les valeurs de r varient entre 0.51 et
0.97, alors que celles de RM SE varient entre 0.004 et 0.02. Elles reflètent un bon accord
entre les prédictions et les mesures.
Les figures 6.4A-C présentent une comparaison des prédictions du modèle avec des expé-
riences de conductivité effectuées à différentes concentrations et températures. Les résul-
tats montrent que la conductivité augmente avec l’ajout d’acide. En effet, la spéciation
des acides en phase aqueuse favorisent le milieu réactionnel en ions chargés. De plus, la
conductivité est plus sensible à la concentration de l’acide sulfurique qu’à celle de l’acide
phosphorique. Dans (Lide et al., 2005 ; Wolf, 1966), il a été montré que la conductivité
spécifique des ions à base de souffre est très importante par rapport à ceux à base de
phosphate. Les indices de performances sont calculés et confirment toujours le bon accord
entre les mesures et les prédictions.
L’influence des impuretés est également évaluée. En effet, les figures 6.4D-E montrent
la variation du pH de la phase liquide en fonction des concentrations en magnésium, en
aluminium et en silice. Les concentrations en soufre et en phosphate et la température
sont fixées dans ces expériences. Il est clair que l’effet du dopage du phosphate tricalcique
avec du magnésium, de l’aluminium et de la silice sur le pH de la phase liquide produite
est négligeable. Ceci est logique vu que le pH dépend surtout de l’activité de l’ion H +
qui n’est pas favorisé par ces éléments ajoutés. Le coefficient r est légèrement faible en
raison de points aberrantes, et les faibles valeurs de la RM SE montrent toujours la bonne
performance du modèle. Par ailleurs, la densité est plus influencée par l’ajout de ces impu-
retés, elle augmente avec l’augmentation de leurs concentrations (figures 6.4F et 6.5A-B).
(Budz et al., 2007 ; Feldmann et George, 2013 ; Guan et al., 2010 ; Khamar, 2012)
ont tous fait les mêmes constats que ce soit par expérimentation ou par modélisation
thermodynamique.
126
6.5. Validation du modèle développé
127
Chapitre 6. Validation expérimentale du modèle thermodynamique
Figure 6.3 – Comparaison des prédictions du modèle développé avec les mesures
128
6.5. Validation du modèle développé
Figure 6.4 – Comparaison des prédictions du modèle développé avec les mesures
129
Chapitre 6. Validation expérimentale du modèle thermodynamique
Figure 6.5 – Comparaison des prédictions du modèle développé avec les mesures
Les comparaisons entre les prédictions et les mesures expérimentales mettent en évidence
les bonnes performances du modèle développé. Cependant, une validation statistique per-
mettrait de confirmer ces performances. Plusieurs tests statistiques peuvent être utilisés
pour la validation. On peut notamment citer le test T de Student, le test F de Fisher-
Snedecor, ou encore le test T de Hotelling, qui sont les plus utilisés dans la littérature
(Dusart, 2012). Ces tests imposent que les mesures et les prédictions suivent une loi
normale. Dans les paragraphes qui suivent, la normalité des mesures et des prédictions
est analysée à l’aide du test de normalité de Kolmogorov-Smirnov, puis, le test F de
Fisher-Snedecor est utilisé pour la validation du modèle.
130
6.5. Validation du modèle développé
la figure 6.6, où on présente les échantillons centrés réduits 2 en fonction des distributions
cumulées. Les résultats montrent que pour tous les échantillons étudiés Dn < Dn∗ . En
conclusion, l’hypothèse nulle H0 ne peut pas être rejetée, ce qui assure la normalité des
mesures et des prédictions utilisées, et permet l’utilisation du test F pour la validation
statistique.
2. Une variable centrée réduite est la transformée d’une variable aléatoire par une application affine,
de telle sorte que sa moyenne soit nulle et son écart type égal à un.
131
Chapitre 6. Validation expérimentale du modèle thermodynamique
S̃n21
Z̃ = (6.4)
S̃n22
avec :
n1
1 X
S̃n21 = (Xi − m1 )2 (6.5)
n1 i=1
n2
1 X
S̃n22 = (Yi − m2 )2 (6.6)
n2 i=1
Le tableau 6.7 résume les résultats du test de Fisher-Snedecor pour la validation statistique
du modèle développé. En effet, il contient les moyennes, les écarts-types et les variances
des sorties du modèle ainsi que des mesures expérimentales. Ces paramètres sont utilisés
pour le calcul de la statistique Z̃ du test. De plus, afin d’évaluer les performances du
modèle, les coefficients de Pearson ainsi que les RM SE sont recalculés pour l’ensemble
des échantillons. Leurs valeurs permettent d’assurer un bon accord entre les prédictions
du modèle et les expériences. Enfin, le tableau contient les valeurs de la statistique Z̃
et ses valeurs critiques lues sur la table de Fisher avec un niveau de risque de 5%. Pour
l’ensemble des cas étudiés, la statistique Z̃ a des valeurs inférieures aux valeurs critiques
Z ∗ , ce qui ne permet pas de rejeter l’hypothèse nulle H0 et de valider par conséquent les
prédictions du modèle développé par rapport aux mesures expérimentales.
132
Tableau 6.7 – Validation du modèle developpé par le test F de Fisher-Snedecor.
Sortie du modèle pH (exp) pH (mod) Densité (exp) Densité (mod) Conductivité (exp) Conductivité (mod)
133
6.5. Validation du modèle développé
Chapitre 6. Validation expérimentale du modèle thermodynamique
Expérience Excès massique (S) Titre massique (P) Température (°C) Indice de saturation
4 1% 18% 70 1.76
7 1% 20% 70 1.57
9 2% 20% 70 1.71
10 3% 16% 76 1.94
13 1% 18% 76 1.74
17 2% 20% 76 1.68
20 1% 18% 82 1.72
24 3% 18% 82 1.73
26 2% 20% 82 1.66
134
6.5. Validation du modèle développé
Figure 6.7 – Résultats de diffraction aux rayons X des échantillons solides. (A) : expé-
rience 4. (B) : expérience 7. (C) : expérience 9.
135
Chapitre 6. Validation expérimentale du modèle thermodynamique
Figure 6.8 – Résultats de diffraction aux rayons X des échantillons solides. (A) : expé-
rience 10. (B) : expérience 13. (C) : expérience 17.
136
6.5. Validation du modèle développé
Figure 6.9 – Résultats de diffraction aux rayons X des échantillons solides. (A) : expé-
rience 20. (B) : expérience 24. (C) : expérience 26.
137
Chapitre 6. Validation expérimentale du modèle thermodynamique
Les calculs thermodynamiques sont basés sur l’utilisation du modèle de Pitzer, et sur la
résolution des équations de bilans de matière, de charge et d’enthalpie, et des équations
d’équilibres chimiques. La résolution fait appel à une base de données thermodynamiques
développée, dont les performances sont testées sur une centaine de mesures expérimen-
tales, effectuées sur une dizaine de systèmes thermodynamiques binaires et ternaires, mis
en jeu lors de l’attaque du minerai de phosphate par les acides sulfuriques et phospho-
riques.
L’utilisation de cette interface est relativement aisée à condition de bien comprendre les
termes, les concepts et les méthodes employés. La figure 6.10 présente une vue d’en-
semble des caractéristiques de l’interface graphique à renseigner pour la réalisation des
simulations. La zone (A) concerne les conditions opératoires qui doivent être fournies par
l’utilisateur, et qui correspondent au fonctionnement du procédé. Elles sont constituées
principalement du profil de phosphate mis en jeu, qui doit être caractérisé par sa composi-
tion chimique en éléments (phosphate, calcium, fluore, silice, magnésium, aluminium, fer,
sodium, potassium, chlore), de la concentration d’acide sulfurique utilisé pour l’attaque
du minerai, de la concentration d’acide phosphorique de retour utilisé pour la dissolution,
et finalement de la température de la réaction. Pour ces conditions opératoires, la com-
position chimique doit être en (moles/kg d’eau), et la température en (°C). La zone (B)
affiche les propriétés thermodynamiques de la phase liquide à l’équilibre, ainsi que sa com-
position chimique en (moles/kg d’eau) pour les espèces ioniques ou moléculaires les plus
significatives. La zone (C) est consacrée aux résultats et doit fournir un graphe en barre
avec les indices de saturation des phases minérales qui ont un potentiel de précipiter lors
de la réaction d’attaque. Pour des raisons d’organisation, seuls les indices de saturation
positive s’affichent pour montrer clairement les phases minérales qui seront précipitées.
138
6.7. Conclusion
Il est finalement à noter que la case sous «Thermodynamic Database » permet de spécifier
l’utilisation de la base de données optimisée (développée dans le cadre de ce travail). Le
bouton «Run the simulation» est mis en place pour lancer les calculs après avoir fourni
toutes les conditions opératoires nécessaires. Le bouton «Clear all» sert à initialiser le
simulateur et lancer de nouvelles simulations.
6.7 Conclusion
Dans ce chapitre, une expérimentation de la fabrication de l’acide phosphorique par voie
humide est effectuée. L’objectif est de comparer les prédictions du modèle développé dans
le chapitre 2 avec des expériences qui n’ont pas servi à son identification. Deux plans
d’expériences complets ont été alors mis au point. Le premier plan met en jeu trois fac-
teurs (c’est-à-dire l’excès de l’acide sulfurique par rapport à sa quantité stoechiométrique,
le pourcentage massique de l’acide phosphorique et la température), avec trois niveaux
différents. Quant au second, il a pour objectif de comparer les prédictions du modèle
aux résultats des expériences en présence d’impuretés. Plusieurs analyses chimiques ont
été réalisées par rapport aux deux phases produites (liquide et solide). Ces analyses ont
d’abord servi pour confirmer la cohérence des mesures par rapport à des travaux de lit-
térature, et ont été ensuite comparées aux prédictions du modèle. En effet, la validation
du modèle développé s’est basée (i) sur le calcul des indices de performance, à savoir le
coefficient de Pearson et la racine carrée de l’erreur quadratique moyenne, (ii) et sur l’uti-
lisation du test statistique F de Fisher. La normalité des échantillons a été vérifiée à l’aide
du test de Kolmogorov-Smirnov. Les deux techniques de validation, ont permis de valider
le modèle développé par rapport aux mesures expérimentales. Enfin, les prédictions du
potentiel de précipitation du gypse par le modèle sont validées par les analyses DRX des
phases solides produites.
Il convient cependant de souligner que bien que le modèle développé soit validé, il peut
toujours être amélioré, notamment par la prise en compte de plus de mesures expérimen-
tales et d’autres propriétés thermodynamiques, lors de l’optimisation de ses paramètres
inconnus.
139
140
Chapitre 6. Validation expérimentale du modèle thermodynamique
141
6.7. Conclusion
Chapitre 6. Validation expérimentale du modèle thermodynamique
142
Conclusion générale et perspectives
Une étude bibliographique sur les grandeurs thermodynamiques les plus importantes, ainsi
que les modèles électrolytiques les plus utilisés a été réalisée. L’objectif de cette partie
est de faire le lien entre les grandeurs thermodynamiques et la modélisation des solutions
étudiées, et de décrire les modèles électrolytiques afin de choisir le modèle thermodyna-
mique le plus approprié pour la modélisation de ces systèmes. Cette étude a permis de
sélectionner le modèle de Pitzer (2018) comme modèle pour nos cas d’études.
143
Conclusion générale et perspectives
été utilisées pour comparer les prédictions du modèle aux mesures expérimentales. Les
résultats ont montré un bon accord entre les mesures utilisées et les prédictions, et ceci
a été confirmé à l’aide du calcul du coefficient de Pearson r dont les valeurs sont très
proches de 1.
Le modèle ainsi validé a d’abord été exploité pour l’optimisation d’un procédé de fa-
brication d’acide phosphorique. L’objectif était de déterminer les conditions opératoires
(c’est-à-dire, le débit d’acide sulfurique et ses taux de distribution sur les sections de
la cuve d’attaque, ainsi que la chaleur à évacuer afin de contrôler la température de la
réaction) qui minimisent les pertes en phosphate et améliorent le rendement de la réac-
tion. Un modèle d’optimisation multicritère a été formulé à partir des bilans de matière,
de charge et de chaleurs, et d’équations d’équilibres chimiques. Le modèle thermodyna-
mique développé a servi pour prédire deux fonctions objectifs (conflictuels) du problème
d’optimisation. La méthode -contrainte a été utilisée pour transformer le problème d’opti-
misation multicritère en plusieurs problèmes d’optimisation mono-critère, et sa résolution
est effectuée à l’aide d’un algorithme "Branch-and-Reduce" sous l’environnement GAMS.
Les résultats d’optimisations (Front de Pareto) ont été classées en utilisant la méthode
MAUT d’aide à la décision. Une solution optimale a été alors sélectionnée et comparée
à des mesures réalisées au niveau industriel. La comparaison a montré que la solution
obtenue permettait de réduire les pertes syncristallisées jusqu’à des valeurs de 0.6%, et les
pertes inattaquées de 0.12%, ce qui correspond à environ 54 tonnes de phosphate économi-
sées par jour (cas d’une usine à Jorf Lasfar au Maroc qui contient 16 unités de production
à 1400 tonnes de cadence par jour chacune). Cette comparaison a permis de montrer
que l’approche d’optimisation muticritère est très prometteuse, et permet d’obtenir des
résultats très proches de ceux obtenus industriellement, et qui demandaient notamment
d’investir beaucoup de main d’œuvre, de consommer beaucoup de réactifs et de temps, et
surtout pouvaient engendrer des problèmes liés à la sécurité.
144
Une deuxième exploitation du modèle développé concerne la prédiction du potentiel de
plusieurs minéraux à précipiter pendant la digestion du minerai de phosphate. Ces miné-
raux précipitent à cause de réactions secondaires et conduisent à des phénomènes d’encras-
sement. Il a été montré dans cette partie que dans les conditions actuelles de production,
la malladrite, la hiératite, la magadiite, le quartz, et l’hexafluorosilicate de magnésium,
sous leurs différentes formes cristallines, sont les minéraux les plus susceptibles de précipi-
ter et de provoquer des phénomènes d’encrassement. Par ailleurs, une étude de sensibilité
globale a permis d’étudier la sensibilité des indices de saturation (qui quantifient le poten-
tiel de précipitation) de ces minéraux en fonction des conditions opératoires du procédé.
Il a été montré que pour limiter la précipitation de certains de ces minéraux, la concen-
tration d’acide sulfurique et la température sont à optimiser. De plus, le traitement de
l’eau de procédé utilisée serait une solution efficace qui devrait être mise en œuvre par les
industriels pour minimiser ces problèmes d’encrassement.
Bien que les résultats obtenus soient promoteurs, le modèle développé peut être davantage
amélioré. En effet, sa précision et sa fiabilité peuvent être améliorées de plusieurs façons.
D’une part, l’ajout de plus de mesures d’équilibre permettrait d’augmenter le nombre
de paramètres estimables, et par conséquent d’identifier ceux qui ont été fixés à partir
de la littérature et/ou de travaux antérieurs. Une autre amélioration importante pour la
145
Conclusion générale et perspectives
précision consiste à utiliser la méthode de calibration bayésienne (Smith, 2013) qui est
fondamentalement différente des méthodes classiques de calibration, où la différence entre
les données observées et les sorties du modèle est minimisée. La calibration bayésienne est
une méthode itérative où les distributions d’incertitude des paramètres du modèle sont
mises à jour de manière cohérente avec les données observées.
146
Annexe A : Valeurs des paramètres
inconnus du modèle de Pitzer
(0) 3
β + 4.211 [3.75 ;4.67] 0.134 [0.12 ;0.15] 2.49.10 [2.21.103 ;2.76.103 ]
H ,H3 P O4
(0)
β − -2.965 [-3.32 ;-2.61] -1.005 [-1.13 ;-0.88] 203.473 [179.06 ;227.89]
H + ,H2 P O4
(0)
β − -186.465 [-242.4 ;-130.5] 2.203 [1.54 ;2.86] 0.198 [0.14 ;0.26]
H + ,H5 P2 O8
(0)
β − 0.183 [0.17 ;0.20] -0.088 [-0.10 ;-0.08] -1.25.104 [-1.37.104 ;-1.14.104 ]
H3 P O4 ,H2 P O4
(0)
β − 4.739 [4.17 ;5.31] -0.161 [-0.18 ;-0.14] 0.925 [0.81 ;1.04]
H3 P O4 ,H5 P2 O8
(0) 5
β − − 168.178 [126.1 ;210.2] 1.941 [1.46 ;2.43] 4.49.10 [3.36.105 ;5.61.105 ]
H2 P O4 ,H5 P2 O8
(0)
β − − -4.093 [-4.91 ;-3.27] 0.236 [0.19 ;0.28] - -
H2 P O4 ,H5 P2 O8
(0)
βH P O ,H P O - - - - -134.046 [-135.39 ;-132.71]
3 4 3 4
(1)
β + -2.495 [-3.24 ;-1.75] -0.311 [-0.40 ;-0.22] -4.96.104 [-6.45.104 ;-3.47.104 ]
H ,H3 P O4
(1)
β − 16.109 [12.4 ;19.8] 1.404 [1.08 ;1.73] - -
H + ,H2 P O4
(1)
β − 56.786 [45.4 ;68.1] - - 9.14.104 [7.31.104 ;1.09.104 ]
H + ,H5 P2 O8
(1)
β − -8.458 [-9.30 ;-7.61] 0.064 [0.06 ;0.07] -1871.020 [-2058.12 ;-1683.92]
H3 P O4 ,H2 P O4
(1) 4
β − -3.712 [-4.19 ;-3.23] - - -3.05.10 [-3.44.104 ;-2.65.104 ]
H3 P O4 ,H5 P2 O8
ψH3 P O4 ,H3 P O4 ,H3 P O4 6.826 [5.46 ;8.19] - - - -
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
(0)
β − 2− -2.471 [-3.00 ;-1.942] 0.007 [0.004 ;0.010] 0.01 -
HSO4 ,SO4
(0)
β − 1.046 [0.518 ;1.573] -0.002 [-0.003 ;-0.001] 0.01 -
HSO4 ,H +
(0)
β 2− -8.001 [-9.665 ;-6.336] 0.01 - 1684.0 [1221.1 ;2147.04]
SO4 ,H +
(1)
β − 2− -254.8 [-312.8 ;-196.8] 0.405 [0.341 ;0.468] 3.9.104 [3.1.104 ;4.8.104 ]
HSO4 ,SO4
(1)
β − -296.9 [-489.1 ;-104.8] 0.453 [0.341 ;0.564] 4.8.104 [2.7.104 ;7.0.104 ]
HSO4 ,H +
(1)
β 2− 8.039 [4.879 ;11.19] -0.023 [-0.032 ;-0.013] 0.01 -
SO4 ,H +
ψ + − 2− 1.696 [1.142 ;2.249] -0.003 [-0.004 ;-0.001] -277.209 [-335.9 ;-218.4]
H ,HSO4 ,SO4
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
147
Annexe A : Valeurs des paramètres inconnus du modèle de Pitzer
(0)
β + - - -0.007 [-0.0078 ;-0.0062] - -
H ,HF
(0)
β − 0.208 [0.19 ;0.23] - - - -
HF2 ,HF
(0)
β − - - 0.013 [0.0118 ;0.0142] - -
F ,HF
(0)
β − -0.392 [-0.44 ;-0.35] . - - -
HF2 ,HF
(0)
βHF,HF - - - - -84.628 [-93.94 ;-75.32]
ψHF,HF,HF -3.081 [-4.01 ;-2.16] - - - -
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
(0)
β + - - 0.072 [0.065 ;0.079] - -
H ,H2 SiF6
(0)
β 2− - - -1.338 [-1.46 ;-1.21] -9.071.104 [-1.00.105 ;-8.07.104 ]
H + ,SiF6
(0)
β 2− 12.05 [10.7 ;13.4] 0.077 [0.067 ;0.086] -3.746 [-4.08 ;-3.40]
H2 SiF6 ,SiF6
(0)
βH SiF ,H SiF - - 0.005 [0.004 ;0.0055] - -
2 6 2 6
(1)
β + -4.595 [-5.01 ;-4.18] -7.095 [-8.15 ;-6.030] -7.728.105 [-8.65.105 ;-6.80.105 ]
H ,H2 SiF6
(1)
β + 4.901 [4.31 ;5.48] - - -9.676.104 [-18.2.104 ;-1.11.104 ]
H ,H2 SiF6
(1) 4
β 2− -63.290 [-72.7 ;-53.7] - - -8.341.10 [-9.42.104 ;-7.25.104 ]
H + ,SiF6
(1)
β 2− -0.149 [-0.16 ;-0.12] - - - -
H2 SiF6 ,SiF6
(1) −6 −6
βH SiF ,H SiF 1.395.10 [0.01 ;2.78].10 0.040 [0.036 ;0.0436] - -
2 6 2 6
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
(0) −4
β 2− 0.205 [0.124 ;0.286] -951.0 - -2.34.10 -
M g 2+ ,SO4
(1) 3 −1
β 2− 3.354 [2.102 ;4.606] -5.78.10 - -1.48.10 -
Mg 2+ ,SO4
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
(0)
β − 0.351 - - - -9.32.10−4 -
Cl ,M g 2+
(1)
β − 1.65 - -2.78.10−7 [-3.55 ;-2.01].10−7 -1.09.10−2 -
Cl ,M g 2+
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
148
Tableau 7 – Paramètres d’interactions de Pitzer pour le système CaSO4 − H2 O − H2 O.
a C.I b C.I c C.I
(0)
β − 2− -2.471 [-3.00 ;-1.934] 0.007 [0.004 ;0.010] 0.01 -
HSO4 ,SO4
(0)
β − 1.046 [0.518 ;1.573] -0.002 [-0.003 ;-0.001] 0.01 -
HSO4 ,H +
(0)
β 2− -8.001 [-9.665 ;-6.336] 0.01 - 1684.0 [1221.1 ;2147.04]
SO4 ,H +
(0)
β − 0.920 [0.81 ;1.02] -0.030 [-0.033 ;-0.026] -3272.216 [-3632.1 ;-2912.2]
HSO4 ,CaSO4
(0)
β − 1.010 [0.91 ;1.11] 0.026 [0.023 ;0.028] - -
HSO4 ,Ca2+
(0)
β 2− -1.801 [-2.01 ;-1.58] 0.010 [0.0088 ;0.0112] - -
SO4 ,CaSO4
(0)
β 2− -1.037 [-1.19 ;-0.88] 0.066 [0.0561 ;0.0759] -1075.927 [-1194.2 ;-957.5]
SO4 ,Ca2+
(0)
β + 0.697 [0.60 ;0.78] 0.033 [0.028 ;0.037] - -
H ,CaSO4
(1)
β − 2− -254.8 [-312.8 ;-196.8] 0.405 [0.341 ;0.468] 3.9.104 [3.1.104 ;4.8.104 ]
HSO4 ,SO4
(1)
β − -296.9 [-489.1 ;-104.8] 0.453 [0.341 ;0.564] 4.8.104 [2.7.104 ;7.0.104 ]
HSO4 ,H +
(1)
β 2− 8.039 [4.879 ;11.19] -0.023 [-0.032 ;-0.013] 0.01 -
SO4 ,H +
(1)
β + -7.917 [-8.78 ;-7.04] - - - -
H ,CaSO4
(1)
β + - - - - -912.578 [-1022.0 ;-803.0]
H ,Ca2+
(1)
β -25.563 [-28.8 ;-23.1] - - - -
CaSO4 ,Ca2+
(1)
β − 3.43 [1.211 ;5.656] - - -0.01.10−3 -
Ca2+ ,HSO4
(1)
β 2− 3.546 - - - 5.77.10−3 -
Ca2+ ,SO4
ψ + − 2− 1.696 [1.142 ;2.249] -0.003 [-0.004 ;-0.001] -277.209 [-335.9 ;-218.4]
H ,HSO4 ,SO4
ψCaSO4 ,CaSO4 ,CaSO4 0.763 [0.69 ;0.83] - - 1647.465 [1400.3 ;1894.5]
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
(0)
β 3+ 0.612 - 6.82134.10−4 - 1.315.10−8 [0.151 ;2.479].10−8
Al ,Cl−
(1)
β 3+ 1.642 - 2.697.10−2 - 5.201.10−8 [2.812 ;9.841].10−8
Al ,Cl−
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
(0) −8 −8 −4
β 0.315 - 3.01.10 [2.57 ;5.03].10 -3.27.10 -
Ca+2 ,Cl−
(1)
β 1.614 - - - 7.63.10−3 -
Ca2+ ,Cl−
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
149
Annexe A : Valeurs des paramètres inconnus du modèle de Pitzer
(0)
β − 2− -2.471 [-3.00 ;-1.934] 0.007 [0.004 ;0.010] 0.01 -
HSO4 ,SO4
(0)
β − 1.046 [0.518 ;1.573] -0.002 [-0.003 ;-0.001] 0.01 -
HSO4 ,H +
(0)
β 2− -8.001 [-9.665 ;-6.336] 0.01 - 1684.0 [1221.1 ;2147.04]
SO4 ,H +
(0)
β 2− 0.205 [0.124 ;0.286] -951.0 - -2.34.10−4 -
M g 2+ ,SO4
(0)
β − -0.0003 - -3.17.10−7 [-2.77.10−7 ;-3.57.10−7 ] 6.65.10−7 [4.47.10−7 ;8.83.10−7 ]
HSO4 ,K +
(0)
β − 0.474 [0.110 ;0.546] - - - -
HSO4 ,M g 2+
(0)
β 2− 3.17.10−2 - - - 9.28.10−4 -
K + ,SO4
(1) 4
β − 2− -254.8 [-312.8 ;-196.8] 0.405 [0.341 ;0.468] 3.9.10 [3.1.104 ;4.8.104 ]
HSO4 ,SO4
(1)
β − -296.9 [-489.1 ;-104.8] 0.453 [0.341 ;0.564] 4.8.104 [2.7.104 ;7.0.104 ]
HSO4 ,H +
(1)
β 2− 8.039 [4.879 ;11.19] -0.023 [-0.032 ;-0.013] 0.01 -
SO4 ,H +
(1) 3 −1
β 2− 3.354 [2.102 ;4.606] -5.78.10 - -1.48.10 -
M g 2+ ,SO4
(1)
β − 0.173 - - - - -
HSO4 ,K +
(1)
β − 1.729 - - - - -
HSO4 ,M g 2+
(1) 4
β 2− 0.756 - -1.514.10 - 0.1091 -
+
K ,SO4
ψ + − 2− 1.696 [1.142 ;2.249] -0.003 [-0.004 ;-0.002] -277.209 [-335.9 ;-218.4]
H ,HSO4 ,SO4
ψ + − -0.0265 [-0.046 ;-0.0174] - - - -
H ,HSO4 ,K +
ψ + − -0.0178 - 2.13.10−7 [0.07.10−7 ;4.19.10−7 ] - -
H ,HSO4 ,M g 2+
ψ + + −2 0.197 [0.068 ;0.326] - - - -
H ,K ,SO4
ψ − −2 -0.0677 - - - - -
HSO4 ,K + ,SO4
ψ − −2 -0.042 - 5.44.10−8 [3.11.10−8 ;7.77.10−8 ] - -
HSO4 ,M g 2+ ,SO4
ψ + 2+ −2 -0.048 - - - - -
K ,M g ,SO4
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
Tableau 11 – Paramètres d’interactions de Pitzer pour le système F eSO4 −H2 SO4 −H2 O.
a C.I b C.I c C.I
(0)
β − 2− -2.471 [-3.00 ;-1.934] 0.007 [0.004 ;0.010] 0.01 -
HSO4 ,SO4
(0)
β − 1.046 [0.518 ;1.573] -0.002 [-0.003 ;-0.001] 0.01 -
HSO4 ,H +
(0)
β 2− -8.001 [-9.665 ;-6.336] 0.01 - 1684.0 [1221.1 ;2147.04]
SO4 ,H +
(0)
β − 0.427 [0.354 ;0.500] - - -1.05.10−8 [-2.07.10−8 ;-0.03.10−8 ]
F e2+ ,HSO4
(0)
β 2− 0.256 - - - - -
F e2+ ,SO4
(1) 4
β − 2− -254.8 [-312.8 ;-196.8] 0.405 [0.341 ;0.468] 3.9.10 [3.1.104 ;4.8.104 ]
HSO4 ,SO4
(1)
β − -296.9 [-489.1 ;-104.8] 0.453 [0.341 ;0.564] 4.8.104 [2.7.104 ;7.0.104 ]
HSO4 ,H +
(1)
β 2− 8.039 [4.879 ;11.19] -0.023 [-0.032 ;-0.013] 0.01 -
SO4 ,H +
(1)
β − 3.48 [1.716 ;5.784] - - - -
F e2+ ,HSO4
(1)
β 2− 3.063 - - - - -
F e2+ ,SO4
ψ + − 2− 1.696 [1.142 ;2.249] -0.003 [-0.004 ;-0.002] -277.209 [-335.9 ;-218.4]
H ,HSO4 ,SO4
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
150
Tableau 12 – Paramètres d’interactions de Pitzer pour le système N aSO4 − N aCl − H2 O.
a C.I b C.I c C.I
(0)
β − 2− -2.471 [-3.00 ;-1.934] 0.007 [0.004 ;0.010] 0.01 -
HSO4 ,SO4
(0)
β − 1.046 [0.518 ;1.573] -0.002 [-0.003 ;-0.001] 0.01 -
HSO4 ,H +
(0)
β 2− -8.001 [-9.665 ;-6.336] 0.01 - 1684.0 [1221.1 ;2147.04]
SO4 ,H +
(0) −8 −8 −4
β − + 0.177 [0.122 ;0.232] -1.39.10 [-2.77 ;-0.02].10 -3.08.10 -
Cl ,H
(0)
β − 7.53.10−2 - 9598.4 - -5.87.10−2 -
Cl ,N a+
(0)
β − 0.0454 - - - - -
HSO4 ,N a+
(0)
β 2− 2.73.10−2 [1.03 ;4.43].10−2 - - 9.89.10−4 [0.796 ;1.18].10−3
N a+ ,SO4
(1) 4
β − 2− -254.8 [-312.8 ;-196.8] 0.405 [0.341 ;0.468] 3.9.10 [3.1.104 ;4.8.104 ]
HSO4 ,SO4
(1)
β − -296.9 [-489.1 ;-104.8] 0.453 [0.341 ;0.564] 4.8.104 [2.7.104 ;7.0.104 ]
HSO4 ,H +
(1)
β 2− 8.039 [4.879 ;11.19] -0.023 [-0.032 ;-0.013] 0.01 -
SO4 ,H +
(1) −8 −8 −4
β − + 0.294 - 1.17.10 [0.01 ;2.33].10 1.419.10 -
Cl ,H
(1)
β − 0.276 [0.165 ;0.387] 1.377.104 - -6.9512.10−2 -
Cl ,N a+
(1)
β − 0.398 [0.189 ;0.607] - - - -
HSO4 ,N a+
(1)
β 2− 0.956 - 2.663.103 - 1.158.10−2 -
N a+ ,SO4
ψ + − 2− 1.696 [1.142 ;2.249] -0.003 [-0.004 ;-0.002] -277.209 [-335.9 ;-218.4]
H ,HSO4 ,SO4
ψCl− ,H + ,N a+ -0.004 [-0.008 ;-0.001] - - - -
ψCl− ,HSO4− ,H + 0.0130 - 1.17.10−8 [0.01 ;2.33].10−8 - -
ψCl− ,HSO4− ,N a+ -0.006 - - - - -
ψ − + −2 - - - - - -
Cl ,N a ,SO4
ψ + − -0.0129 - 0.17.10−8 [0.012 ;0.328].10−8 - -
H ,HSO4 ,N a+
ψ − + −2 -0.0094 - - - - -
HSO4 ,N a ,SO4
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
151
Annexe A : Valeurs des paramètres inconnus du modèle de Pitzer
(0)
β − 2− -2.471 [-3.00 ;-1.934] 0.007 [0.004 ;0.010] 0.01 -
HSO4 ,SO4
(0)
β − 1.046 [0.518 ;1.573] -0.002 [-0.003 ;-0.001] 0.01 -
HSO4 ,H +
(0)
β 2− -8.001 [-9.665 ;-6.336] 0.01 - 1684.0 [1221.1 ;2147.04]
SO4 ,H +
(0)
β 2− 0.205 [0.124 ;0.286] -951.0 - -2.34.10−4 -
M g 2+ ,SO4
(0)
β − 0.474 [0.110 ;0.838] - - - -
HSO4 ,M g 2+
(0)
β − 0.351 - - - -9.32.10−4 -
Cl ,M g 2+
(0) −8 −8 −8
β − + 0.177 [0.122 ;0.232] -1.395.10 [-2.77.10 ;-0.02.10 ] -3.08.10−4 -
Cl ,H
(1)
β − 2− -254.8 [-312.8 ;-196.8] 0.405 [0.341 ;0.468] 3.9.104 [3.1.104 ;4.8.104 ]
HSO4 ,SO4
(1)
β − -296.9 [-489.1 ;-104.8] 0.453 [0.341 ;0.564] 4.8.104 [2.7.104 ;7.0.104 ]
HSO4 ,H +
(1)
β 2− 8.039 [4.879 ;11.19] -0.023 [-0.032 ;-0.013] 0.01 -
SO4 ,H +
(1) 3 −1
β 2− 3.354 [2.102 ;4.606] -5.78.10 - -1.48.10 -
M g 2+ ,SO4
(1)
β − 1.729 - - - - -
HSO4 ,M g 2+
(1) −7 −7 −7 −2
β − 1.65 - -2.78.10 [-3.55.10 ;-2.01.10 ] -1.09.10 -
Cl ,M g 2+
(1)
β − + 0.294 - 1.17.10−8 [0.01.10−8 ;2.33.10−8 ] 1.419.10−4 -
Cl ,H
ψ + − 2− 1.696 [1.142 ;2.249] -0.003 [-0.004 ;-0.002] -277.209 [-335.9 ;-218.4]
H ,HSO4 ,SO4
ψCl− ,H + ,M g2+ -0.011 - - - - -
ψCl− ,HSO4− ,H + 0.0130 - 1.17.10−8 [0.01 ;2.33].10−8 - -
ψ − 2+ −2 -0.008 - 32.63 - - -
Cl ,M g ,SO4
ψ + − -0.0178 - 2.13.10−7 [0.07 ;4.18].10−7 - -
H ,HSO4 ,M g 2+
−8
ψ − 2+ −2 -0.042 - 5.44.10 [3.11.10−8 ;7.77.10−8 ] - -
HSO4 ,M g ,SO4
ψH3 P O4 ,H3 P O4 ,H3 P O4 AAAAAAA AAAAAAAAA AAAA AAAAAAAAA AAAAAA [-1.00.105 ;-8.07.104 ]
152
Annexe B : Modèle d’optimisation
multicritère
— Fonctions objectifs
f2 (x) ≥ (2)
Avec :
1X j
f1 (x) = log aCa2+ (x).aHP O42− (x).a2H2 O (x).ks−1 (3)
3 j 1
1X j
f2 (x) = 2 −1
log aCa2+ (x).aSO42− (x).aH2 O (x) .ks2 (4)
3 j
— Contraintes d’égalité :
— Bilans de matière :
(S)tot = mSO42− (x) + mHSO4− (x) + mCaSO4 (x) (7)
(P )tot = mH3 P O4 (x) + mH2 P O4− (x) + mHP O42− (x) + 2.mH5 P2 O8− (x) (8)
153
Annexe B : Modèle d’optimisation multicritère
— Bilans d’électroneutralité
mH + (x) + 2.mCa2+ (x) = mHSO4− (x) + 2.mSO42− (x) + mH2 P O4− (x) (12)
+[Link] O42− (x) + mH5 P2 O8− (x) + mF − (x) + mHF2 − (x) + mSiF6− (x)
— Constantes d’équilibre
mH + (x)mSO42− (x) γH + (x).γSO42− (x)
k1 = . (13)
mHSO4− (x) γHSO4− (x)
154
— Modèle de Pitzer pour calculer γi
zi2 X X 0
ln(γi ) = f +2 mj (x)λi,j (I) + zi2 mj (x)mk (x)λj,k (I) (21)
2 j j
X
+3 mi (x)mj (x)mk (x)ψi,j,k
j,k
4Aφ I
ln 1 + I 0.5 (22)
f =−
b
1X dλi,j ni (x)
I= mi (x)zi2 , λ0i,j = , mi (x) = (23)
2 i dI nw (x)
1 1 3
Aφ = (2πN0 dw /1000) 2 (e2 DkT ) 2 (24)
3
— Bilans de chaleur
QH (x) + QAG (x) + QR+D (x) = Qcool (x) + Qout (x) + Qloss (x) (25)
— Contraintes d’inégalité
1 % ≤ Sc (x) ≤ 3 % (32)
T ≤ 355 K (33)
0 ≤ mi (x) (34)
155
Annexe B : Modèle d’optimisation multicritère
156
Bibliographie
Abutayeh, M. & Campbell, S. W. (2009). Predicting the citrate soluble loss of the
dihydrate process. Industrial & Engineering Chemistry Research, 48 (18), 8670-
8677. [Link]
Al-Fariss, T., Razik, S. A., Ghasem, N. M., Abdelaleem, F. & Elnashaie, S. S.
(1993). Computer program for the mass and heat balance calculations for the wet
process phosphoric acid (hemi-hydrate) and its applications to Saudi phosphate
rock. Fertilizer research, 35 (3), 169-175. [Link]
Al-Thyabat, S. & Zhang, P. (2015). In-line extraction of REE from Dihydrate (DH)
and HemiDihydrate (HDH) wet processes. Hydrometallurgy, 153, 30-37. https :
//[Link]/10.1016/[Link].2015.01.010
Amin, M., Ali, M., Kamal, H., Youssef, A. & Akl, M. (2010). Recovery of high grade
phosphoric acid from wet process acid by solvent extraction with aliphatic alcohols.
Hydrometallurgy, 105 (1-2), 115-119. [Link]
007
Andalibi, M. R., Kumar, A., Srinivasan, B., Bowen, P., Scrivener, K., Ludwig,
C. & Testino, A. (2018). On the mesoscale mechanism of synthetic calcium–
silicate–hydrate precipitation : a population balance modeling approach. Journal
of Materials Chemistry A, 6 (2), 363-373. [Link]
Andalibii, M. R., Bowen, P., Carino, A. & Testino, A. (2020). Global uncertainty-
sensitivity analysis on mechanistic kinetic models : From model assessment to
theory-driven design of nanoparticles. Computers & Chemical Engineering, 140,
106971. [Link]
Archer, D. G. & Rard, J. A. (1998). Isopiestic Investigation of the Osmotic and Activity
Coefficients of Aqueous MgSO4 and the Solubility of MgSO4-7H2O (cr) at 298.15
K : Thermodynamic Properties of the MgSO4+ H2O System to 440 K. Journal of
Chemical & Engineering Data, 43 (5), 791-806. [Link]
Awbery, J. & Griffiths, E. (1940). The thermal capacity of pure iron. Proceedings
of the Royal Society of London. Series A. Mathematical and Physical Sciences,
174 (956), 1-15. [Link]
Azaroual, M. [Mohamed], Kervevan, C., Lassin, A., André, L., Amalhay, M.,
Khamar, L. & Guendouzi, M. E. (2012). Thermo-kinetic and physico-chemical
modeling of processes generating scaling problems in phosphoric acid and fertilizers
production industries. Procedia Engineering, 46, 68-75. [Link]
10.1016/[Link].2012.09.447
157
BIBLIOGRAPHIE
Azimi, G. & Papangelakis, V. (2010). The solubility of gypsum and anhydrite in simu-
lated laterite pressure acid leach solutions up to 250 C. Hydrometallurgy, 102 (1-4),
1-13. [Link]
Bailey, C., Taylor, G., Cross, M. & Chow, P. (1999). Discretisation procedures for
multi-physics phenomena. Journal of Computational and Applied Mathematics,
103 (1), 3-17. [Link]
Baker, L. E., Pierce, A. C. & Luks, K. D. (1982). Gibbs energy analysis of phase
equilibria. Society of Petroleum Engineers Journal, 22 (05), 731-742. [Link]
org/10.2118/9806-PA
Bakir, T. (2006). Estimation d’un procédé de cristallisation en batch (thèse de doct.).
Université Claude Bernard-Lyon I.
Ball, J. W. & Nordstrom, D. K. (1991). WATEQ4F–User’s manual with revised ther-
modynamic data base and test cases for calculating speciation of major, trace and
redox elements in natural waters (rapp. tech.).
Baptiste, B. (2021). Choix d’un modèle thermodynamique et simulation. emse.
Barbier, E., Coste, M., Genin, A., Jung, D., Lemoine, C., Logette, S. & Muhr,
H. (2009). Simultaneous determination of nucleation and crystal growth kinetics
of gypsum. Chemical Engineering Science, 64 (2), 363-369. [Link]
1016/[Link].2008.10.036
Becker, P. (1983). Phosphates and phosphoric acid.
Becker, P. (1989). Phosphates and phosphoric acid : raw materials, technology, and
economics of the wet process. revised and expanded. (T. 6). Marcel Dekker, Inc.
Belotti, P., Lee, J., Liberti, L., Margot, F. & Wächter, A. (2009). Branching
and bounds tighteningtechniques for non-convex MINLP. Optimization Methods
& Software, 24 (4-5), 597-634. [Link]
Benyahia, B. (2009). Modélisation, expérimentation et optimisation multicritère d’un
procédé de copolymérisation en émulsion en présence d’un agent de transfert de
chaıne (thèse de doct.). Institut National Polytechnique de Lorraine.
Benyahia, B. (2018). Applications of a plant-wide dynamic model of an integrated
continuous pharmaceutical plant : design of the recycle in the case of multiple
impurities. Computer Aided Chemical Engineering (p. 141-157). Elsevier. https:
//[Link]/10.1016/B978-0-444-63963-9.00006-3
Benyahia, B. et al. (2019). Estimability analysis and parameter identification for a batch
emulsion copolymerization reactor in the presence of a chain transfer agent.
Benyahia, B., Anandan, P. D. & Rielly, C. (2021). Control of Batch and Continuous
Crystallization Processes using Reinforcement Learning. Computer Aided Chemical
Engineering (p. 1371-1376). Elsevier. [Link]
5.50211-4
Benyahia, B., Bandulasena, M. V., Bandulasena, H. H. & Vladisavljević, G. T.
(2021). Experimental and computational analysis of mixing inside droplets for
microfluidic fabrication of gold nanoparticles. Industrial & Engineering Chemistry
Research, 60 (38), 13967-13978.
Benyahia, B., Latifi, A. M., Fonteix, C. & Pla, F. (2013). Emulsion copolymerization
of styrene and butyl acrylate in the presence of a chain transfer agent. Part 2 :
158
BIBLIOGRAPHIE
159
BIBLIOGRAPHIE
process. 2020 6th IEEE Congress on Information Science and Technology (CiSt),
389-394. [Link]
Bouchkira, I., Latifi, A. M., Khamar, L. & Benjelloun, S. (2021c). Process mode-
ling and multi-criteria optimization of an industrial phosphoric acid wet-process.
50, 499-504. [Link]
Bouchkira, I., Latifi, A. M., Khamar, L. & Benjelloun, S. (2022). Modeling and
multi-objective optimization of the digestion tank of an industrial process for ma-
nufacturing phosphoric acid by wet process. Computers & Chemical Engineering,
156, 107536. [Link]
Braddy, R., McTigue, P. & Verity, B. (1994). Equilibria in moderately concentrated
aqueous hydrogen fluoride solutions. Journal of fluorine chemistry, 66 (1), 63-67.
[Link]
Bromley, L. A. (1973). Thermodynamic properties of strong electrolytes in aqueous
solutions. AIChE Journal, 19 (2), 313-320. [Link]
Budz, J. A., Jones, A. & Mullin, J. W. (2007). Effect of selected impurities on the
continuous precipitation of calcium sulphate (gypsum). Journal of Chemical Tech-
nology & Biotechnology, 36, 153-161.
Bullough, W., Canning, T. & Strawbridge, M. (1952). The solubility of ferrous
sulphate in aqueous solutions of sulphuric acid. Journal of Applied Chemistry,
2 (12), 703-707. [Link]
Cardenas, C., Farrusseng, D., Daniel, C. & Aubry, R. (2022). Modeling of equili-
brium water vapor adsorption isotherms on activated carbon, alumina and hopca-
lite. Fluid Phase Equilibria, 113520. [Link]
Cardenas, C., Marsteau, S., Sigot, L., Vallières, C. & Latifi, A. M. (2020).
Analysis of an Industrial Adsorption Process based on Ammonia Chemisorption :
Modeling and Simulation. Computer Aided Chemical Engineering (p. 625-630).
Elsevier. [Link]
Cardenasi, C., Latifi, A. M., Vallières, C., Marsteau, S. & Sigot, L. (2021).
Analysis of an industrial adsorption process based on ammonia chemisorption :
modeling and simulation, 107474. [Link]
107474
Cauvain & Young. (2009). Bakery food manufacture and quality : water control and
effects. John Wiley & Sons.
Chang, C.-K. & Lin, S.-T. (2019). Extended Pitzer–Debye–Hückel Model for Long-Range
Interactions in Ionic Liquids. Journal of Chemical & Engineering Data, 65 (3),
1019-1027. [Link]
Chen, C.-C., Bokis, C. P. & Mathias, P. (2001). Segment-based excess Gibbs energy
model for aqueous organic electrolytes. AIChE Journal, 47 (11), 2593-2602. https:
//[Link]/10.1002/aic.690471122
Chen, C.-C., Britt, H. I., Boston, J. & Evans, L. (1982). Local composition model for
excess Gibbs energy of electrolyte systems. Part I : Single solvent, single completely
dissociated electrolyte systems. AIChE Journal, 28 (4), 588-596. [Link]
10.1002/aic.690280410
160
BIBLIOGRAPHIE
Chen, C.-C. & Song, Y. (2004). Generalized electrolyte-NRTL model for mixed-solvent
electrolyte systems. AIChE Journal, 50 (8), 1928-1941. [Link]
aic.10151
Cherif, M., Mgaidi, A., Ammar, M. N., Abderrabba, M. & Fürst, W. (2000).
Modelling of the equilibrium properties of the system H3PO4–H2O : representation
of VLE and liquid phase composition. Fluid Phase Equilibria, 175 (1-2), 197-212.
[Link]
Cherif, M., Mgaidi, A., Ammar, M. N., Abderrabba, M. & Fürst, W. (2002a).
Representation of VLE and liquid phase composition with an electrolyte model :
application to H3PO4–H2O and H2SO4–H2O. Fluid phase equilibria, 194, 729-738.
[Link]
Cherif, M., Mgaidi, A., Ammar, M. N., Abderrabba, M. & Fürst, W. (2002b).
Representation of VLE and liquid phase composition with an electrolyte model :
application to H3PO4–H2O and H2SO4–H2O. Fluid phase equilibria, 194, 729-738.
[Link]
Cherif, M., Mgaidi, A., Ammar, M. N., Abderrabba, M. & Fürst, W. (2002c).
Representation of VLE and liquid phase composition with an electrolyte model :
application to H3PO4–H2O and H2SO4–H2O. Fluid phase equilibria, 194, 729-738.
[Link]
Chidambaram, S., Karmegam, U., Sasidhar, P., Prasanna, M. V., Manivannan,
R., Arunachalam, S., Manikandan, S. & Anandhan, P. (2011). Significance
of saturation index of certain clay minerals in shallow coastal groundwater, in and
around Kalpakkam, Tamil Nadu, India. Journal of earth system science, 120 (5),
897-909.
Cho, H. J., Yeo, Y.-K., Park, W. H. & Moon, B. K. (1996). Modeling and simula-
tion of a wet hemihydrate phosphoric acid process. Korean Journal of Chemical
Engineering, 13 (6), 585-595. [Link]
Chow, R., Blindt, R., Chivers, R. & Povey, M. (2005). A study on the primary
and secondary nucleation of ice by power ultrasound. Ultrasonics, 43 (4), 227-230.
[Link]
Christides, S. & Barr, A. (1984). One-dimensional theory of cracked Bernoulli-Euler
beams. International Journal of Mechanical Sciences, 26 (11-12), 639-648. https:
//[Link]/10.1016/0020-7403(84)90017-1
Christoffersen, J., Christoffersen, M. R. & Johansen, T. (1996). Kinetics of
growth and dissolution of fluorapatite. Journal of crystal growth, 163 (3), 295-303.
[Link]
Christov, C. (2000). Thermodynamic study of the Na-Cu-Cl-SO4-H2O system at the
temperature 298.15 K. The Journal of Chemical Thermodynamics, 32 (3), 285-295.
[Link]
Christov, C. (2001). Thermodynamic study of the K- Mg- Al- Cl- SO4- H2O system
at the temperature 298.15 K. Calphad, 25 (3), 445-454. [Link]
S0364-5916(01)00063-3
Christov, C., Dickson, A. G. & Moller, N. (2007). Thermodynamic Modeling of
Aqueous Aluminum Chemistry and Solid-Liquid Equilibria to High Solution Concen-
tration and Temperature. I. The Acidic H-Al-Na-K-Cl-H 2 O System from 0 to 100
161
BIBLIOGRAPHIE
162
BIBLIOGRAPHIE
El Guendouzi, M. & Marouani, M. (2003). Water activities and osmotic and acti-
vity coefficients of aqueous solutions of nitrates at 25 C by the hygrometric me-
thod. Journal of solution chemistry, 32 (6), 535-546. [Link]
1025365900350
El Guendouzis, M. & Dinane, A. (2000). Determination of water activities, osmotic
and activity coefficients in aqueous solutions using the hygrometric method. The
Journal of Chemical Thermodynamics, 32 (3), 297-310. [Link]
jcht.1999.0574
El Moussaouiti, M., Boistelle, R., Bouhaouss, A. & Klein, J. (1997). Crystalli-
zation of calcium sulphate hemihydrate in concentrated phosphoric acid solutions.
Chemical Engineering Journal, 68 (2-3), 123-130. [Link]
8947(97)00116-2
Elguendouzi, M., Rifai, A. & Skafi, M. (2015). Properties of fluoride in wet phos-
phoric acid processes : Fluorosilicic acid in an aqueous solution of H2SiF6–H20
at temperatures ranging from 298.15 K to 353.15 K. Fluid Phase Equilibria, 396,
43-49. [Link]
Elmisaoui, S., Latifi, A. M., Khamar, L. & Salouhi, M. (2021). Analysis of the dis-
solution mechanism in the phosphoric acid manufacturing process : modelling and
simulation. Computer Aided Chemical Engineering (p. 891-897). Elsevier. https:
//[Link]/10.1016/B978-0-323-88506-5.50138-8
Elmore, K., Hatfield, J., Dunn, R. & Jones, A. (1965). Dissociation of phosphoric
acid solutions at 25. The Journal of Physical Chemistry, 69 (10), 3520-3525. https:
//[Link]/10.1021/j100894a045
Elyamani, Y., Skafi, M., Rifai, A. & Guendouzi, M. E. (2021). Solubility behaviors of
(Na+, K+ or NH4+) hexafluoridosilicates in H2SO4, HF, H2SiF6, H3PO4 acidic
aqueous solutions at T= 353.15 K. Journal of Fluorine Chemistry, 247, 109796.
[Link]
Epstein, N. (1997). Elements of particle deposition onto nonporous solid surfaces paral-
lel to suspension flows. Experimental Thermal and Fluid Science, 14 (4), 323-334.
[Link]
Feldmann, T. & Demopoulos, G. P. [George P.]. (2013). Influence of Impurities on
Crystallization Kinetics of Calcium Sulfate Dihydrate and Hemihydrate in Strong
HCl-CaCl2 Solutions. Industrial & Engineering Chemistry Research, 52, 6540-6549.
Förster, M., Augustin, W. & Bohnet, M. (1999). Influence of the adhesion force
crystal/heat exchanger surface on fouling mitigation. Chemical Engineering and
Processing : Process Intensification, 38 (4-6), 449-461. [Link]
S0255-2701(99)00042-2
Förster, M. [Markus] & Bohnet, M. (2000). Modification of molecular interactions
at the interface crystal/heat transfer surface to minimize heat exchanger fouling.
International Journal of Thermal Sciences, 39 (7), 697-708. [Link]
1016/S1290-0729(00)00229-5
Foti, C., Gianguzza, A. & Sammartano, S. (1997). A comparison of equations for
fitting protonation constants of carboxylic acids in aqueous tetramethylammonium
chloride at various ionic strengths. Journal of solution chemistry, 26 (6), 631-648.
[Link]
163
BIBLIOGRAPHIE
164
BIBLIOGRAPHIE
165
BIBLIOGRAPHIE
166
BIBLIOGRAPHIE
Langman, M. (1986). Towards estimation and confidence intervals. British medical jour-
nal (Clinical research ed.), 292 (6522), 716. [Link]
716
Lapidot, T., Matar, O. K. & Heng, J. Y. (2019). Calcium sulphate crystallisation in
the presence of mesoporous silica particles : Experiments and population balance
modelling. Chemical Engineering Science, 202, 238-249.
Lassis, M., Mizane, A., Dadda, N. & Rehamnia, R. (2015). Dissolution of Djebel Onk
phosphate ore using sulfuric acid. Environmental Nanotechnology, Monitoring &
Management, 4, 12-16. [Link]
Lei, M., Vallières, C., Grévillot, G. & Latifi, A. M. (2013). Thermal swing ad-
sorption process for carbon dioxide capture and recovery : modeling, simulation,
parameters estimability, and identification. Industrial & Engineering Chemistry
Research, 52 (22), 7526-7533. [Link]
Leikam, D. F. & Achorn, F. P. (2005). Phosphate fertilizers : Production, characteris-
tics, and technologies. Phosphorus : Agriculture and the environment, 46, 23-50.
[Link]
Letcher, T. M. (2007). Thermodynamics, solubility and environmental issues. Elsevier.
Li, D. [Dongdong], Zeng, D., Yin, X., Han, H., Guo, L. & Yao, Y. (2016). Phase
diagrams and thermochemical modeling of salt lake brine systems. II. NaCl+ H2O,
KCl+ H2O, MgCl2+ H2O and CaCl2+ H2O systems. Calphad, 53, 78-89. https:
//[Link]/10.1016/[Link].2016.03.007
Li, D. [Dongyue], Li, Z. & Gao, Z. (2019). Quadrature-based moment methods for the
population balance equation : An algorithm review. Chinese Journal of Chemical
Engineering, 27 (3), 483-500. [Link]
Li, J., Polka, H.-M. & Gmehling, J. (1994). A gE model for single and mixed solvent
electrolyte systems : 1. Model and results for strong electrolytes. Fluid Phase Equi-
libria, 94, 89-114.
Li, R., Wang, J. J., Zhou, B., Awasthi, M. K., Ali, A., Zhang, Z., Lahori, A. H. &
Mahar, A. (2016). Recovery of phosphate from aqueous solution by magnesium
oxide decorated magnetic biochar and its potential as phosphate-based fertilizer
substitute. Bioresource technology, 215, 209-214. https : / / doi . org / 10 . 1016 / j .
biortech.2016.02.125
Li, Z. & Demopoulos, G. P. [George P]. (2006). Development of an improved chemical
model for the estimation of CaSO4 solubilities in the HCl- CaCl2- H2O system up
to 100° C. Industrial & engineering chemistry research, 45 (9), 2914-2922. https:
//[Link]/10.1021/ie0508280
Lide, D. R. (2005). CRC handbook of chemistry and physics, internet version 2005.
Liu, A. & Rigopoulos, S. (2019). A conservative method for numerical solution of the
population balance equation, and application to soot formation. Combustion and
Flame, 205, 506-521. [Link]
Liu, J. & Benyahia, B. (2022). Optimal start-up strategies of a combined cooling and
antisolvent multistage continuous crystallization process. Computers & Chemical
Engineering, 107671. [Link]
Liu, J.-L. & Li, C.-L. (2019). A generalized Debye-Hückel theory of electrolyte solutions.
AIP Advances, 9 (1), 015214. [Link]
167
BIBLIOGRAPHIE
Long, B., Chen, S., Li, Z. & Ding, Y. (2020). Solubility of ammonium fluoride in
aqueous sodium fluoride solutions of various concentrations from dilute to satura-
ted at 298.15 K for fluorine recovery from wet-process phosphoric acid. Journal of
Molecular Liquids, 299, 112247. [Link]
Ma, H., Feng, X. & Deng, C. (2018). Water–Phosphorus Nexus for Wet-Process Phos-
phoric Acid Production. Industrial & Engineering Chemistry Research, 57 (20),
6968-6979. [Link]
Mantzaris, N. V., Daoutidis, P. & Srienc, F. (2001). Numerical solution of multi-
variable cell population balance models. III. Finite element methods. Computers
& Chemical Engineering, 25 (11-12), 1463-1481.
Marion, G. M. & Farren, R. E. (1999). Mineral solubilities in the Na-K-Mg-Ca-Cl-SO4-
H2O system : A re-evaluation of the sulfate chemistry in the Spencer-Møller-Weare
model. Geochimica et Cosmochimica Acta, 63 (9), 1305-1318. [Link]
1016/S0016-7037(99)00102-7
Marliacy, P., Solimando, R., Bouroukba, M. & Schuffenecker, L. (2000). Ther-
modynamics of crystallization of sodium sulfate decahydrate in H2O–NaCl–Na2SO4 :
application to Na2SO4· 10H2O-based latent heat storage materials. Thermochi-
mica Acta, 344 (1-2), 85-94. [Link]
Marshall, W. L. & Slusher, R. (1966). Thermodynamics of Calcium Sulfate Dihy-
drate in Aqueous Sodium Chloride Solutions, 0-110° 1, 2. The Journal of Physical
Chemistry, 70 (12), 4015-4027. [Link]
Mateo, J. R. S. C. (2012). Multi-attribute utility theory. Multi Criteria Analysis in the
Renewable Energy Industry (p. 63-72). Springer.
Mazloumi, S. H. & Shojaeian, A. (2019). Modified nonelectrolyte Wilson-NRF : A
new model for strong and weak electrolyte solutions. Journal of Molecular Liquids,
277, 714-725. [Link]
McCleskey, R. B., Nordstrom, D. K., Ryan, J. N. & Ball, J. W. (2012). A new
method of calculating electrical conductivity with applications to natural waters.
Geochimica et Cosmochimica Acta, 77, 369-382. [Link]
2011.10.031
McNaught, A. D. & Wilkinson, A. (1997). IUPAC compendium of chemical termino-
logy Oxford. UK : Blackwell Science.
Messnaoui, B. (2008). Representation of excess properties and liquid composition of
aqueous solutions of the HF+ water system. Journal of solution chemistry, 37 (6),
715-726. [Link]
Messnaoui, B., Bouhaouss, A. & Derja, A. (2012). A steady state model for describing
the phosphate lattice loss. Procedia Engineering, 46, 134-142. [Link]
1016/[Link].2012.09.456
Messnaoui, B. & Bounahmidi, T. (2005). Modeling of excess properties and vapor–
liquid equilibrium of the system H3PO4–H2O. Fluid Phase Equilibria, 237 (1-2),
77-85. [Link]
Messnaoui, B. & Bounahmidi, T. (2006). On the modeling of calcium sulfate solubility
in aqueous solutions. Fluid phase equilibria, 244 (2), 117-127. [Link]
1016/[Link].2006.03.022
168
BIBLIOGRAPHIE
Mohs, A. & Gmehling, J. (2013). A revised LIQUAC and LIFAC model (LIQUAC*/LIFAC*)
for the prediction of properties of electrolyte containing solutions. Fluid Phase
Equilibria, 337, 311-322. [Link]
Møller, N. (1988). The prediction of mineral solubilities in natural waters : A chemical
equilibrium model for the Na-Ca-Cl-SO4-H2O system, to high temperature and
concentration. Geochimica et Cosmochimica Acta, 52 (4), 821-837. https : / / doi .
org/10.1016/0016-7037(80)90287-2
Mullin, J., Chakraborty, M. & Mehta, K. (1970). Nucleation and growth of am-
monium sulphate crystals from aqueous solution. Journal of Applied Chemistry,
20 (12), 367-371. [Link]
Ngo, V. V., Michel, J., Gujisaite, V., Latifi, A. & Simonnot, M.-O. (2014). Para-
meters describing nonequilibrium transport of polycyclic aromatic hydrocarbons
through contaminated soil columns : Estimability analysis, correlation, and opti-
mization. Journal of contaminant hydrology, 158, 93-109. [Link]
[Link].2014.01.005
Onyemelukwe, I., Benyahia, B., Reis, N. M., Nagy, Z. K. & Rielly, C. D. (2018).
The heat transfer characteristics of a mesoscale continuous oscillatory flow crystal-
liser with smooth periodic constrictions. International Journal of Heat and Mass
Transfer, 123, 1109-1119. [Link]
015
Pabalan, R. T. & Pitzer, K. S. (1987). Thermodynamics of concentrated electrolyte
mixtures and the prediction of mineral solubilities to high temperatures for mix-
tures in the system Na-K-Mg-Cl-SO4-OH-H2O. Geochimica et Cosmochimica Acta,
51 (9), 2429-2443. [Link]
Papadopoulos, A. I. & Seferlis, P. (2009). Generic modelling, design and optimiza-
tion of industrial phosphoric acid production processes. Chemical Engineering and
Processing : Process Intensification, 48 (1), 493-506. [Link]
cep.2008.06.011
Park, G.-J. (2007). Analytic methods for design practice. Springer Science & Business
Media.
Parkhurst, D. L. & Appelo, C. (1999). User’s guide to PHREEQC (Version 2) : A
computer program for speciation, batch-reaction, one-dimensional transport, and
inverse geochemical calculations. Water-resources investigations report, 99 (4259),
312.
Peleg, M., Normand, M. D. & Corradini, M. G. (2012). The Arrhenius equation
revisited. Critical Reviews in Food Science and Nutrition, 52 (9), 830-851. https:
//[Link]/10.1080/10408398.2012.667460
Peng, D.-Y. & Robinson, D. B. (1976). A new two-constant equation of state. Industrial
& Engineering Chemistry Fundamentals, 15 (1), 59-64. [Link]
i160057a011
Peng, Y., Zhu, Z., Braatz, R. D. & Myerson, A. S. (2015). Gypsum crystallization
during phosphoric acid production : modeling and experiments using the mixed-
solvent-electrolyte thermodynamic model. Industrial & Engineering Chemistry Re-
search, 54 (32), 7914-7924. [Link]
169
BIBLIOGRAPHIE
170
BIBLIOGRAPHIE
Renner, T. (2007). Quantities, units and symbols in physical chemistry. Royal Society
of Chemistry.
Renom Carrasco, M., Nikitine, C., Hamou, M., de Bellefon, C., Thieuleux, C. &
Meille, V. (2020). Self-Metathesis of Methyl Oleate Using Ru-NHC Complexes :
A Kinetic Study. Catalysts, 10 (4), 435.
Richter, U., Brand, P., Bohmhammel, K. & Könnecke, T. (2000). Thermodynamic
investigations of aqueous solutions of aluminium chloride. The Journal of Chemical
Thermodynamics, 32 (2), 145-154. [Link]
Rifai, A., El Guendouzi, M., Amalhay, M. & Azaroual, M. (2009). Thermodynamic
properties of aqueous fluorosilicate acid H2SiF6 (aq) at 298.15 Kand 333.15 K.
Proceedings COVAPHOSIII, 5, 125-129.
Rumyantsev, A. V., Hagemann, S. & Moog, H. C. (2004). Isopiestic Investigation of
the Systems Fe2 (SO4) 3–H2SO4–H2O, FeCl3–H2O, and Fe (III)–(Na, K, Mg, Ca)
Cln–H2O at 298.15 K. Zeitschrift für Physikalische Chemie, 218 (9), 1089-1127.
[Link]
Sahinidis, N. V. (2002). Global optimization and constraint satisfaction : The branch-
and-reduce approach. International Workshop on Global Optimization and Constraint
Satisfaction, 1-16. [Link]
Saltelli, A. (2002). Making best use of model evaluations to compute sensitivity indices.
Computer Physics Communications, 145 (2), 280-297. [Link]
S0010-4655(02)00280-1
Salvatori, F., Muhr, H., Plasari, E. & Bossoutrot, J. (2002). Determination of
the nucleation and crystal growth kinetics of precipitation processes involving si-
multaneous gas release. Chemical Engineering Transactions, 1, 293-298. https :
//[Link]/10.1016/S0009-2509(98)00488-6
Sander, B., Fredenslund, A. & Rasmussen, P. (1986). Calculation of vapour-liquid
equilibria in mixed solvent/salt systems using an extended UNIQUAC equation.
Chemical Engineering Science, 41 (5), 1171-1183. [Link]
2509(86)87090-7
Shen, Y., Linnow, K. & Steiger, M. (2020). Crystallization behavior and damage
potential of Na2SO4–NaCl mixtures in porous building materials. Crystal Growth
& Design, 20 (9), 5974-5985. [Link]
Shenn, L., Sippola, H., Li, X., Lindberg, D. & Taskinen, P. (2019). Thermodynamic
Modeling of Calcium Sulfate Hydrates in the CaSO4–H2O System from 273.15 to
473.15 K with Extension to 548.15 K. Journal of Chemical & Engineering Data,
64 (6), 2697-2709. [Link]
Shenn, L., Sippola, H., Li, X., Lindberg, D. & Taskinen, P. (2020). Thermodyna-
mic Modeling of Calcium Sulfate Hydrates in a CaSO4–H2SO4–H2O System from
273.15 to 473.15 K up to 5 m Sulfuric Acid. Journal of Chemical & Engineering
Data, 65 (5), 2310-2324. [Link]
Shiah, I.-M. & Tseng, H.-C. (1996). Experimental and theoretical determination of
vapor pressures of NaCl KCl, NaBr KBr and NaCl CaCl2 aqueous solutions at 298
to 343 K. Fluid phase equilibria, 124 (1-2), 235-249. [Link]
3812(96)03013-0
171
BIBLIOGRAPHIE
Simoes, M. C., Hughes, K. J., Ingham, D. B., Ma, L. & Pourkashanian, M. (2017).
Temperature Dependence of the Parameters in the Pitzer Equations. Journal of
Chemical & Engineering Data, 62 (7), 2000-2013. [Link]
jced.7b00022
Sippola, H. & Taskinen, P. (2014). Thermodynamic properties of aqueous sulfuric acid.
Journal of Chemical & Engineering Data, 59 (8), 2389-2407. [Link]
1021/je4011147
Skafi, M., Elyamani, Y. & El Guendouzi, M. (2019). Solubility of Mixed Ternary
Systems of Hexafluoridosilicate Salts (M= Na+, K+, or NH4+) in Hexafluoridosi-
licic Acid Aqueous Solutions at T= 353.15 K. Journal of Chemical & Engineering
Data, 64 (8), 3290-3299. [Link]
Slack, A. (1977). Commercial ammonia processes. Fertilizer Science and Technology
Series.
Slack, A. & James, G. (1973). Ammonia. Fertilizer Science and Technology Series.
Marcel Deckker. Inc. New York.
Smith, R. C. (2013). Uncertainty quantification : theory, implementation, and applica-
tions (T. 12). Siam.
Sobol, I. M. (1990). On sensitivity estimation for nonlinear mathematical models. Ma-
tematicheskoe modelirovanie, 2 (1), 112-118.
Spencer, R. J., Møller, N. & Weare, J. H. (1990). The prediction of mineral solu-
bilities in natural waters : A chemical equilibrium model for the Na- K-Ca- Mg-
Cl- SO4-H2O system at temperatures below 25° C. Geochimica et Cosmochimica
Acta, 54 (3), 575-590. [Link]
Steiger, M. & Asmussen, S. (2008). Crystallization of sodium sulfate phases in porous
materials : the phase diagram Na2SO4–H2O and the generation of stress. Geochi-
mica et Cosmochimica Acta, 72 (17), 4291-4306. [Link]
2008.05.053
Steiger, M., Linnow, K., Ehrhardt, D. & Rohde, M. (2011). Decomposition reac-
tions of magnesium sulfate hydrates and phase equilibria in the MgSO4–H2O and
Na+–Mg2+–Cl—SO42—H2O systems with implications for Mars. Geochimica et
Cosmochimica Acta, 75 (12), 3600-3626. [Link]
Tian, P., Ning, P., Cao, H. & Li, Z. (2012). Determination and Modeling of Solubility
for CaSO4· 2H2O–NH4+–Cl—-SO42—-NO3—-H2O System. Journal of Chemical
& Engineering Data, 57 (12), 3664-3671. [Link]
Tsang, T. & Rao, A. (1990). A moving finite element method for the population balance
equation. International journal for numerical methods in fluids, 10 (7), 753-769.
Tuncer, E., Serdyuk, Y. V. & Gubanski, S. M. (2002). Dielectric mixtures : elec-
trical properties and modeling. IEEE Transactions on Dielectrics and Electrical
Insulation, 9 (5), 809-828. [Link]
Waerstad, K. R. (1985). Quantitative X-ray Diffraction Analysis of Calcium Sulfates and
Quartz in Wet-Process Phosphoric Acid Filter Cakes. Advances in X-ray Analysis,
29, 211-216. [Link]
Wang, J., Petit, C., Zhang, X. & Cui, S. (2016). Phase equilibrium study of the
AlCl3–CaCl2–H2O system for the production of aluminum chloride hexahydrate
172
BIBLIOGRAPHIE
from Ca-rich flue ash. Journal of Chemical & Engineering Data, 61 (1), 359-369.
[Link]
Wang, M., Kaur, H. & Chen, C.-C. (2017). Thermodynamic modeling of HNO3-H2SO4-
H2O ternary system with symmetric electrolyte NRTL model. AIChE Journal,
63 (7), 3110-3117. [Link]
Wang, P., Anderko, A. & Young, R. D. (2002). A speciation-based model for mixed-
solvent electrolyte systems. Fluid Phase Equilibria, 203 (1-2), 141-176. [Link]
org/10.1016/S0378-3812(02)00178-4
Wang, P., Pitzer, K. S. & Simonson, J. M. (1998). Thermodynamic properties of
aqueous magnesium chloride solutions from 250 to 600 K and to 100 MPa. Journal
of Physical and Chemical Reference Data, 27 (5), 971-991. https : / / doi . org / 10 .
1063/1.556026
Wang, W., Zeng, D., Chen, Q. & Yin, X. (2013). Experimental determination and mo-
deling of gypsum and insoluble anhydrite solubility in the system CaSO4–H2SO4–
H2O. Chemical Engineering Science, 101, 120-129. [Link]
2013.06.023
Wangi, W., Zeng, D., Wang, J., Li, H. & Wu, L. (2014). Experimental determination
and modeling of the solubility of CaSO4· 2H2O and CaSO4 in the quaternary
system CaSO4+ MgSO4+ H2SO4+ H2O. Industrial & Engineering Chemistry
Research, 53 (32), 12839-12847. [Link]
Watkinson, A. & Wilson, D. (1997). Chemical reaction fouling : A review. Experimental
Thermal and Fluid Science, 14 (4), 361-374. https : / / doi . org / 10 . 1016 / S0894 -
1777(96)00138-0
White, H. J., Parker, V. B. & Garvin, D. (1987). Codata Thermodynamic Tables :
Selections for Some Compounds of Calcium and Related Mixtures : A Prototype
Set of Tables.
Wilson, G. M. (1964). Vapor-liquid equilibrium. XI. A new expression for the excess
free energy of mixing. Journal of the American Chemical Society, 86 (2), 127-130.
[Link]
Wolf, A. V. (1966). Aqueous solutions and body fluids ; their concentrative properties
and conversion tables.
Wood, J. R. (1975). Thermodynamics of brine-salt equilibria—I. The systems NaCl-KCl-
MgCl2-CaCl2-H2O and NaCl-MgSO4-H2O at 25 C. Geochimica et Cosmochimica
Acta, 39 (8), 1147-1163. [Link]
Wright, M. R. (2007). An introduction to aqueous electrolyte solutions. John Wiley &
Sons. [Link]
Wu, S., McLean, K. A., Harris, T. J. & McAuley, K. B. (2011). Selection of optimal
parameter set using estimability analysis and MSE-based model-selection criterion.
International Journal of Advanced Mechatronic Systems, 3 (3), 188-197. https://
[Link]/10.1504/IJAMECHS.2011.042615
Yang, H., Zhao, Z., Zeng, D. & Yin, R. (2016). Isopiestic Measurements of Water
Activity for the H2SO4–H3PO4–H2O System at 298.15 K. Journal of Solution
Chemistry, 45 (11), 1580-1587. [Link]
Yangx, L., Cao, J. & Luo, T. (2018). Effect of Mg2+, Al3+, and Fe3+ ions on crys-
tallization of type α hemi-hydrated calcium sulfate under simulated conditions of
173
BIBLIOGRAPHIE
174
Résumé
Ce travail de thèse s’intéresse au développement d’un modèle thermodynamique et à son exploitation
dans la modélisation et l’optimisation des procédés de fabrication d’acide phosphorique. Il est décrit
par des équations de bilans de matière et de charge, d’équilibres chimiques et du modèle de Pitzer,
et met en jeu plusieurs paramètres inconnus à identifier à partir de mesures expérimentales. Dans cet
objectif, une base de données expérimentales est utilisée. Elle contient des mesures de spéciation d’acides
phosphorique et sulfurique, des mesures de solubilité de dix minéraux, et des mesures de l’activité de
l’eau de huit systèmes binaires. Ces mesures sont effectuées dans des conditions de température allant
de 298K à 353K et des concentrations allant de 0 à 20 mol/kg d’eau. Une méthode fondée sur l’analyse
des sensibilités globales a ensuite été développée et utilisée pour évaluer l’estimabilité des paramètres
inconnus. Les paramètres les plus estimables sont identifiés, et les valeurs de ceux qui ne le sont pas sont
fixées à partir de précédents travaux ou de la littérature. Le modèle identifié a ensuite été validé à l’aide
de mesures qui n’ont pas servies pour son identification. La validation est basée sur le calcul des indices
de performance, et sur l’utilisation du test statistique F de Fisher-Snedecor. Le modèle ainsi identifié
et validé est exploité dans l’optimisation muliticritère d’une unité industrielle de production d’acide
phosphorique, notamment pour la minimisation des pertes chimiques en phosphate et l’amélioration du
rendement. Il est ensuite utilisé dans l’analyse des problèmes d’encrassement qui altèrent les performances
des unités de production. Enfin, le modèle est mis à profit pour la modélisation et la simulation de
la cristallisation des sulfates de calcium lors de la production de l’acide phosphorique. L’ensemble des
résultats obtenus montrent que le modèle développé et validé peut être utilisé comme outil prédictif précis
pour la conception et le fonctionnement optimal des procédés de fabrication d’acide phosphorique actuels,
voire pour le développement de nouveaux procédés plus compacts, plus intégrés et plus performants.
Abstract
This PhD work deals with the development of a thermodynamic model and its use in the modeling
and optimization of phosphoric acid manufacturing processes. The model is based on mass and charge
balance equations, along with chemical equilibria equations and Pitzer model equations. It involves several
unknown parameters that should be identified from experiments. To this end, an experimental database
is used. It contains speciation measurements of phosphoric and sulfuric acids, solubility measurements
of ten minerals, and water activity of eight binary systems. These measurements are performed under
temperature conditions ranging from 298K to 353K and concentrations ranging from 0 to 20 mol/kg water.
A global sensitivity based estimability analysis is then developed and used to evaluate the estimability
of the unknown parameters from the available data. The estimable parameters are then identified, and
the values of the non-estimable ones are taken from previous works or from the literature. The identified
model is then validated by means of additional measurements that did not serve for its identification. The
validation is based on the computation of the performance indices, and on the use of the Fisher-Snedecor
F statistical test. The identified and validated model is exploited in a multi-objective optimization of
an industrial phosphoric acid production unit, especially for minimizing chemical phosphate losses and
performance improvement. It is then used in the study and analysis of fouling problems that degrade
the performance of production units. Finally, the model is exploited in the modeling and simulation of
the crystallization of calcium sulfates during the production of phosphoric acid. All the results obtained
in this work show that the developed and validated model can be used as an accurate predictive tool in
the design and optimal operation of current phosphoric acid manufacturing processes, and even in the
development of new processes which are more compact, integrated and efficient.
175
Keywords: Modelling, Optimisation, Thermodynamics, Phosphoric acid, Experimental measurements,
Global estimability analysis.
176