Question

方程的求解

Answer 1

6.1.3.1 数值解法

方程（6-1-25）余值是a_k（k=1，2，...，m）的函数，可以写成如下形式：

水文地球化学

式中：a=［a₁，a₂，…，a_m］，a的值越接近于方程组（6-1-25）的解，|g_k（a）|就越小。这样，求解方程组（6-1-25）的问题即转化为求出一组a_k值（a_k≥0，这是由a_k的物理意义所决定的），使|g_k（a）|尽可能地小。达到这种目的途径有多种，由于Newton-Raphon迭代法具有二次收敛性，使求解有较快的收敛速度，因此我们采用这种方法

把式（6-1-26）在一组初值处展开为泰勒级数有：

水文地球化学

在Δa_k（k=1，2，…，m）趋近于方程组（6-1-25）的解的情况下，|+Δa₁，+Δa，…，+Δa_m）|即趋近于0，忽略高阶无穷小项0（Δa²），便可由式（6-1-27）得到：

水文地球化学

方程（6-1-28）为一含有m个未知数（Δa₁，Δa₂，...，Δa_m）的线性方程组。在给定一组初值（k=1，2，…，m）的情况下，可由方程组（6-1-28）求得基本组分浓度a_k的绝对校正量Δa_k。

由式（6-1-26）可得：

水文地球化学

由式（6-1-28）、（6-1-29）和（6-1-30）可以看出，与其求a_k的绝对校正量Δa_k，还不如求a_k的相对校正量Δa_k/a_k来得方便。若令：

{X}表示下述的未知数列向量：

水文地球化学

［Q］为具有如下元素的m×m阶对称矩阵：

水文地球化学

{B}为常数项：

水文地球化学

则线性方程组（6-1-28）可写成下面的矩阵形式：

水文地球化学

方程组（6-1-35）可采用高斯全主元消去法求解

假定根据第r次迭代值已求出了与其相应的相对校正量则第 r+1 次迭代值可由下式求出：

水文地球化学

以第r+1次迭代值为基础，可求得a_k的第r+1次相对校正量，代入式（6-1-36）又可求得第r+2次迭代值。这样一步一步地校正下去，当满足某一预先给定的误差限后，即可用所求得的a_k值作为主要组分的浓度。将此结果代入式（6 1 23），便可求得所有衍生组分的浓度。在后面的计算实例中，当δ=max（x₁，x₂，…，x_m）＜10 ^-6时，即认为迭代已达到了所要求的精度。

6.1.3.2 活度系数的校正

在前面的讨论中，我们假定所有水溶组分的活度系数均为1，但在天然水溶液中，由于离子间的相互作用，必须考虑活度系数的影响。在考虑活度系数的情况下，式（6-1-22）应具有下述的形式：

水文地球化学

式中大括弧标记组分的活度，活度与浓度仅差一比例常数——活度系数，即：

水文地球化学

这里α_k、β_j分别表示第k种基本组分和第j种衍生组分的活度系数，可根据具体情况使用第二章所介绍的不同形式的Debye-Hükel公式来计算，计算中水溶液的离子强度由下式定义：

水文地球化学

由式（6-1-38）～（6-1-39）可得：

水文地球化学

这里，F_j为活度系数校正量，可用下式计算：

水文地球化学

当对活度系数进行校正时，式（6-1-26）变为：

水文地球化学

活度系数的校正虽然没有改变求解线性方程组（6-1-35）的形式，但其内容已经发生了改变。例如，给定一组a_k的初值后，向量{B}中的各元素需用（6-1-43）式而不是（6-1-26）式来确定；方程的系数矩阵［Q］中的各元素则由下面的关系式来计算：

水文地球化学

确定线性方程组（6-1-35）的系数矩阵［Q］和常数项{B}的另一个问题就是式（6-1-43）～（6-1-45）中化学反应平衡常数K_j的计算，在已知组分的标准生成自由能资料或化学反应的标准焓变和标准熵变的情况下，使用第二章所介绍的方法可方便地进行反应平衡常数的计算。同时，前人在化学反应的平衡常数计算方面也积累了大量的经验值和经验公式，在对其进行对比分析的情况下也可直接引用或根据具体条件进行计算。

6.1.3.3 迭代步骤

由上面的讨论可知，水溶组分活度系数的计算需要已知a_k、c_j的值，而a_k与c_j值的计算又需要已知组分的活度系数，这种似乎无法解决的问题可用迭代法求解，具体的迭代步骤如下：

（1）计算水溶液中由基本组分形成衍生组分的各化学反应的平衡常数K_j；

（2）用各基本组分的化学分析浓度作为a_k（k=1，2，…，m）的初次估计值，取衍生组分的浓度c_j（j=1，2，…，n）等于0，由式（6-1-40）计算水溶液的离子强度；

（3）以上述I的计算值为基础，使用Debye Hükel公式计算各组分的活度系数，并用式（6-1-42）计算活度系数校正量F_j；

（4）使用式（6-1-43）计算方程组（6-1-35）的常数项，使用式（6-1-44）和（6-1-45）计算方程组（6-1-35）的系数矩阵，在此基础上使用高斯全主元消去法求解线性方程组（6-1-35），进而求得a_k（k=1，2，…，m），将a_k代入式（6-1-41）计算c_j（j=1，2，…，n）；

（5）根据第四步所求得的a_k、c_j值，用式（6-1-40）计算水溶液的离子强度，并分别用Debye-Hükel公式和式（6-1-42）计算各组分的活度系数及活度系数校正量F_j；

（6）重复第四、第五步，直到相邻两次离子强度的计算结果Iⁿ、Iⁿ⁺¹满足条件：

水文地球化学

这时所求得的a_k、c_j值即为水溶组分的平衡分布。

在实际计算中，由于已测得了水溶液的pH值，因此基本组分H⁺和OH^-的浓度通常是由下述两式计算的：

水文地球化学

这样实际所要求解的方程个数就不是m个，而是（m-2）个。