ВведениеПолнота описания технологических процессовпищевых производств как сложных системопределяется степенью детализации свойств изакономерностей всех подсистем и связей между ними.Дополнение новых подсистем, уточнение их связис другими подсистемами, а также оценка влиянияна эффективность функционирования всей системыпозволяет выявить ее новые свойства, повыситькачество оценки ее состояния, точность прогнозови эффективность управления технологическимпроцессом как системой в целом.Глубина и интенсивность протекания отде-льных процессов, а также показатели качестваполуфабрикатов и готовой продукции зависятот реологических свойств, главным из которыхявляется вязкость. Пищевые массы представляютсобой структурированные дисперсные системы,обладающие аномалией вязкости. Реологическиесвойства пищевых масс определяют выбор способа ихдальнейшей переработки, способность к сохранениюили изменению сформированной структуры [1–3].Реологические свойства пищевого сырья иполуфабрикатов влияют на термогидродинамикупроцессов их перемешивания, транспортированиепо трубам, формование, нагрев, охлаждение и пр.При проектировании оборудования пищевыхпроизводств используются расчетные методики,базирующиеся на следующих допущениях: изо-термичность и ламинарность режима, осреднениеосновных параметров процесса (скорость, вязкость,температура, концентрация). При этом не учитываютсянестационарность реологических и теплофизическихсвойств, а также их пространственная неодно-родность [3–5].Сделанные допущения оказываются справедливыдля многих процессов, обладающих высокойинерционностью, интенсивным перемешиванием ипротекающих в аппаратах, допускающих упрощениедо простейших геометрических форм (например,1 Voronezh State Technical University , Voronezh, Russia2 Military Educational and Scientific Centre of the Air Force N.E . Zhukovskyand Y.A. Gagarin Air Force Academy, Voronezh, Russia3 Voronezh State University of Engineering Technologies , Voronezh, Russia4 Voronezh State University , Voronezh, RussiaReceived: April 21, 2021 Accepted in revised form: May 13, 2021Accepted for publication: July 15, 2021*е-mail: zhuraa1@rambler.ru© A.A. Khvostov, G.O. Magomedov, V.I. Ryazhskih,A.V. Kovalev, A.A. Zhuravlev, M.G. Magomedov, 2021Abstract.Introduction. Carreau’s rheological model can describe the three-dimensional flows of non-Newtonian media. However, itrequires modeling parameters for the viscosity of the medium at the limiting values of shear rates, which cannot be achievedby instrumental methods. The present article introduces a novel method that can identify the parameters of Carreau’s modelusing a regularization algorithm.Study objects and methods. The study featured fondant mass produced according to the traditional formulation for CreamyFondant unglazed candies. Standard methods were used to describe the properties of the raw materials and semi-finishedproducts, as well as methods of mathematical processing, modeli ng, and optimization.Results and its discussion. The research produced an algorithm based on A.N. Tikhonov’s regularization method of the parametricidentification of Carreau’s rheological model. The calculation residual was minimized by the viscometric measurements andthe CFD model, which provided the calculation of the hydrodynamic flow regime at the limiting values of shear rates. TheCFD model of a steady non-isothermal flow of a nonlinear viscou s medium through a cylindrical capillary was based on theequations of conservation of mass, energy, and momentum. The rheological parameters of Carreau’s model were illustratedby the case of fondant mass. The error for the viscosity predic tion did not exceed 14.07%.Conclusion. The parametric identification algorithm made it possible to evaluate the rheological parameters of structured liquidmedia with Carreau’s rheological law in cases that lack experimental information on the behavior of the medium at limitingshear rates. The algorithm eliminated the computational problems typical of Ostwald and de Ville’s rheological model, whichusually arise when solving practical problems of three-dimensional flows of non-Newtonian media with limiting viscosity values.Keywords. Regularization, identification, rheological model, Carreau flui d model, hydrodynamics, CFD-modelFor citation: Khvostov AA, Magomedov GO, Ryazhskih VI, Kovalev AV, Zhuravlev AA, Magomedov MG. Carreau’sRheological Model and A.N. Tikhonov’s Regularization Method: Parametric Identification Based on a CFD model. FoodProcessing: Techniques and Technology. 2021;51(3):615–627. (In Russ.). https://doi.org/10.21603/2074-9414-2021-3-615-627.617Хвостов А. А. [и др.] Техника и технология пищевых производств. 2021. Т. 51. № 3 С. 615–627цилиндрические и плоскощелевые каналы). Однакопри расчетах низкоинерционных процессов ваппаратах со сложной геометрией (спиральные каналы,каналы с застойными зонами, области резких суженийи расширений и пр.) необходимы нестационарныепространственные модели, учитывающие существен-ную геометрическую неоднородность полейскоростей, давлений, температур и скоростейсдвига [6, 7]. В этой связи гидромеханическиепроцессы необходимо рассматривать сопряженнос процессами теплообмена, использовать моделинеизотермического течения сплошной среды, а такжезависимость теплофизических свойств среды оттемпературы, а вязкость – как нелинейную функциютемпературы и скорости сдвига [8, 9].Основой для гидродинамических расчетовслужит реологическая модель сплошной среды,создание которой возможно в рамках структурно-параметрического синтеза. Это требует проведенияструктурной идентификация реологической моделии оценки ее параметров.Для описания реологического поведенияструктурированных жидкообразных пищевыхсред широко используется уравнение Оствальда-де Виля, устанавливающее связь между вязкостьюμ (Па·с) и скоростью сдвига γ (с–1) в виде степенногоуравнения [1, 2]: 11 **nK nγμ γ μγ−−  = =    (1)где K – консистентная переменная, Па·сn; n – индекстечения; γ* – скорость сдвига приведения, с–1,μ* – вязкость при скорости сдвига приведения, Па·с.Несмотря на то, что уравнение Оствальда-де Виля не имеет теоретического обоснованияи является всего лишь удачным эмпирическимприближением, степенная зависимость (1) показаласвою состоятельность при описании реологическихкривых μ (γ ) широкого спектра пищевых масс,расплавов полимеров, биологических сред, продуктовнефтепереработки и др. [1–5, 10, 11].Однако уравнение Оствальда-де Виля имеетсущественные недостатки. Во-первых, прииспользовании степенного уравнения (1) возникаютсущественные вычислительные проблемы прирасчетах в окрестности предельных значенийскорости сдвига γ →0 и γ →∞. В этом случае модельнекорректна. Нестационарные пространственныемодели течения сплошной среды используютполе скоростей, распределенное по сечению вдиапазоне от нулевого значения на стенке донекоторого максимального значения. При этом вточке экстремума первая производная скорости попространственной координате (скорость сдвига) можетобращаться в ноль, что приводит, в соответствиис (1), к сингулярности ввиду бесконечного значениявязкости ( ( )0limγμ γ→= ∞ ). С другой стороны, в отсутствиипроскальзывания на стенке канала скорость сдвигастремится к бесконечности. При этом вязкостьобращается в ноль ( lim ( ) 0γμ γ→∞= ).Бесконечная вязкость при нулевой скорости сдвигаприводит к ошибочному результату, когда в расчетнойобласти встречается область нулевой скорости сдвига.Например, при течении в канале любой геометрии привыполнении условий прилипания скорость течениясреды будет изменяться от нуля на стенке канала,проходить через экстремум и снижаться до нуля напротивоположной стенке. Это приведет к появлениюнулевой скорости сдвига и бесконечной вязкости вточке экстремума. При этом замена нулевых значенийсдвиговой вязкости в точке экстремума некоторымиконечными значениями приводит к искажениюрасчетных значений вязкости в этой области.Предсказанное распределение скоростей в этом случаебудет более плоским в центре, чем экспериментальныйпрофиль [11]. Другим примером является процессзаполнение формы нелинейно-вязкой средой приее литье под давлением. Отливка сопровождаетсявысокими скоростями потока в начале процессазаполнения и низкими значениями скорости потокав момент его завершения. Здесь прогнозируемоезначение вязкости по модели Оствальда-де Виля (1)в конце заполнения будет слишком высоким. Этоприведет к завышенным значениям необходимогодавления отливки.Вторая проблема заключается в невозможностиописания полной реологической кривой μ (γ ) однимуравнением (1), что имеет место при изменении γ нанесколько десятичных порядков. Проблема решаетсяразбиением реологической кривой μ (γ ) на несколькоотдельных участков, в пределах которых соблюдаетсялинейная зависимость ln μ от lnγ, которая достигаетсяподбором значения скорости сдвига приведения γ*.Согласно основным положениям физико-химической механики дисперсных систем,сформулированными академика П. А. Ребиндером,структурированные жидкообразные средыхарактеризуются двумя уровнями постоянной(ньютоновской) вязкости [1, 2]. Один из них0 μ– наибольшая вязкость практически неразрушеннойструктуры при γ →0 (первое ньютоновское плато);при скорости сдвига γ →∞ структурированнаяжидкообразная среда обладает наименьшей вязкостьюпредельно разрушенной структуры μ∞ (второеньютоновское плато). В переходной области (областилавинного разрушения структуры) ( ) 0μ μ γ μ∞≥  ≥ .Преодоление вычислительных проблем,связанных с получением неустойчивого решенияпри использовании уравнения Оствальда-де Виля,достигается путем его регуляризации за счет введенияв уравнение (1) дополнительных параметров:предельных вязкостей 0μ и μ∞ . На сегодняпредложено несколько реологических уравнений,учитывающих конечные значения вязкости припредельных значениях скорости сдвига.Трехпараметрические уравнения Штейнера,Ферри, де Хавена и Сиско, учитывающие одну из618Khvostov A.A. et al. Food Processing: Techniques and Technology, 2021, vol. 51, no. 3, pp. 615–627двух предельных вязкостей – 0μ или μ∞, являютсяматематическим описанием неполных реологическихкривых структурированных жидкообразных сред [12,13]. Четырехпараметрические уравнения Кросса,Карро и его модификация (Карро-Яшида) имеютсходную между собой структуру и обеспечиваютвысокую точность описания полных реологическихкривых нелинейно-вязких сред за счет одновременногоучета двух ньютоновских вязкостей 0μ и μ∞ [2,12–14].Наиболее распространенной является модельКарро, коэффициенты которой имеют некотороетеоретическое обоснование:( ) ( )12 20 1nμ μ μ μ λγ−∞ ∞  = + − +  (2)где 0μ – наибольшая вязкость практическинеразрушенной структуры, Па·с; μ∞ – наименьшаявязкость предельно разрушенной структуры, Па·с;γ – скорость сдвига, с–1, n – показатель нелинейностиреологической кривой μ (γ ); λ – время релаксации, с.Модель Карро учитывает предельныеньютоновские состояния структурированнойжидкообразной среды ( ) 00limγμ γ μ→= и lim ( )γμ γ μ∞→∞= иописывает ее нелинейно-вязкое поведение в интервалескоростей сдвига от γ →0 до γ →∞. Степеньотклонения от ньютоновского характера течениясреды характеризуется показателем нелинейности n,смысл которого тождественен индексу теченияуравнения Оствальда-де Виля (1). Обратной величинойк времени релаксации λ является критическая скоростьсдвига, при которой происходит резкое снижениевязкости, связанное с лавинообразным разрушениемструктуры при сдвиговом течении [2].Модель Карро (2) может быть использована дляустранения вычислительных проблем, возникающихпри решении широкого круга задач ламинарныхпространственных течений нелинейно-вязких средв изотермической и неизотермической постановках.Однако, ввиду существенной нелинейности моделиКарро, получение аналитического решения такихзадач, даже с учетом общепринятых допущений(ламинарный, установившийся, изотермическийрежим течения, пренебрежение сжимаемостьюжидкости, инерционными и массовыми силами),не представляется возможным.Известны также реологические уравнения,включающие гиперболические и экспоненциальныефункции скорости и напряжения сдвига (уравненияПрандтля-Эйринга, Пауэлла-Эйринга, Рейнера,Михайлова-Лихтхайма и пр.), а также уравнения,содержащие шесть и более коэффициентов. Попричине редкого использования таких уравненийони в данной работе не рассматриваются.Основным методом определения вязкостиявляется сдвиговая вискозиметрия (ротационная икапиллярная), широко используемая в производствен-ных условиях и научно-исследовательской практикедля изучения вязкостных свойств ньютоновских иненьютоновских пищевых и биологических сред,сырья и продуктов нефтехимии и пр. [1–3, 10–12].Серийно выпускаемые ротационные вискозиметры,укомплектованные различными измерительнымисистемами (коаксиальные цилиндры, конус-плоскость, конус-конус и пр.) с возможностьюреализации режимов измерения с контролируемымнапряжением или скоростью сдвига (CS- илиCR-режим соответственно), позволяют определятьвязкость в интервале низких и средних скоростейсдвига от 10–3 до 102 с–1. Использование в капиллярныхвискозиметрах капилляров разного диаметра и длины,а также измерения в режимах постоянного расходаили давления делают возможным изучение вязкостныхсвойств в интервале средних и высоких скоростейсдвига от 10–1 до 105 с–1.Существующие методы сдвиговой вискозиметриии их приборное оформление не позволяют оценитьвязкость в предельных случаях при γ →0 и γ →∞,что объясняется рядом причин.Верхняя граница скорости сдвига определяетсяследующими эффектами:1) переход от ламинарного к турбулентному режимутечения и возникновение неустойчивости теченияпотока материала в измерительной системе;2) значительные тепловыделения в слое испытуемогоматериала вследствие диссипации механическойработы при высоких скоростях сдвига иневозможность обеспечения изотермичности течения;3) появление эффекта Вайсенберга.Невозможность проведения измерений наротационных и капиллярных вискозиметрах вусловиях γ →0 объясняется следующими причинами:1) весьма существенная продолжительностьреологических измерений (несколько десятковминут), с которой связана стабильность (как физико-химическая, так и структурная) исследуемогоматериала в условиях эксперимента;2) предъявление высоких, зачастую невыполнимых,требований к техническим средствам задания иизмерения малых линейных и угловых перемещенийкрутящего момента за длительное время;3) наличие трудно поддающихся устранениюпаразитных сопротивлений (трение в подшипниковыхопорах и механических передачах подвижныхчастей вискозиметра), которые вносят в результатыизмерений систематические ошибки.Частичному разрешению проблемы измеренияпредельной вязкости при γ →0 способствуетприменение в реометрической практике сдвиговыхпластометров с реализацией испытаний на ползу-честь [3]. Однако следует с осторожностью относитьсяк результатам, получаемым на пластометрах попричине отсутствия теоретического обоснованияприборной инвариантности методов ротационной(капиллярной) вискозиметрии и методов испытанияна ползучесть. При этом возникает необходимостьв статистической проверке корреляции между619Хвостов А. А. [и др.] Техника и технология пищевых производств. 2021. Т. 51. № 3 С. 615–627результатами вискозиметрии и испытаний напластометрах.Другой применяемый на практике подход –экстраполяция результатов вискозиметрии на уровниγ →0 и γ →∞ . Недостатком экстраполяционногоподхода к определению значений предельныхвязкостей является высокая доля субъективизма,связанного с интуицией и предпочтениямиэкспериментатора, приводящая к значительнымпогрешностям предсказания вязкости в предельныхслучаях.Таким образом, методы и приборное оформлениесдвиговой вискозиметрии неньютоновских средне дают информации о полной реологическойкривой μ (γ ) и не позволяют провести идентификациюреологических параметров модели Карро (2).Кроме того, общепринятые методы обработкивискозиметрических данных дают информациюо вязкости среды по осредненным по объемуизмерительной ячейки значениям скорости илинапряжения сдвига.Решение проблемы видится в уточнениирезультатов натурных вискозиметрическихэкспериментов посредством дополнительнойинформации о гидродинамическом поведенииматериала в измерительной системе вискозиметра.Источником такой информации может служитьвычислительный реологический экспериментс привлечением методов вычислительнойгидродинамики или CFD-подхода (Computationalfluid dynamics). Суть CFD-подхода заключается вчисленном решении для выбранной реологическоймодели уравнений неразрывности, сохраненияимпульса и энергии при соответствующих начальныхи граничных условиях [15]. CFD-модель даетпространственную картину потока продукта ввискозиметрической системе ротационного иликапиллярного типа с учетом большого числаэффектов (тепловых, пристенных, входовых, кинети-ческих пр.), возникающих в ходе натурногореологического эксперимента, учет которых вреальном эксперименте не всегда возможен.Для эффективного применения CFD-модели-рования имеется ряд коммерческих программныхпродуктов, таких как ANSYS, FlowVision, COMSOLMultiphysics, и программ с открытым кодом.Например, OpenFOAM, Salome и Code Saturn.В этой связи сдвиговая вискозиметриякак физический метод определения вязкостипереходит на качественно новый уровень –приборно-вычислительную вискозиметрию.Последняя представляет собой совокупностьметодов экспериментальной вискозиметрии,вычислительных экспериментов на CFD-моделях,а также математического аппарата для структурногои параметрического анализа и синтеза реологическихсвойств и моделей.Целью работы явилась разработка математи-ческого метода идентификации параметровреологической модели Карро с использованиемалгоритма регуляризации по результатамсдвиговой вискозиметрии и вычислительногоэксперимента, проведенного на CFD-моделитечения структурированной жидкообразной средыв измерительной системе вискозиметра.Объекты и методы исследованияОбъектом исследований явилась помаднаяконфетная масса, приготовленная традиционнымспособом по рецептуре неглазированных помадныхконфет белого сорта «Сливочная помадка».Для приготовления образцов применяли сырье,соответствующее требованиям нормативных доку-ментов: сахарный песок (ГОСТ 33222-2015), патока(ГОСТ 33917-2016), молоко сгущенное (ГОСТ 31688-2012), масло сливочное (ГОСТ 32261-2013), ванилин(ГОСТ 16599-71).В работе использовали органолептические, физико-химические и вискозиметрические методы анализа.Органолептические показатели помадной конфетноймассы определяли по ГОСТ 5897-90; массовую долювлаги – по ГОСТ 5900-2014; органолептические ифизико-химические показатели помадной конфетноймассы соответствовали ГОСТ 4570-2014.Ротационную вискозиметрию конфетной массыпроводили на ротационном вискозиметре Rhеоtest RN4.1 с использованием цилиндрической измерительнойсистемы S3 в режиме CR. Диапазон измененияскорости сдвига от 0,47 до 20 с–1. Обработку данныхпроводили согласно ГОСТ 1929-87.Капиллярную вискозиметрию осуществлялина автоматическом капиллярном вискозиметреАКВ-2ЖВ. Для исключения входовых потерь давленияиспользовали два капилляра одного диаметра(0,0032 м) и разной длины (0,1 и 0,15 м). Обработкуданных проводили согласно ГОСТ 7163-84.В течение всей продолжительности экспериментаконфетную массу термостатировали при температуре60 °С. Для этого использовали циркуляционныйтермостат Thermovisc B100F8TFT.Предметом исследований явилась реологическаямодель помадной конфетной массы, а также мате-матический метод ее параметрической иденти-фикации с использованием алгоритма регуляризацииА. Н. Тихонова на основе CFD-модели течениянелинейно-вязкой жидкообразной среды визмерительной системе вискозиметра.В работе использовали методы и программныесредства для математической обработки экспери-ментальных данных, моделирования и оптимизации.Для обработки результатов натурных реологическихэкспериментов и построения графических зави-симостей использовали программу MS Excel2010 и математический пакет Mathcad 15.CFD-моделирование течения помадной массы визмерительной системе вискозиметра проведено сиспользованием пакета конечно-элементного анализаANSYS.620Khvostov A.A. et al. Food Processing: Techniques and Technology, 2021, vol. 51, no. 3, pp. 615–627Результаты и их обсуждениеРезультаты ротационной вискозиметриипредставлены на рисунке 1 в виде зависимостивязкости от скорости сдвига (графическиепиктограммы). На рисунке 2 представлены резуль-таты (графические пиктограммы) капиллярнойвискозиметрии в виде зависимостей перепада давленияна концах капилляра от объемного расхода конфетноймассы для двух капилляров с разным отношениемдлины к диаметру L D после исключения входовыхпотерь давления согласно методике [11].Результаты ротационной и капиллярнойвискозиметрии указывают на то, что конфетнаямасса при сдвиговом течении проявляет аномалиювязкости. Это выражается в непропорциональномуменьшении вязкости при увеличении скоростисдвига и является псевдопластичной средой.Поскольку ротационная и капиллярнаявискозиметрия не позволяют провести натурныйэксперимент при предельных значениях скоростисдвига γ →0 и γ →∞ , что не позволяет оценитьпредельные значения вязкости 0μ и μ∞ , то дляпараметрической идентификации реологическоймодели Карро (2) предлагается использоватьалгоритм регуляризации Тихонова, нашедшийпрактическое применение при решении обратныхзадач теплофизики, нефтегазодобычи, спектрометрии,геофизики, биологии и пр. [16–20].Для случая дискретной задачи идентификациипараметров реологической модели Карро (2)регуляризация заключается в минимизациипараметрической функции [16]:( ) ( ) ( ) 1 2 F F F min α = + α → β β β β (3)где ( ) 0μ ,μ ,λ,n ∞ β = – вектор реологических параметровмодели (2); ( ) 1 F β – регуляризующая функция;( ) 2 F β – стабилизирующая (корректирующая) функция;α – параметр регуляризации (α > 0), контролирующийвклад стабилизирующей функции в минимизациюкритерия оптимальности.Минимизация параметрической функции (3)обеспечивает нахождение регуляризованногорешения μ (γ,β) в виде модели Карро (2) порезультатам ротационной вискозиметрии содновременной коррекцией значений векторареологических параметров ( ) 0μ ,μ ,λ,n ∞ β = поCFD-модели капиллярной вискозиметрии. В качествекорректирующего параметра, входящего в ( ) 2 F β ,предлагается использовать перепад давления вкапилляре, являющийся интегральным показателемпотокового течения и поддающийся прямомуизмерению.Регуляризующая ( ) 1 F β и стабилизирующая ( ) 2 F βфункции представлены в виде суммы квадратовотклонений расчетных значений от соответствующихэкспериментальных:( ) ( ) ( ) 211,nýi i i iiF μ γ μ γ=β =Σ  β −   ; (4)( ) ( ) ( ) 221,mýj j j jjF P Q P Q=β =ΣΔ β − Δ  (5)где i, j – номер измерения при ротационной икапиллярной вискозиметрии; n, m – количествоизмерений при ротационной и капиллярнойвискозиметрии; iγ – скорость сдвига в i-омэксперименте; ( ) ýii μ γ – экспериментальное значениевязкости, полученное в i-ом эксперименте; i ( i , ) μ γ β –расчетное значение вязкости, полученное по моделиКарро (2) при скорости сдвига iγ ; j Q – объемныйрасход среды через капилляр в j-ом эксперименте;ý ( )j j ΔP Q – экспериментальное значение перепададавления, полученное в j-ом эксперименте; j ( j , ) ΔP Q β –расчетное значение перепада давления, полученноепо CFD-модели при объемном расходе среды j Q .01503004506007500 3 6 9 12 15 1801002003004005006007008000,0E+00 2,0E-07 4,0E-07 6,0E-07 8,0E-07 1,0E-06L/D = 31,25 L/D = 46,872000300040005000600070008000900010000100200300400500Вязкость, Па·сСкорость сдвига, с–1Рисунок 1. Результаты ротационной вискозиметрииFigure 1. Rotational viscometry results01503004506007500 3 6 9 12 15 1801002003004005006007008000,0E+00 2,0E-07 4,0E-07 6,0E-07 8,0E-07 1,0E-06L/D = 31,25 L/D = 46,872000300040005000600070008000900010000100200300400500Рисунок 2. Результаты капиллярной вискозиметрииFigure 2. Capillary viscometry resultsПерепад давления, ПаОбъемный расход, м3/с1503004506007500100200300400500600700800L/D = 31,25 L/D= 46,8701503004506007500 3 6 9 12 15 1801002003004005006007008000,0E+00 2,0E-07 4,0E-07 6,0E-07 8,0E-07 1,0E-06L/D = 31,25 L/D = 46,87100020003000400050006000700080009000100000,0001 0,001 0,01 0,1 1 1001002003004005000,0001 0,001 0,01 0,1 1 1010010000,30,4L L621Хвостов А. А. [и др.] Техника и технология пищевых производств. 2021. Т. 51. № 3 С. 615–627Связь между реологической кривой μ (γ ) ввиде уравнения Карро (2) и перепадом давления вкапилляре может быть формализована с помощьюCFD-модели, реализованной в вычислительномпакете конечно-элементного анализа ANSYS [15].С учетом (4) и (5) имеем критерий минимизациипараметрической функции:( ) ( ) ( ) ( ) ( ) 2 21 1, , minn mý ýi i i i j j j ji jF P Q P Q α μ γ μ γ α= == Σ  −   + ΣΔ − Δ  → β β β β( ) ( ) ( ) ( ) 2 21 1, , minmý ýi i i i j j j jjP Q P Q μ γ μ γ α=Σ  −   + ΣΔ − Δ  → β β β (6)Вектор реологических параметров β, достав-ляющий минимум параметрической функции (6),зависит от параметра регуляризации α, отражающегобаланс между предпочтительным выбором результатовротационной вискозиметрии и коррекциейрезультатов по CFD-модели. При α ≈ 0 минимизацияпараметрической функции (6) эквивалентнаминимизации функции (4) по данным ротационнойвискозиметрии без учета результатов капиллярнойвискозиметрии. Когда параметр регуляризации αслишком велик, минимизация параметрическойфункции (6) эквивалентна поиску минимумафункции (5). Получаемые в этом случае значениявектора реологических параметров β обеспечиваютрешение, наиболее точно предсказывающеерезультаты капиллярной вискозиметрии без учетаданных, полученных на ротационном вискозиметре.Таким образом, оптимальному значению параметрарегуляризации α * соответствует вектор реологическихпараметров модели Карро (2) ( ) * 0, , ,n α μ μ λ ∞ β = ,при котором одновременно достигается удовлетво-рительная аппроксимация результатов ротационнойи капиллярной вискозиметрии.Параметр регуляризации может быть выбранна основании принципов невязки, максимумаправдоподобия, методов отношений, перекрестнойзначимости (cross-validation), по минимуму функциичувствительности [16, 18, 21]. В работе для выборапараметра регуляризации α предлагается использоватьметод L-кривой [17, 22]. Под ней понимаютпараметрически заданную кривую от параметра α ,определяемую соотношениями:( ) ( ( )) 1 x ln F α α = β , ( ) ( ( )) 2 y ln F α α = β (7)где α β – вектор реологических параметров дляфиксированного значения параметра регуляри-зации α, рассчитанный из условия минимизации (6).Метод L-кривой заключается в нахождениизначения параметра регуляризации α *, прикотором достигается максимальная кривизна k (α )параметрически заданной L-кривой:( ) ( ) ( ) ( ) ( )( ( ) ( ) )2 2 3maxx y x ykx yαα α α ααα α′ ′′ − ′′ ′= →′ + ′(8)где x′(α ), y′(α ), x′′(α ), y′′(α ) – первая и втораяпроизводная (7) по переменной α .Для проведения параметрической идентификацияреологической модели Карро предлагается алгоритмрегуляризации на основе результатов сдвиговойвискозиметрии и вычислительного эксперимента,проведенного на CFD-модели потока среды.На первом этапе для ряда заданных значенийпараметра регуляризации α осуществляют поискзначений вектора реологических параметров α β ,минимизирующего функцию (6). Для этого могутбыть использованы известные алгоритмы численноймногомерной оптимизации [23].На втором этапе осуществляют выбороптимального значения параметра регуляризацииα * по методу L-кривой. Для этого используют условиемаксимизации критерия (8).При CFD-моделировании в вычислительномпакете конечно-элементного анализа ANSYS былаиспользована методика, включающая следующиеэтапы: создание геометрической модели расчетнойобласти моделирования; математическая постановказадачи с выбором физического интерфейса иуправляющих уравнений; определение физическихсвойств частей объекта; задание начальных играничных условий для различных элементовобъекта; задание параметров и генерация конечно-элементной сетки; выбор и задание параметроврешающего устройства; проведение вычислительногоэксперимента, визуализация и анализ результатов [15].При составлении математической модели течениянелинейно-вязкой среды по цилиндрическомукапилляру длиной L и диаметром D принятыследующие допущения: движение среды устано-вившееся, осесимметричное, ламинарное, среданесжимаема, источники тепла отсутствуют,гравитационные силы и диссипативные члены неучитываются. Расчетная область моделированиязадана границами: Ω1 – входное сечение капилляра;Ω2 – цилиндрическая поверхность капилляра (стенка);Ω3 – выходное сечение капилляра.Уравнения сохранения массы, энергии и импульсапри заданных условиях и допущениях имеют вид[15, 24]:( ) ( ( ) )0ÒpPc T Tρ ⋅ μρ λ ∇ = ∇⋅ − + ∇ ∇     ⋅∇ = ∇⋅ ∇  ∇⋅ = u u I u + uuu(9)где ρ – плотность, кг/м3; u – скорость, м/с;P – давление, Па; I – единичный тензор; μ –вязкость, Па·с; p c – удельная теплоемкость,Дж/кг·К; T – температура, К; λ – коэффициенттеплопроводности, Вт/м·К.Реологическое уравнение нелинейно-вязкой средыимеет вид реологической модели Карро (2).Система уравнений (9) дополнена граничнымиусловиями:622Khvostov A.A. et al. Food Processing: Techniques and Technology, 2021, vol. 51, no. 3, pp. 615–627Ω1: u = −u0n , -n⋅q = ρΔHu⋅n ;Ω2: u = 0 , ( ) î ñ −n⋅q =α (T−T расшT dPpасш 2,5 10 4      α 1104  10 вх T , вх P ос 0shkT kT     – Т); (10) ˜Ω3: [ ] àò ì −PI +K n = P n , -n⋅q = 0где 0 u – скорость входного потока, м/с; n – векторнормали к соответствующей границе расчетнойобласти; q – плотность теплового потока; ΔH –изменение энтальпии, Дж/кг; dPpасш 2,5 10 4      α 1104  10 вх T , вх P ос T0shkT kT     – температура ˜термостатирующей жидкости (воды), омывающейкапилляр, К; α – коэффициент теплоотдачи отконфетной массы к термостатирующей жидкостибез учета тепловой инерционности стенки капилляраВт/м2·К; атм P 1 1 атмвх вхT Pp расшT PH c dT  T dP      pасш 2,5 10 4      α 1104  10 вх T , вх P ос T0shkT kT     – атмосферное давление, Па. ˜Изменение энтальпии ΔH зависит от разноститемператур и давлений:атм P 1 1 атмвх вхT Pp расшT PH c dT  T dP      pасш 2,5 10 4      α 1104  10 вх T , вх P ос T0shkT kT     , (11) ˜где 1 1 атмрасш T dP pасш 2,5 10 4      α 1104  10 вх T , вх P ос T0shkT kT     – температурный коэффициент объемного ˜расширения, К–1; pасш 2,5 10 4      α 1104  10 вх T , вх P ос T0shkT kT     , ˜ pасш 2,5 10 4     α 1104  10 вх T вх P ос T0shkT kT     – температура и дав ˜л ениево входном потоке, К и Па соответственно.Удельная теплоемкость, плотность и коэффициенттеплопроводности конфетной массы, в зависимостиот температуры, определяются уравнениями:p 3782,1 18,708 c = − + T ; (12)ρ = 1969,9 −1,952T ; (13)λ = 0,6205 − 0,0008T (14)Верификация CFD-модели (9) с учетом граничныхусловий (10) проведена в вычислительном пакетеANSYS по следующим исходным данным:геометрические размеры капилляра (диаметр идлина), а также значения объемных расходовконфетной массы через капилляр соответствовалипараметрам, использованных в капиллярнойвискозиметрии; температура конфетной массыи термостатирующей жидкости T = 333,15 К и1 1 атмвх вхT Pp расшT Pc dT  T dP    pасш 2,5 10 4      α 1104  10 вх T , вх P ос T0shkT kT     î ñ 293,15 ˜ T = К соответственно; коэффициенттеплоотдачи от конфетной массы к термостатирующейжидкости α = 5 Вт/м2·К; температурный коэффициентатм объемного расширения P 1 1 атмвх вхT Pp расшT PH c dT  T dP      pасш 2,5 10 4      α 1104  10 вх T , вх P = атм P 1 1 атмвх вхT Pp расшT PH c dT  T dP      pасш 2,5 10 4      α 1104  10 вх T , вх P ос Tshpàñø 2,5 α = ⋅ − К–1.Вычислительный эксперимент проведен всоответствии с представленным выше алгоритмомв интервале изменения параметра регуляризации1 атм P 1 1 атмвх вхT Pp расшT PH c dT  T dP      pасш 2,5 10 4      α 1104  10 вх T , вх P ос T0shkT kT     1⋅10−4 ≤α ≤ 10 . Значения вектора реологичес ˜к ихпараметров α β (при фиксированном α) определялисьиз условия минимизации параметрическойфункции (6). Поскольку целевая функция (6) можетиметь несколько экстремумов, был использованкомбинированный метод случайного поиска,заключающийся в локализации области существованияглобального минимума целевой функции методомсканирования на равномерной сетке переменных ипоследующего уточнения в локализованной областизначения вектора реологических параметров α βметодом статистического градиента [23].На рисунке 3 представлены зависимостирегуляризующей ( ) 1 F β (а) и стабилизирующей( ) 2 F β (b) функций от параметра регуляризации α.Соотношение между регуляризующей истабилизирующей функциями для фиксированныхзначений параметра регуляризации α представленона рисунке 4. Построенный в логарифмическихкоординатах график имеет характерный L-образныйвид с ярко выраженными горизонтальной ивертикальной ветвями. Горизонтальная ветвьL-кривой соответствует решениям, при которыхвеличина регуляризующей невязки ( ) 1 F β болеечувствительна к изменениям параметра регуляризации.Вертикальная ветвь соответствует решениям, вкоторых величина стабилизирующей невязки ( ) 2 F βа bРисунок 3. Зависимость регуляризующей (a) и стабилизирующей (b) функций от параметра регуляризацииFigure 3. Effect of the regularization parameter on the regular izing (a) and stabilizing (b) functions01503004506007500 3 6 9 12 15 1801002003004005006007008000,0E+00 2,0E-07 4,0E-07 6,0E-07 8,0E-07 1,0E-06L/D = 31,25 L/D = 46,87100020003000400050006000700080009000100000,0001 0,001 0,01 0,1 1 1001002003004005000,0001 0,001 0,01 0,1 1 1011010010001000 100000,00,10,20,30,40,0001 0,001 0,01 0,1 1 10Регуляризующая функцияПараметр регуляризации01503004506007500 3 6 9 12 15 1801002003004005006007008000,0E+00 2,0E-07 4,0E-07 6,0E-07 8,0E-07 1,0E-06L/D = 31,25 L/D = 46,87100020003000400050006000700080009000100000,0001 0,001 0,01 0,1 1 1001002003004005000,0001 0,001 0,01 0,1 1 1011010010001000 100000,00,10,20,30,40,0001 0,001 0,01 0,1 1 10Стабилизирующая функцияПараметр регуляризации623Хвостов А. А. [и др.] Техника и технология пищевых производств. 2021. Т. 51. № 3 С. 615–627наиболее чувствительна к изменениям параметрарегуляризации. Таким образом, L-кривая показываетоптимум между параметрической идентификацийпо результатам ротационной вискозиметрии икоррекцией по CFD-модели по данным капиллярнойвискозиметрии.Оптимальное значение параметра регуляризацииα * соответствует вершине сглаженного угла междугоризонтальной и вертикальной ветвями L-кривой,который показывает величину кривизны k (α ) кривой(рис. 5) [17].Решением экстремальной задачи (8) являетсяоптимальное значение параметра регуляризацииα * = 0,027, при котором достигается максимальнаякривизна k (α ) неявно заданной функции. Приоптимальном значении α * = 0,027 коэффициентыреологической модели Карро равны: 0 μ = 16054,5Па·с, μ 1,04 ∞ = Па·с, λ = 155,6 с, n = 0,301 .На рисунке 6 в логарифмических координатахпредставлены экспериментальные данные сдвиговойвискозиметрии помадной массы (графическиепиктограммы) и полная реологическая кривая,построенная по уравнению Карро (2). Видно, чтов диапазоне скоростей сдвига от 0,47 до 669,7 с–1,который соответствует диапазону экспериментальныхисследований, наблюдается удовлетворительноесогласие между экспериментальными значениямивязкости и расчетными, предсказанными по моделиКарро. Максимальная относительная ошибкапредсказания вязкости по уравнению Карро составила14,07 % для ротационной вискозиметрии и 9,17 %для капиллярной вискозиметрии.Рисунок 4. L-криваяFigure 4. L-curve100020003000400050000,0001 0,001 0,01 0,1 1 1001002000,0001 0,001 0,01 0,1 1 1011010010001000 100000,00,10,20,30,40,0001 0,001 0,01 0,1 1 101101001000100001000000,0001 0,001 0,01 0,1 1 10 100 1000 10000 100000 1000000Ротационная вискозиметрия Капиллярная вискозиметрияСтабилизирующая функцияРегуляризующая функцияРисунок 5. Зависимость кривизны L-кривойот параметра регуляризацииFigure 5. Effect of regularization parameter on L-curve100020003000400050000,0001 0,001 0,01 0,1 1 1001002000,0001 0,001 0,01 0,1 1 1011010010001000 100000,00,10,20,30,40,0001 0,001 0,01 0,1 1 101101001000100001000000,0001 0,001 0,01 0,1 1 10 100 1000 10000 100000 1000000Ротационная вискозиметрия Капиллярная вискозиметрияКривизна L-кривойПараметр регуляризацииРисунок 6. Полная реологическая кривая помадной конфетной массыFigure 6. Complete rheological curve for fondant mass01503004506007500 3 6 9 12 15 1801002003004005006007008000,0E+00 2,0E-07 4,0E-07 6,0E-07 8,0E-07 1,0E-06L/D = 31,25 L/D = 46,87100020003000400050006000700080009000100000,0001 0,001 0,01 0,1 1 1001002003004005000,0001 0,001 0,01 0,1 1 1011010010001000 100000,00,10,20,30,40,0001 0,001 0,01 0,1 1 101101001000100001000000,0001 0,001 0,01 0,1 1 10 100 1000 10000 100000 1000000Ротационная вискозиметрия Капиллярная вискозиметрияВязкость, Па·сСкорость сдвига, с–1624Khvostov A.A. et al. Food Processing: Techniques and Technology, 2021, vol. 51, no. 3, pp. 615–627Характер расположения экспериментальныхточек вдоль линейного участка (в логарифмическихкоординатах) полной реологической кривойсвидетельствует об инвариантности методовротационной и капиллярной вискозиметрии [2, 3,10, 11].Полная реологическая кривая μ (γ ) омаднойконфетной массы (рис. 6) имеет вид, характерныйдля структурированных жидкообразных сред,проявляющих аномалию вязкости и обладающихдвумя уровнями постоянной (ньютоновской) вязкости.Для объяснения аномально-вязкого поведенияжидкообразных структур предложено большое числотеорий, использующие различные механические ифизические гипотезы и представления для описанияреологического поведения структурированныхжидкостей. Е. А. Кирсанов и В. Н. Матвеенко всвоей работе приводят описание десяти наиболееизвестных и распространенных концепций, принятыхдля объяснения реологического поведения дисперсныхсистем [2]. Г. В. Виноградов и А. Я. Малкин всеизвестные теории аномалии вязкости условно делят натри большие группы («Реология полимеров», 1977):кинетическая, структурная и гидродинамическаятеории.Кинетическая теория Г. Эйринга и егопоследователей связывает явление аномалиивязкости с влиянием напряжения или скоростисдвига на высоту потенциального барьера, которыйпрепятствует переходу молекулярно-кинетическихединиц из одного равновесного положения в другое. Воснове структурной теории аномалии вязкости лежатпредставления П. А. Ребиндера о влиянии напряженияили скорости сдвига на процессы тиксотропногоразрушения и восстановления структуры системы безучета природы взаимодействия частиц, образующихэту структуру. Согласно гидродинамическим теориямуменьшение вязкости при увеличении скорости илинапряжения сдвига объясняется изменением формымакромолекул и гидродинамическими эффектами,возникающими при их движении в жидкости.Несмотря на то, что в основе этих теорий лежатразличные явления и механизмы объяснения аномалиивязкости, они не являются противоречащими иливзаимоисключающими. Напротив, они дополняютдруг друга и позволяют рассматривать аномально-вязкое поведение структурированных жидких средна различных уровнях.Помадная масса представляет собой высоко-дисперсную концентрированную структури-рованную систему, состоящую из раздробленнойдисперсной фазы, окруженной дисперсионнойсредой. Частицы твердой фазы, взаимодействуядруг с другом через прослойки дисперсионнойсреды, образуют пространственную упорядоченнуюструктуру преимущественно коагуляционного типаотличительными особенностями являются невысокаяпрочность, пластичность, способность к ползучестии тиксотропия [1].Согласно воззрениям П. А. Ребиндера аномально-вязкое стационарное течение жидкообразных средпредставляется как итог двух конкурирующихпроцессов в каждый момент времени: разрушениесвязей в составе пространственной структурнойсетки и тиксотропного восстановления частиэтих разорванных связей в результате теплового(броуновского) движения частиц структурыпри одновременном действии на них внешнегосдвигающего усилия. При отсутствии внешнейсдвигающей силы пространственная структурапомадной массы равнопрочна во всех направлениях.При одномерном сдвиге приложение внешнейсдвигающей силы, не превышающей пределапрочности пространственной структуры, вызываетупорядочение пространственной сетки структуры. Приэтом выделяются связи между частицами продольные,ориентированные по направлению действия силы,и связи поперечные, перпендикулярные к данномунаправлению. Поперечные связи, препятствующиесдвигу слоев жидкообразной среды, разрушаются.Однако при малой скорости сдвига разорванныепоперечные связи успевают восстановиться подвлиянием броуновского движения. Структурированнаясреда в этом случае течет как ньютоновская связкостью 0μ наибольшая вязкость практическинеразрушенной структуры. Дальнейшее увеличениеприкладываемого напряжения сдвига (скоростисдвига) приводит к увеличению числа разрушенныхпоперечных связей. При превышении пределапрочности их число оказывается больше числасвязей, восстановленных под действием броуновскогодвижения. Частицы, освобожденные от поперечныхсвязей, преодолевая силы отталкивания соседнихчастиц, встраиваются в продольные цепочки, упрочняяих. Влияние броуновского движения на частицыцепочки заметно убывает, течение происходит сменьшей вязкостью. При дальнейшем увеличениисдвигающего напряжения (скорости сдвига) плотностьупаковки частиц в продольных цепочках повышаетсяс одновременным увеличением ее прочности. Принекотором критическом значении напряжения сдвигадостигнутая цепочками прочность настолько велика,что в нее не могут встроиться частицы из поперечныхсвязей. Кроме того, число частиц, оставшихся впоперечных связях, мало. Поэтому переход их впродольные цепочки является маловероятным. С этогомомента структурированная жидкость вновь течет какньютоновская, но уже с вязкостью μ∞ – наименьшаявязкость предельно разрушенной структуры. Инымисловами, жидкообразная структура с вязкостью μ∞образована предельно упорядоченными частицами,образующими продольные цепочки. После снятиявнешнего сдвигающего воздействия прослойкимежду частицами продольных цепочек расширяются.При этом продольные цепочки удлиняются иискривляются. Связи между частицами цепочкиослабляются, что облегчает работу броуновскомудвижению. Под действием оставшихся частиц,625Хвостов А. А. [и др.] Техника и технология пищевых производств. 2021. Т. 51. № 3 С. 615–627образующих поперечные связи, а также броуновскогодвижения происходит разрушение уже продольныхсвязей и полное тиксотропное восстановлениеструктуры жидкообразной среды.Концепция П. Я. Ребиндера, хотя и связываетаномалию вязкости с изменением структуры системыпри ее деформировании, не учитывает природувзаимодействия частиц системы и не позволяет ваналитическом виде установить зависимость вязкостиот скорости или напряжения сдвига.Формализация зависимости вязкости структу-рированной жидкообразной системы от условий еетечения возможна с позиций кинетической теорииЭйринга (Н. В. Михайлов и А. М. Лихтгейм (1955 г.)Г. М. Бартенев (1955 г.), Е. Е. Бибик (1981 г.)).Согласно Эйрингу вязкость неньютоновской жидкостиубывает с ростом напряжения сдвига в соответствиис законом0 4 shkTτμτων=   (15)где τ – напряжение сдвига; 0 ν – частота перескоковмолекулярно-кинетической единицы жидкостив единицу времени при отсутствии внешнегосдвигающего усилия; ω – эффективный объеммолекулярно-кинетической единицы; k – постояннаяБольцмана; T – температура.При напряжении сдвига, стремящегося кбесконечно малому значению, справедливаэквивалентность pасш 2,5 10 4     α 1104  10 вх T , вх P ос T0shkT kT     ˜ 10 вх T , вх P ос T0shkT kT     ˜ . В этом случае (15)запишется как00 4μ μ kTν ω= = (16)Это подтверждает тот факт, что вязкостьструктурированной жидкообразной среды прискорости сдвига, стремящейся к нулю, не зависит отнапряжения сдвига и принимает свое максимальноеконечное значение (наибольшая вязкость практическинеразрушенной структуры).ВыводыВ работе показано, что степенное реологическоеуравнение Оствальда-де Виля является физическинекорректным при предсказаниях значений вязкости впредельных случаях. При решении задач ламинарныхпространственных течений нелинейно-вязких средв изотермической и неизотермической постановкахиспользование степенного уравнения Оствальда-деВиля сопряжено с вычислительными проблемами,которые приводят к существенным погрешностям.Для устранения вычислительных проблем предложеноиспользовать реологическое уравнение Карро,учитывающее предельные ньютоновские состояниянелинейно-вязкой среды.Анализ показал, что существующие методы иприборное оформление сдвиговой вискозиметрииненьютоновских сред не позволяют оценитьпредельные значения вязкости. В этой связи дляпараметрической идентификации реологическоймодели Карро предложено использоватьмодифицированный алгоритм регуляризацииТихонова. Согласно алгоритму регуляризованноерешение (в виде вектора реологических параметров)получается в результате достижения баланса междупредпочтительным выбором результатов ротационнойвискозиметрии и коррекцией результатов поCFD-модели течения структурированнойжидкообразной среды в измерительной системевискозиметра.Для реализации алгоритма регуляризацииформализована CFD-модель течения нелинейно-вязкой среды с уравнением Карро в цилиндрическомкапилляре. CFD-модель составлена на основеуравнений сохранения массы, энергии и импульсаи дополнена соотношениями, определяющимиграничные условия и температурные зависимоститеплофизических свойств среды.На примере помадной массы кондитерскогопроизводства, проявляющей неньютоновскиесвойства, показана процедура идентификациипараметров реологического уравнения Карро. Полу-ченное уравнение предсказывает вязкость помадноймассы с ошибкой, не превышающей 14,07 %,и позволяет устранить вычислительные проблемы,характерные для реологической модели Оствальда-деВиля, возникающие при решении задач ламинарныхпространственных течений нелинейно-вязких средв изотермической и неизотермической постановках.Критерии авторстваАвторы в равной степени участвовали в подготовкеи написании статьи.Конфликт интересовАвторы заявляют об отсутствии конфликтаинтересов.ContributionAll the authors bear equal responsibility for thecontent of the article.Conflict of interestThe authors declare that there is no conflict of interestregarding the publication of this article.