Variants of application of the least squares method in Szyszkowski and Rosin–Rammler approximations

全文:

Введение. Особенности экспериментов

В настоящее время повышенное внимание уделяется поверхностным явлениям, так как именно эти явления вызывают повышенный интерес научных и инженерно-технических работников. Возникает необходимость важнейшие фундаментальные закономерности представлять в виде количественных соотношений, что позволяет глубже раскрывать физический смысл явлений и определять границы применимости получаемых соотношений. Закономерности протекания поверхностных явлений служат теоретической основой получения материалов с заданными свойствами. Основой изучения закономерностей является эксперимент и его математическая обработка.

Опишем особенности проведения экспериментов, из которых получаются данные для зависимостей Шишковского и Розина–Раммлера [1–3]. Для выявления этих зависимостей и определения их параметров используется метод наименьших квадратов (МНК) [4]. Эмпирическая зависимость Шишковского применяется для определения и изучения поверхностного натяжения раствора в зависимости от концентрации в нём поверхностно активного вещества (ПАВ) при средних концентрациях ПАВ и исследования закономерностей адсорбции ПАВ из водного раствора на твёрдых адсорбентах. В свою очередь зависимость Розина–Раммлера применяется для описания интегральной кривой накопления в координатах «время – относительные весовые доли» и позволяет получить дифференциальное распределение всех весовых долей по размерам частиц. Это распределение хорошо описывает результат воздействия на материал процесса осаждения взвешенных частиц в седиментационном анализе или результат механического разрушения материала [5].

Что касается получения зависимости Шишковского, то здесь при проведении эксперимента часто используется так называемый капельный метод. Суть его заключалась в измерении на сталагмометре СТ-2 объёма капли исследуемого раствора ПАВ и определении величины поверхностного натяжения с дальнейшим расчётом и построением изотерм адсорбции. Для повышения точности эксперимента проводится 100 итераций для каждого раствора ПАВ с последующим расчётом среднего значения объёма капли. По полученным точкам строится аппроксимирующая зависимость Шишковского, которая даёт связь коэффициента поверхностного натяжения (КПН) с концентрацией ПАВ [6–9].

Использование зависимости Розина–Раммлера возможно для данных, подготовленных методом непрерывного взвешивания седиментационного осадка. Целью является получение кривых распределения, анализ которых позволяет установить фракционный состав системы. В этом методе площадка весов находится в жидкости на заданной глубине. При этом взаимного влияния частиц и влияния стенок нет. Частицы сферические, одинаковой плотности. Их постоянная скорость определяется из закона Стокса и зависит от диаметра частицы, её плотности, вязкости и плотности жидкости [10, 11]. В процессе эксперимента идёт регистрация интегрального веса осадка через определённые промежутки времени. Полученная кривая осаждения аппроксимируется зависимостью Розина–Раммлера, по которой определяются интегральная и дифференциальная кривые по размерам частиц.

Особенность обработки данных экспериментов изложена отдельно для зависимостей Розина–Раммлера и Шишковского.

О применении метода наименьших квадратов

Пусть имеются экспериментальные (с индексом «е») данные

$\left\{\left(w_{i,e};z_{i,e}\right)\right\};i=\mathrm{1,}\dots ,N,$          (1)

где $w_{i,e}$– вес или КПН; $z_{i,e}$– время или концентрация ПАВ.

Предполагается, что данные должны аппроксимировать некоторую двухпараметрическую зависимость

$w=w\left(z,m_1,m_2\right)$                    (2)

с неизвестными параметрами $m_1$ и $m_2$, которые и требуется определить. Даже если мы знаем значения истинной зависимости (2), равенства

$w_{i,e}=w\left(z_{i,e},m_1,m_2\right);i=\mathrm{1,}\dots ,N,$ (3)

могут выполняться лишь приближённо в силу погрешностей эксперимента. Поэтому для экспериментальных данных (1) может оказаться, что такие параметры $m_1$ и $m_2$, обеспечивающие строгое выполнение всех равенств (3), подобрать невозможно. Тогда мы будем добиваться не строгого выполнения равенств (3), а наименьшей погрешности в этих равенствах. Математически это можно записать так, что нам нужно найти такие значения параметров $m_1$ и $m_2$, чтобы суммарная погрешность в равенствах (3) была наименьшей. Для этого вычисляются все отклонения (невязки) в равенствах (3) и суммируются квадраты этих отклонений, после чего находится минимум целевой функции

$S_0\left(m_1,m_2\right)={\sum }_{i=1}^N{\left[w_{i,e}-w\left(z_{i,e},m_1,m_2\right)\right]}^2.$ (4)

Ясно, что если для данных (1) есть параметры $m_1$ и $m_2$, обеспечивающие точное выполнение равенств (3), то значение целевой функции (4) при этих параметрах даст 0. Подбор таких параметров $m_1$ и $m_2$ в аппроксимирующей зависимости (2), которые доставляют минимум сумме (4), и носит название «метод наименьших квадратов» (МНК) [4], и запись задачи имеет общий вид:

$\left(m_1,m_2\right)=\arg \left\{\min \left(S_0\left(m_1,m_2\right)\right)\right\}.$     (5)

Здесь и далее считается, что аналитические выражения функций, образующих двухпараметрическую зависимость (2), позволяют дифференцировать эти функции по параметрам $m_1$ и $m_2$.

Непосредственное решение задачи (5) основывается на нелинейном программировании [12]. Кроме этого, используется необходимое условие минимума (НУМ):

$\frac{\partial S_0}{\partial m_1}=\mathrm{0,}\frac{\partial S_0}{\partial m_2}=0.$        (6)

Если искомые параметры входят в (2) линейно, или если после преобразований удаётся прийти к линейным выражениям, то из (6) получаются более удобные в применении алгебраические соотношения. Иногда в качестве таких преобразований удобно логарифмировать или экспонировать равенство (2), что будет вести к изменению целевой функции, но также обеспечивающей наименьшее отклонение в преобразованных равенствах (3).

Рассмотрим использование таких преобразований, когда получается аналитическое выражение для искомых параметров, но при этом возникает особенность типа $ln\left(0\right)$ или получается трансцендентное уравнение, которое их связывает.

Варианты реализации при аппроксимации соотношением Шишковского

Эмпирическое уравнение Шишковского описывает влияние концентраций ПАВ в растворе на КПН раствора

$\sigma ={\sigma }_0-B\cdot \ln \left(1+AC\right),$         (7)

где $\sigma$ – КПН раствора с ПАВ; ${\sigma }_0$– известный КПН чистого раствора; $A>0$, $B>0$– искомые коэффициенты аппроксимации, условие $B>0$ будет доказано, условие $A>0$ очевидно; C – концентрация ПАВ в растворе.

Заметим, что для равенства (7) можно использовать следующие эквивалентные записи:

$\ln \left({\sigma }_0-\sigma \right)=\ln \left(B\cdot \ln \left(1+AC\right)\right),$       (8)

$\exp \left(\frac{{\sigma }_0-\sigma }{B}\right)=1+AC,$          (9)

$\frac{{\sigma }_0-\sigma }{B}=\ln \left(1+AC\right).$            (10)

Пусть имеются данные из эксперимента: $i=\mathrm{1,}\dots , N; C_i>\mathrm{0,}$ здесь ${\sigma }_i$ – КПН, $C_i$ – концентрация ПАВ в растворе. Тогда сумма квадратов отклонений (невязок) для равенства (7) будет:

$S(A,B)={\sum }_{i=1}^N{\left[{\sigma }_0-{\sigma }_i-B\cdot \ln \left(1+A{C}_i\right)\right]}^2.$ (11)

Эквивалентные равенства (8)–(10) приводят к следующим целевым функциям:

$S_1\left(z,A\right)=$

${\sum }_{i=1}^N{\left[\ln \left({\sigma }_0-{\sigma }_i\right)-z-\ln \left(\ln \left(1+A{C}_i\right)\right)\right]}^2,z=\ln B,$ (12)

$S_2(A,B)={\sum }_{i=1}^N{\left[\exp \left(\frac{{\sigma }_0-{\sigma }_i}{B}\right)-1-A{C}_i\right]}^2,$ (13)

$S_3(A,B)={\sum }_{i=1}^N{\left[\frac{{\sigma }_0-{\sigma }_i}{B}-\ln \left(1+A{C}_i\right)\right]}^2.$ (14)

Рассмотрим варианты нахождения неизвестных параметров A и B.

Вариант 1. Методы нелинейного программирования применяются непосредственно к (11) или к аналогам (12)–(14), ищется минимум и определяются значения A и B, при которых этот минимум достигается:

$\left(A,B\right)=\arg \left\{\min \left[S\left(A,B\right)\right]\right\}.$ (15)

Вариант 2. Для целевой функции (11) из НУМ $\frac{\partial S_1}{\partial z}=0$ получаем:

$B=\frac{{\displaystyle {\sum }_{i=1}^N\left[\left({\sigma }_0-{\sigma }_i\right)\cdot \ln \left(1+A{C}_i\right)\right]}}{{\displaystyle {\sum }_{i=1}^N{\left[\ln \left(1+A{C}_i\right)\right]}^2}}.$    (16)

Другое условие минимума $\frac{\partial S}{\partial A}=0$ приводит к равенству: ${\sum }_{i=1}^N\frac{\left({\sigma }_0-{\sigma }_i\right)C_i}{1+A{C}_i}-B{\sum }_{i=1}^N\frac{C_i\cdot \ln \left(1+A{C}_i\right)}{1+A{C}_i}=\mathrm{0,}$

подставив в которое выражение B из (16), приходим к трансцендентному уравнению относительно неизвестного параметра A:

${\sum }_{i=1}^N\frac{\left({\sigma }_0-{\sigma }_i\right)C_i}{1+A{C}_i}-\frac{{\displaystyle {\sum }_{i=1}^N\left[\left({\sigma }_0-{\sigma }_i\right)\cdot \ln \left(1+A{C}_i\right)\right]}}{{\displaystyle {\sum }_{i=1}^N{\left[\ln \left(1+A{C}_i\right)\right]}^2}}\times$

$\times {\sum }_{i=1}^N\frac{C_i\cdot \ln \left(1+A{C}_i\right)}{1+A{C}_i}=\mathrm{0,}$

решив которое затем определяем параметр B из (16).

Вариант 3. НУМ $\frac{\partial S_1}{\partial z}=0$ для целевой функции (12) приводит к:

$z=\frac{1}{N}\left[{\sum }_{i=1}^N\ln \left({\sigma }_0-{\sigma }_i\right)-{\sum }_{i=1}^N\mathrm{lnln}\left(1+A{C}_i\right)\right],$ (17)

Или

$B=\exp \left\{\frac{1}{N}\left[\sum _{i=1}^N\ln \left({\sigma }_0-{\sigma }_i\right)-\right.\right.$

$-{\sum }_{i=1}^N\left.\left.\mathrm{lnln}\left(1+A{C}_i\right)\right]\right\},$ (18)

а другое НУМ $\frac{\partial S_1}{\partial A}=0$ для этой целевой функции позволяет записать равенство:

${\sum }_{i=1}^N\frac{\left[\ln \left({\sigma }_0-{\sigma }_i\right)-z-\text{ln}\left(\ln \left(1+A{C}_i\right)\right)\right]·C_i}{\left(1+A{C}_i\right)\text{ln}(1+A{C}_i)}=0.$ (19)

Подстановка z из (17) в (19) даёт трансцендентное уравнение относительно A:

${\sum }_{i=1}^N\ln \frac{{\sigma }_0-{\sigma }_i}{\ln \left(1+A{C}_i\right)}\frac{C_i}{\left(1+A{C}_i\right)\ln \left(1+A{C}_i\right)}$=

=$\frac{1}{N}{\sum }_{i=1}^N\ln \frac{{\sigma }_0-{\sigma }_i}{\ln \left(1+A{C}_i\right)}{\sum }_{i=1}^N\frac{C_i}{\left(1+A{C}_i\right)\ln \left(1+A{C}_i\right)},$

после решения которого будет определён параметр B из (18). Заметим, равенство (18) показывает, что $B>0$.

Вариант 4. Для целевой функции (13) НУМ $\frac{\partial S_2}{\partial A}=0$ позволяет выразить параметр A:

$A=\frac{{\displaystyle {\sum }_{i=1}^N\left[\exp \left(\frac{{\sigma }_0-{\sigma }_i}{B}\right)-1\right]C_i}}{{\displaystyle {\sum }_{i=1}^N{C}_i^2}},$                       (20)

из другого НУМ $\frac{\partial S_2}{\partial B}=0$ можно получить соотношение:

${\sum }_{i=1}^N\left[\exp \left(\frac{{\sigma }_0-{\sigma }_i}{B}\right)-1-A{C}_i\right]\exp \left(\frac{{\sigma }_0-{\sigma }_i}{B}\right)({\sigma }_0-{\sigma }_i)=0$,

подставив в которое выражение для параметра A из (20), придём к трансцендентному уравнению относительно параметра B:

${\sum }_{i=1}^N\left[\exp \left(\frac{{\sigma }_0-{\sigma }_i}{B}\right)-1\right]\exp \left(\frac{{\sigma }_0-{\sigma }_i}{B}\right)\left({\sigma }_0-{\sigma }_i\right)$=

$=\frac{{\displaystyle {\sum }_{i=1}^N\left[\exp \left(\frac{{\sigma }_0-{\sigma }_j}{B}\right)-1\right]C_i}}{{\displaystyle {\sum }_{i=1}^N{C}_i^2}}{\sum }_{i=1}^N{C}_i\left({\sigma }_0-{\sigma }_i\right)\exp \left(\frac{{\sigma }_0-{\sigma }_i}{B}\right)$.

После нахождения параметра B параметр A находится из (20).

Вариант 5. Применение НУМ $\frac{\partial S_3}{\partial B}=0$ к целевой функции (14) позволяет получить:

$B=\frac{{\displaystyle {\sum }_{i=1}^N{\left({\sigma }_0-{\sigma }_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^N\left[\ln \left(1+A{C}_i\right)\left({\sigma }_0-{\sigma }_i\right)\right]}}.$                       (21)

Другое НУМ $\frac{\partial S_3}{\partial A}=0$ приводит к соотношению:

$\underset{i=1}{\overset{N}{}}\left[\frac{{\sigma }_0-{\sigma }_i}{B}-\ln \left(1+A{C}_i\right)\right]\frac{C_i}{1+A{C}_i}=0$,

из которого после подстановки B из (21) получаем трансцендентное уравнение относительно параметра A:  ${\sum }_{i=1}^N\frac{\left({\sigma }_0-{\sigma }_i\right)C_i}{1+A{C}_i}=\frac{{\displaystyle {\sum }_{i=1}^N{\left({\sigma }_0-{\sigma }_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^N\left[\ln \left(1+A{C}_i\right)\left({\sigma }_0-{\sigma }_i\right)\right]}}\underset{i=1}{\overset{N}{}}\frac{C_i\ln \left(1+A{C}_i\right)}{1+A{C}_i}$

после решения которого, параметр B находится из (21).

Таким образом, в вариантах 2–5 один из неизвестных параметров – A или B – искомой зависимости (7) находится из трансцендентного уравнения, и после этого по явной формуле в виде алгебраического выражения находится другой. При этом используется не численный поиск минимума некоторой функции, как в (15), а более простой способ – численное решение трансцендентного уравнения, которое получается из НУМ целевой функции.

О вычислении адсорбции

После применения метода наименьших квадратов и определения параметров зависимости уравнения Шишковского, которому подчиняются результаты рассматриваемого эксперимента, можно построить изотермы поверхностного натяжения (рис. 1), т. е. зависимость КПН раствора с ПАВ от концентрации ПАВ в растворе [13–16].

 

Рис. 1. Изотерма поверхностного натяжения

Fig. 1. Surface tension isotherm

 

Простое термодинамическое соотношение между поверхностной концентрацией (адсорбцией), обозначаемой Γ, и изменением КПН с активной концентрацией растворяемого вещества σ было выведено Гиббсом:

$\text{Γ}=-\frac{C}{RT}\frac{d\sigma }{dC}$

где R – универсальная газовая постоянная; T – температура. Тогда величина адсорбции может быть вычислена путём дифференцирования уравнения Шишковского (7):

$\text{Γ}=-\frac{C}{RT}\times \frac{d\sigma }{dC}=\frac{C}{RT}\times \frac{BA}{1+AC}.$ (22)

Изотерма адсорбции изображена на рис. 2 [17, 18].

 

Р

Рис. 2. Изотерма адсорбции

Fig. 2. Adsorption isotherm

 

Для нахождения величины предельной адсорбции с целью описания возможного поведения изотермы необходимо взять предел от адсорбции Γ при концентрации С, стремящейся к бесконечности:

${\text{Γ}}_{\infty }=\underset{C\to \infty }{\lim }\Gamma =\frac{C}{RT}\times \frac{BA}{1+AC}=\frac{B}{RT}.$(23)

Обработка экспериментов

Опытные значения данного раздела взяты из [19]. Эксперименты по определению поверхностного натяжения чистых растворов изопропилового спирта и растворов ПАВ после адсорбции на твёрдых адсорбентах проводились при температуре 23 °С и давлении 775 мм рт. ст. Проведены три эксперимента: один с раствором чистого изопропилового спирта и ещё два с растворами изопропилового спирта с добавлением измельчённого твёрдого адсорбента (активированный уголь; смесь глины с двуокисью кремния). Результаты экспериментов представлены в табл. 1–3.

 

Таблица 1. Экспериментальные значения поверхностного натяжения чистых растворов изопропилового спирта

Table 1. Experimental values of surface tension of isopropyl alcohol pure solutions

i	0	1	2	3	4	5
Ci, моль/л/mole/l	0	0,2	0,3	0,4	0,5	0,6
σi, мН/м/mN/m	72,27	46,51	45,28	43,19	41,88	41,54

 

Таблица 2. Экспериментальные значения поверхностного натяжения растворов изопропилового спирта после адсорбции на угле

Table 2. Experimental values of surface tension of isopropyl alcohol solutions after adsorption on carbon

i	0	1	2	3	4	5
Ci, моль/л/mole/l	0	0,2	0,3	0,4	0,5	0,6
σi, мН/м/mN/m	72,27	48,25	47,63	46,01	45,25	44,35

 

Таблица 3. Экспериментальные значения поверхностного натяжения растворов изопропилового спирта после адсорбции на смеси глины с двуокисью кремния

Table 3. Experimental values of surface tension of isopropyl alcohol solutions after adsorption on a mixture of clay and silicon dioxide

i	0	1	2	3	4	5
Ci, моль/л/mole/l	0	0,2	0,3	0,4	0,5	0,6
σi, мН/м/mN/m	72,27	47,98	46,23	44,86	43,73	43,45

 

Аппроксимация полученных экспериментальных значений КПН реализована путём применения описанного МНК. Параметры уравнения Шишковского в первом эксперименте равны A=889,54; B=4,92; во втором – A=3184,94; B=3,67; и в третьем – A=1481,79; B=4,31.

После определения и аппроксимации полученных значений поверхностного натяжения всех исследуемых растворов строились соответствующие изотермы (рис. 3).

 

Рис. 3. Изотермы поверхностного натяжения экспериментов

Fig. 3. Surface tension isotherms of experiments

 

Результаты расчёта адсорбции и данных для построения изотерм адсорбции по формуле (22) представлены в табл. 4–6.

Изотермы адсорбции растворов приведены на рис. 4. Предельная адсорбция, вычисленная по формуле (23), для рассматриваемых растворов равна 2,0001∙10^–9; 1,4903∙10^–9 и 1,7531∙10^–9, соответственно.

Варианты реализации при аппроксимации соотношением Розина–Раммлера

Рассмотрим процесс седиментации, когда на весы оседают сферические частицы разного размера, но одинаковой плотности [20–23]. Их скорость осаждения определяется законом Стокса. Пусть для кривой накопления имеются экспериментальные данные $\left\{\left(P_{i, e}, t_{i,e}\right)\right\}; i=\mathrm{1,}\dots , N$, где $P_{i, e}$ – вес осадка в момент времени $t_{i,e}$; $t_{i,e}<t_{i+\mathrm{1,}e}$; $P_{i,e}\le P_{i+\mathrm{1,}e}$ (рис. 5).

 

Рис. 4. Изотермы адсорбции исследуемых растворов

Fig. 4. Adsorption isotherms of the studied solutions

 

Рис. 5. Кривая седиментации в моменты времени t₁, t₂, t₃…t_N

Fig. 5. Sedimentation curve at times t₁, t₂, t₃…t_N 

 

Таблица 4. Адсорбции растворов чистого изопропилового спирта

Table 4. Adsorption of pure isopropyl alcohol solutions

i	0	1	2	3	4	5
Ci, моль/л/mole/l	0	0,2	0,3	0,4	0,5	0,6
Гi, кмоль/м2/kmole/m2	0	1,9889∙10–9	1,9926∙10–9	1,9945∙10–9	1,9956∙10–9	1,9964∙10–9

 

Таблица 5. Адсорбции растворов изопропилового спирта после эксперимента на угле

Table 5. Adsorption of isopropyl alcohol solutions after the experiment on carbon

i	0	1	2	3	4	5
Ci, моль/л/mole/l	0	0,14	0,19	0,23	0,27	0,31
Гi, кмоль/м2/kmole/m2	0	1,4869∙10–9	1,4878∙10–9	1,4882∙10–9	1,4886∙10–9	1,4888∙10–9

 

Таблица 6. Адсорбции растворов изопропилового спирта после эксперимента на смеси глины с двуокисью кремния

Table 6. Adsorption of isopropyl alcohol solutions after the experiment on a mixture of clay and silicon dioxide

i	0	1	2	3	4	5
Ci, моль/л/mole/l	0	0,16	0,23	0,30	0,37	0,43
Гi, кмоль/м2/kmole/m2	0	1,746∙10–9	1,748∙10–9	1,749∙10–9	1,750∙10–9	1,750∙10–9

 

Перейдём к весовым долям и относительному времени $\left\{\left(y_i, t_i\right)\right\}$:

$y_i=\frac{P_{i,e}-P_{\mathrm{1,}e}}{P_{N,e}-P_{\mathrm{1,}e}};t_i=t_{i,e}-t_{\mathrm{1,}e}.$ (24)

Для аппроксимации весовых долей (24) используют уравнение Розина–Раммлера:

 $y=1-\exp \left(-b{t}^c\right),$               (25)

$y\left(0\right)=\mathrm{0,}y\left(\infty \right)=\mathrm{1,}$

в котором «b» и «c» – искомые параметры аппроксимации. Тогда целевая функция, сумма квадратов невязок, имеет вид:

$S\left(b,c\right)={\sum }_{i=1}^N{\left[1-\exp \left(-b{t}_i^c\right)-y_i\right]}^2.$(26)

Для уравнения Розина–Раммлера (25) можно привести несколько эквивалентных записей:

 $\ln \left(1-y\right)=-b{t}^c,$                   (27)

 $\ln \left(-\ln \left(1-y\right)\right)=\ln b+c\times \ln t,$         (28)

${\left(1-y\right)}^{1/t^c}=\exp \left(-b\right).$             (29)

Приведённые формы записи (27)–(29) уравнения Розина–Раммлера (25) позволяют привести другие целевые функции МНК:

$S_1\left(b,c\right)=\underset{i=1}{\overset{N}{}}{\left[\ln \left(1-y_i\right)+b{t}_i^c\right]}^2,$      

$S_2\left(x,c\right)=\underset{i=1}{\overset{N}{}}{\left\{\ln \left[-\ln \left(1-y_i\right)\right]-x-c\times \ln \left(t_i\right)\right\}}^2,$

$x=\ln b,$           (30)

 $S_3\left(z,c\right)=\underset{i=1}{\overset{N}{}}{\left[{\left(1-y_i\right)}^{\frac{1}{t_i^c}}-z\right]}^2,$         

$z=\exp \left(-b\right)$.               (31)

Рассмотрим варианты нахождения параметров «b» и «c».

Вариант 1. Методы нелинейного программирования применяются для поиска минимума непосредственно к целевой функции (26) или к аналогам (27)–(29), ищутся значения параметров «b» и «c», при которых этот минимум достигается:

 $\left(b,c\right)=\arg \left\{\min S\left(b,c\right)\right\}.$     (32)

Вариант 2. Из НУМ $\frac{\partial S_1}{\partial b}=0$ получим:

$b=\frac{{\displaystyle {\sum }_{i=1}^N{t}_i^c\text{ln}(1-y_i)}}{{\displaystyle {\sum }_{i=1}^N{\left(t_i^c\right)}^2}}.$(33)

Другое НУМ $\frac{\partial S_1}{\partial c}=0$ приводит к соотношению: $\underset{i=1}{\overset{N}{}}\left[t_i^c\text{ln}{t}_i\text{ln}(1-y_i)\right]+b\underset{i=1}{\overset{N}{}}{\left(t_i^c\right)}^2\text{ln}{t}_i=\mathrm{0,}$

из которого, используя (33), получаем трансцендентное уравнение относительно параметра c:

$\underset{i=1}{\overset{N}{}}{\left(t_i^c\right)}^2·\underset{i=1}{\overset{N}{}}\left[t_i^c\ln t_i\ln \left(1-y_i\right)\right]+$        

$\underset{i=1}{\overset{N}{}}{t}_i^c\ln \left(1-y_i\right)·\underset{i=1}{\overset{N}{}}{\left(t_i^c\right)}^2\ln t_i=\mathrm{0,}$ (34)

решение которого позволяет из (33) найти и параметр b.

Вариант 3. Применение НУМ $\frac{\partial S_2}{\partial x}=0$ к целевой функции (30) позволяет получить: 

$x=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left[\ln \left(-\ln \left(1-y_i\right)\right)-c·\text{ln}{t}_i\right]$ (35)

или:

$b=\exp \left\{\frac{1}{N}\underset{i=1}{\overset{N}{}}\left[\ln \left(-\ln \left(1-y_i\right)\right)-c·\ln t_i\right]\right\}.$(36)

Другое условие минимума $\frac{\partial S_2}{\partial c}=0$ приводит к равенству: $\underset{i=1}{\overset{N}{}}\left[\ln \left(-\ln \left(1-y_i\right)\right)-c·\text{ln}{t}_i\right]·\text{ln}{t}_i-x\underset{i=1}{\overset{N}{}}\ln t_i=\mathrm{0,}$

с помощью (35), исключая параметр x аналогично [24], получаем явное выражение для нахождения параметра c:

$c=\frac{N{{\displaystyle \sum }}_{i=1}^N\left[\ln \left(-\ln \left(1-y_i\right)\right)\right]·\text{ln}{t}_i-{{\displaystyle \sum }}_{i=1}^N\ln t_i·{{\displaystyle \sum }}_{i=1}^N\left[\ln \left(-\ln \left(1-y_i\right)\right)\right]}{N{{\displaystyle \sum }}_{i=1}^N{\left(\text{ln}{t}_i\right)}^2-{\left({\displaystyle {\sum }_{i=1}^N\text{ln}{t}_i}\right)}^2}.$(37)

После нахождения параметра c параметр b также определяется по явной формуле (36).

Вариант 4. Для целевой функции (31) НУМ $\frac{\partial S_3}{\partial z}=0$ позволяет выразить:

$z=\frac{1}{N}\underset{i=1}{\overset{N}{}}{\left(1-y_i\right)}^{\frac{1}{t_i^c}}$   (38)

или:

$b=-\ln \left\{\frac{1}{N}\underset{i=1}{\overset{N}{}}{\left(1-y_i\right)}^{\frac{1}{t_i^c}}\right\}.$      (39)

Второе НУМ $\frac{\partial S_3}{\partial b}=0$ даёт соотношение: $\underset{i=1}{\overset{N}{}}\left[{\left(1-y_i\right)}^{\frac{1}{t_i^c}}-z\right]{\left(1-y_i\right)}^{\frac{1}{t_i^c}}·\text{ln}(1-y_i)\frac{\ln t_i}{t_i^c}=\mathrm{0,}$

подставив в которое z из (38), приходим к трансцендентному уравнению относительно параметра c:

$\underset{i=1}{\overset{N}{}}{\left(1-y_i\right)}^{\frac{2}{t_i^c}}·\text{ln}(1-y_i)\frac{\ln t_i}{t_i^c}$

$\frac{1}{N}\underset{i=1}{\overset{N}{}}{\left(1-y_i\right)}^{\frac{1}{t_i^c}}·\underset{i=1}{\overset{N}{}}{\left(1-y_i\right)}^{\frac{1}{t_i^c}}·\ln \left(1-y_i\right)\frac{\ln t_i}{t_i^c},$(40)

решив которое, параметр b найдём из (39).

Отметим, что для целевой функции (26) ни одно из НУМ не приводит к линейному уравнению относительно какого-либо из искомых параметров, удаётся свести лишь к квадратному уравнению. Поэтому ввиду громоздкости выражений мы не стали получать решение для функции (26). А в варианте 3 оба НУМ приводят к линейным уравнениям относительно обоих искомых параметров, что позволило выписать решение в явном виде. В вариантах 2 и 4 параметр «c» искомой зависимости (25) находится из трансцендентного уравнения (34) или (40), и после этого по явной формуле в виде алгебраического выражения (33) или (39), соответственно, находится другой параметр «b». При этом используется не численный поиск минимума некоторой функции как в (32), а более простой способ – численное решение трансцендентного уравнения, которое получается из НУМ целевой функции.

Решение проблемы с особенностью

Формулы в вариантах 2–4 имеют особенность типа $\ln \left(t_i\right)$, или $\ln \left(1-y_i\right)$, или $\ln \left[-\ln \left(1-y_i\right)\right]$ при $y_i=\mathrm{0,} y_i=\mathrm{1,} t_i=0$, что во время счёта приведёт к аварийному останову. Для преодоления этой проблемы предлагается алгоритм:

• Выбирается вариант для расчёта параметров b и c, которые вычисляются методом (32). Это будет образец для сравнения.
• Выбирается малое число $\epsilon >0$, например, $\epsilon ={10}^{-3}$.
• Значения 0 и 1 в данных (24) заменяются на $\epsilon$ и $1-\epsilon$.
• По новым данным проводятся вычисления для выбранного варианта.
• Значение $\epsilon$ уменьшается, $\epsilon =\epsilon /2$.
• Значения 0 и 1 в данных (24) заменяются на $\epsilon$ и $1-\epsilon$.
• По новым данным проводятся повторные вычисления.
• Точность вычислений определяется путём сравнения найденных параметров b и c для $\epsilon$ и $\epsilon /2$.
• Если точность вычислений не удовлетворяет, то $\epsilon$ уменьшается и производится переход на пункт 3.

Для проверки алгоритма этого пункта рассмотрен модельный пример

$\left\{\left(y_i,t_i\right)\right\};i=\mathrm{1,}\dots ,14$;

где

$t_1=\mathrm{0,}{t}_{i+1}=t_i+\mathrm{10,}{t}_{14}=130;$ $y_i=y\left(t_i\right)·(1+\mathrm{0,01}·RNd)$ ,

где $RNd\in \left[\mathrm{0,}1\right]$ – случайное число, $y\left(t_i\right)$ соответствует формуле (25) с коэффициентом $b=\mathrm{0,2}$ и $c=\mathrm{0,2}$. Расчеты показали, что при $\epsilon <{10}^{-3}$ ошибки в найденных коэффициентах меньше 1 %.

О дифференциальном распределении

Дифференциальное распределение в зависимости от радиуса частицы r получается с использованием формулы Стокса и перехода от рассмотрения зависимости от времени оседания к зависимости от радиуса частицы. Имеем:

$\left({\rho }_{\text{част}}-{\rho }_{\text{ж}}\right)·V·g=6\pi r·w·\eta ,$                 (41)

$w=\frac{L}{t},$                           (42)

здесь ${\rho }_{част}$ и ${\rho }_{ж}$ – плотность материала частиц и жидкости; V – объем частицы; g – ускорение свободного падения; w – скорость частицы; $\eta$ – вязкость жидкости; L – высота столба жидкости.

Из (41) и (42) следует:

$t=\frac{k^2}{r^2}$                            (43)

$k=\sqrt{\frac{9·L·\eta }{2g\left({\rho }_{част}-{\rho }_{ж}\right)}}$,

где k – постоянная эксперимента (константа Стокса). Тогда интегральная кривая распределения по размерам частиц получается как зависимость относительного веса частиц дисперсной фазы от радиуса при подстановке (43) в (25):

$Q\left(r\right)=y\left(t\left(r\right)\right)=1-\exp \left(-b\frac{k^{2c}}{r^{2c}}\right).$          (44)

Имеем:

$\frac{dt}{dr}=-k^2·r^{-3},$                            (45)

тогда из (25), (43), (45) получается дифференциальное распределение в зависимости от r :

$\frac{dQ}{dr}=\frac{dy}{dt}\frac{dt}{dr}=-\exp \left(-b\frac{k^{2c}}{r^{2c}}\right)·2bc·\frac{k^{2c}}{r^{2c+1}}.$ (46)

Найдём $r_{*}$ – радиус, которому соответствует максимальная доля частиц:

$r_{\text{*}}=\underset{r}{\arg }\left\{\max \left(\frac{dQ}{dr}\right)\right\}=$

$={\underset{}{\arg }}_r\left(\frac{d^{2}Q}{d{r}^2}=0\right)=k{\left(\frac{2bc}{2c+1}\right)}^{\frac{1}{2c}}.$                        (47)

Тогда максимальная доля частиц будет:

${\left.\frac{dQ}{dr}\right|}_{r_{\text{*}}}=\max \left(\frac{dQ}{dr}\right)=\exp \left(-1-\frac{1}{2c}\right)\times \frac{2c+1}{k}{\left(\frac{2c+1}{2bc}\right)}^{\frac{1}{2c}}.$(48)

Анализ седиментации суспензии

Проведём седиментационную обработку результатов эксперимента по данным табл. 1.2.1 из [3], где использовалась суспензия оксида алюминия А1₂O₃ в дистиллированной воде. Для получения кривой осаждения суспензии определяли изменение веса чашечки, помещённой в суспензию, по мере накопления на ней осадка с помощью торсионных весов ВТ-500.

 

Таблица 7. Экспериментальные значения накопленного веса от времени осаждения

Table 7. Experimental values of accumulated weight versus deposition time

i	1	2	3	4	5	6	7	8	9	10	11	12	13	14
ti,e, мин/min	0,25	0,5	0,75	1	1,5	2	2,5	3	5	7	10	15	30	45
Рi,е, мг/mg	31,5	50,5	59,5	65,5	68,5	72,5	76,5	79,5	80,5	81,5	81,5	81,5	82,5	82,5

 

Плотность оксида алюминия – 3,99 г/см³, высота столба суспензии до измерительной чашечки – 0,2 м. Эксперимент продолжается до тех пор, пока измерения веса осадка за 5 минут принимают очень близкие значения. Данные эксперимента представлены в табл. 7.

В соответствии с формулой (24) переходим к весовым долям и относительному времени. Находим параметры зависимости Розина–Раммлера по варианту 3 МНК с явными формулами нахождения параметров (37) и (36). Однако $y_{13}=y_{14}=1$, непосредственное применение указанных формул приведёт к аварийному останову, поэтому, следуя алгоритму, предложенному для решения проблем с особенностью, положим $y_{13}=1-{10}^{-4}$, $y_{14}=1-{10}^{-5}$, что позволяет обойти эту проблему.

На рис. 6 представлены экспериментальные значения в весовых долях и кривая осаждения, представленная зависимостью Розина–Раммлера, аппроксимирующей экспериментальные значения по МНК. Найденные параметры таковы b=1,2138, c=0,576.

 

Рис. 6. Кривая седиментации из МНК и экспериментальные данные

Fig. 6. Sedimentation curve from the least squares method (LSM) and experimental data

 

Результаты седиментационного анализа, отражающие распределение частиц по размерам в системе, представляют в виде интегральной и дифференциальной функций (рис. 7, 8), вычисленных по формулам (44) и (46) соответственно. Максимальная доля частиц рассчитывается по (48).

 

Рис. 7. Интегральная кривая распределения частиц по размерам

Fig. 7. Integral particle size distribution curve

 

Рис. 8. Дифференциальная кривая распределения частиц по размерам

Fig. 8. Differential particle size distribution curve

 

Радиус, которому соответствует максимальная доля частиц эксперимента, определяем по формуле (47), он равен $r_{*}=1.14·{10}^{-4}$.

Заключение

Приведены варианты формул, полученных на основе МНК, позволяющие аппроксимировать экспериментальные данные соотношениями Шишковского и Розина–Раммлера с использованием как известных способов на основе нелинейного программирования, так и способов на основе решения трансцендентного уравнения. Приведены примеры применения полученных формул на данных реальных экспериментов.

Для известных алгебраических формул, которые позволяют найти решение безитерационно, но имеют особенность, предложен алгоритм обхода этой особенности.

作者简介

Vladislav Galkin

National Research Tomsk Polytechnic University

Email: vlg@tpu.ru

Associate Professor

俄罗斯联邦, Tomsk

Yuriy Volkov

Sobolev Institute of Mathematics

Email: volkov@math.nsc.ru
ORCID iD: 0000-0002-7298-8578

Dr. Sc., Chief Researcher

俄罗斯联邦, Novosibirsk

Liliya Chekantseva

National Research Tomsk Polytechnic University

编辑信件的主要联系方式.
Email: chlb@tpu.ru

Senior Lecturer

俄罗斯联邦, Tomsk

Vladimir Ivanov

National Research Tomsk Polytechnic University

Email: vai14@tpu.ru

Master's Student

俄罗斯联邦, Tomsk

参考

补充文件