Online tool for finite element analysis of postmortem convective heat transfer of the head

Full Text

Обоснование

Одним из наиболее перспективных и активно развивающихся современных направлений в области термометрической диагностики давности наступления смерти (ДНС) является метод конечно-элементного анализа посмертного теплообмена [1-3]. Данный подход позволяет преодолеть ограничения, присущие феноменологическим уравнениям охлаждения трупа, поскольку непосредственно основывается на физике теплопередачи и принципиально может быть применен к любым условиям охлаждения [4-7].

Метод конечных элементов (МКЭ), используемый при оценке посмертного теплообмена, является валидированным инструментом решения в сложных геометрических объектах и краевых условиях дифференциальных уравнений с частными производными 2-го порядка, например, таких, как уравнение теплопроводности, которое в пространственной области Ω размерности n с системой координат  имеет вид

                                            ,                                         (1)

где c – удельная теплоемкость, ρ – плотность, λ – теплопроводность, r – координаты точки,  - искомая функция, задающая температуру в точке пространства Ω в момент времени τ на промежутке [0,∞);  - функция плотности источников тепла [6].

После развития компьютерных технологий МКЭ получил широкое распространение при решении различных термодинамических задач в технике, а также и медицине. В настоящее время на рынке присутствует значительное количество различных коммерческих пакетов программного обеспечения МКЭ. Высокая стоимость этих продуктов, даже отечественного производства, существенно ограничивает возможности их использования для диагностики ДНС. Бесплатные же версии программных пакетов, выпускаемые в образовательных целях, характеризуются наличием ограничений на количество узлов в генерируемых сетках конечно-элементных моделей (КЭМ) или вовсе допускают постановку только двумерных задач.

Кроме того, использование коммерческих решателей предполагает наличие у пользователей компетенций по самостоятельному построению геометрии КЭМ, заданию ее теплофизических свойств и краевых условий. Это препятствует широкому внедрению МКЭ для диагностики ДНС, сводя его к единичным эпизодам применения в научных целях, к тому же реализуемым коллективами исследователей, помимо судебных медиков включающими в свой состав математиков и программистов [2, 4-7].

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

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

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

С учетом названных ограничений одной из наиболее перспективных для конечно-элементного моделирования областей тела представляется голова, геометрия мозгового отдела которой очень близка к такому телу, как многослойный шар с равномерным распределением слоев однородной структуры [8]. Также поверхность шара имеет с касательной плоскостью только одну общую точку, что позволяет пренебречь потерями тепла путем теплопроводности при контакте головы с опорной поверхностью и не учитывать в КЭМ теплофизические свойства последней.

ЦЕЛЬ

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

МЕТОДЫ

Дизайн исследования

Дизайн исследования представляет собой разработку КЭМ посмертного конвективного теплообмена головы с последующим построением на ее основе онлайн-приложения, предназначенного для определения ДНС путем краниоэнцефальной термометрии трупа.

Геометрия КЭМ

В качестве геометрической модели головы взрослого человека и, прежде всего, ее мозгового отдела использовали многослойный шар радиусом 98 мм, состоящий из пяти равномерно распределенных однородных слоев: кожно-апоневротического лоскута (5 мм), костей свода черепа (5 мм), ликвора субарахноидального пространства (2 мм), а также поверхностного (17 мм) и глубокого (69 мм) слоев головного мозга (86 мм). Два подслоя в блоке головного мозга выделяли в целях моделирования его начального температурного поля у человека, которое в поверхностном слое является неоднородным с монотонным понижением температуры [9, 10]. Непосредственно расчетной моделью являлось сечение шара диаметральной плоскостью с нулевой аппликатой всех его точек (рис. 1). Сгенерированная сетка КЭМ включала 1311 узлов и 9277 конечных элементов, в т.ч. 2552 треугольников и 6725 тетраэдров (см. рис. 1).

Рисунок 1. КЭМ мозгового отдела головы, заданная в онлайн-приложении.

Figure 1. The finite element model of the brain region of the head, set in an online application.

Теплофизические свойства КЭМ

В онлайн-приложении предусматривали возможность произвольного задания теплофизических свойств выделенных в КЭМ анатомических слоев головы. Поскольку на практике точное измерение удельной теплоемкости, теплопроводности и плотности органов и тканей трупа затруднительно, по умолчанию задавали усредненные теплофизические свойства актуальных тканей (см. таблицу), указанные в специальной литературе [1, 6-10].

Таблица. Теплофизические свойства, по умолчанию заданные в КЭМ

Table. Thermophysical properties, set by default in the finite element model

Начальное температурное поле

Для нахождения начального температурного поля КЭМ решали уравнение Лапласа

c краевыми условиями Дирихле, соответствовавшими постоянной температуре на поверхностях выделенных анатомических слоев:

,

где Δt – оператор Лапласа, в декартовой системе координат трехмерного пространства с осями ox, oy и oz определяемый как

,

а символ ∂Ω обозначает границу области Ω.

По умолчанию в качестве условий Дирихле задавали усредненные прижизненные температуры на внешних границах слоев (см. таблицу). Затем, используя слабую формулировку уравнения Лапласа, выводили матричное уравнение [11, 12], из которого находили начальную температуру в каждой точке КЭМ (рис. 2).

Посмертное температурное поле

На следующем этапе находили посмертное температурное поле расчетной области. Сумму выделяемой внутренними источниками тепловой энергии и теплопотерь за счет радиационного теплообмена ввиду их небольшой величины и примерно одинаковых количеств принимали равной нулю [8]. Также предполагали, что в силу однородности выделенных в КЭМ анатомических слоев их теплофизические свойства являются константами. С учетом названных допущений уравнение теплопроводности (1) принимает вид

,

где а – коэффициент температуропроводности, м²/с, равный

,

где c – удельная теплоемкость, Дж/(К∙кг), ρ – плотность, кг/м³; λ – коэффициент теплопроводности, Вт/(м∙К).

Рисунок 2. Начальное температурное поле диаметрального сечения КЭМ.

Figure 2. The initial temperature field of the diametrical section of the finite element model

В качестве начального условия поставленной задачи нестационарного теплопереноса принимали найденное на первом этапе стационарное начальное температурное поле

.

Для моделирования охлаждения в воздушной среде на внешнем ребре модели задавали конвективный теплообмен с воздухом, протекающий по закону Ньютона-Рихмана:

,

где n – нормаль к поверхности охлаждаемой части тела, м; s – индекс, обозначающий поверхность физического тела; α – коэффициент теплоотдачи, Вт/(м²·К); t_a – температура внешней среды, К. Коэффициент теплоотдачи в расчетах задавали постоянным и одинаковым для всей поверхности КЭМ. Актуальный расчетный период принимали равным 24,5 ч с шагом интегрирования в 600 с.

Валидация КЭМ

В ходе валидации онлайн-приложения проверяли адекватность выполняемых им расчетов МКЭ и адекватность геометрии КЭМ и заданных теплофизических свойств.

Адекватность МКЭ тестировали путем оценки соответствия посмертных температурных полей КЭМ при различных условиях конвективной теплоотдачи и внешнего температурного режима с таковыми, полученными с помощью бесплатной версии приложения ELCUT 6.6 (доступного по URL адресу: https://elcut.ru/downloads_r.htm) для КЭМ с аналогичной геометрией, теплофизическими свойствами и краевыми условиями.

Адекватность геометрии КЭМ и заданных теплофизических свойств подтверждали путем оценки соответствия посмертных температурных кривых для центра и поверхности КЭМ, полученных с помощью разработанного онлайн-приложения, таковым, вычисленным по уравнениям с общепризнанной валидностью, основанным на законах охлаждения Ньютона-Рихмана и Marshall-Hoare [8, 13].

В качестве валидного описания посмертной динамики краниоэнцефальной температуры использовали уравнение Henssge:

                                                     ,                                                  (2)

где t –краниоэнцефальная температура, °С; t_а – внешняя температура, °С; τ – ДНС, ч [8, 14].

В качестве эталона динамики охлаждения поверхности головы применяли уравнение, основанное на законе охлаждения Ньютона-Рихмана:

                                                            ,                                                           (3)

где t –поверхностная температура, °С; t₀ – начальная поверхностная температура, °С; k – константа охлаждения, ч^-1 [8, 14].

В уравнении (3) начальную температуру поверхности головы задавали равной 31 °С, а константу k определяли по формуле:

,

где t₁ и t₂ – температура поверхности головы в 23 и 24 ч посмертного периода соответственно, ºС; а Δτ – промежуток времени между указанными измерениями, равный 1 ч [8].

В ходе валидации КЭМ находили и оценивали ее посмертное температурное поле при внешней температуре от 10 до 22 ºС и коэффициенте теплоотдачи от 5 до 25 Вт/(м²·К).

Язык программирования

Код онлайн-приложения, предназначенного для операционных систем Windows и Linux, составляли на языке программирования Python 3. При разработке пользовательского интерфейса также использовали библиотеку React из языка JavaScript. В качестве асинхронной очереди задач на расчет и брокера сообщений применяли инструменты Celery и Redis.

РЕЗУЛЬТАТЫ

На основе разработанной КЭМ было создано приложение «Simple Finite Element Model of Postmortem Convective Heat Transfer of the Head» (свидетельство о государственной регистрации программы для ЭВМ № 2024668362) объемом 1,779 Гб, предназначенное для конечно-элементного онлайн-анализа посмертного конвективного теплообмена головы и доступное по URL-адресу: https://forenscalc.ru/simple-fem.

Для работы с онлайн приложением пользователю необходимо указать координаты диагностической точки из области мозгового отдела головы, температуру внешней среды в посмертном периоде, коэффициент теплоотдачи, толщину и теплофизические свойства анатомических слоев головы (начальную температуру, плотность, теплопроводность и удельную теплоемкость). Результатом работы приложения является визуализация внешнего температурного тренда и кривых охлаждения в диагностической точке и на поверхности головы в первые сутки посмертного периода с генерацией числовых значений указанных температур в виде файла с расширением .csv, доступного в личном кабинете пользователя. Построение графиков приложение выполняет в режиме real-time. Помимо числовых данных в личном кабинете пользователя дополнительно формируются файлы с расширениями .msh и .png, содержащие графическую и табличную информацию о геометрии и свойствах сетки КЭМ, а также начальную термограмму вдоль контура интегрирования от центра КЭМ до диагностической точки.

Сравнение кривых охлаждения, полученных в одинаковых условиях с помощью разработанного онлайн-приложения и программы ELCUT 6.6, установило их практически полное соответствие. Незначительные несовпадения были обусловлены небольшим количеством узлов и элементов в сетке КЭМ бесплатной версии приложения ELCUT 6.6. и, следовательно, меньшей точностью ее вычислений. Дальнейшая отладка КЭМ в различных условиях конвективного теплообмена с отсутствием внешней радиации и контакта с другими физическими телами показала максимальную сходимость ее вычислительных результатов с таковыми уравнений (2) и (3) при задании коэффициента теплоотдачи, равного 6 Вт/(м²·К).

Приложение помимо ввода константных значений внешней температуры допускает задание и любых типов ее изменений, описываемых функциями из библиотеки Python 3. Для этого в соответствующее текстовое поле достаточно ввести составленную в соответствии с синтаксисом Python-библиотеки NumPy функцию внешней температуры от ДНС, измеряемой в секундах.

При отсутствии у судебно-медицинского эксперта информации о биометрических и теплофизических свойствах тканей головы конкретного трупа или их начального температурного поля, а также характеристиках конвективного теплообмена следует проводить расчеты, используя заданные по умолчанию значения, в т.ч. и коэффициента теплоотдачи, равного 6 Вт/(м²·К). В этом случае при постоянной внешней температуре результаты МКЭ для точки начала координат в КЭМ будут соответствовать таковым уравнения C. Henssge.

Пример

Внешняя температура в посмертном периоде изменялась по синусоидальной кривой, описываемой уравнением

,

где t – ДНС, с (рис. 3). Необходимо вычислить динамику краниоэнцефальной температуры трупа при усредненных биометрических и теплофизических параметрах тканей головы и конвективном теплообмене с коэффициентом теплоотдачи 6 Вт/(м²·К).

Для решения задачи в соответствующих текстовых полях онлайн-приложения «Simple Finite Element Model of Postmortem Convective Heat Transfer of the Head» необходимо указать декартовы координаты диагностической точки, равные нулю. Графические результаты расчетов представлены на рисунке 3.

Рисунок 3. Графики охлаждения центра и поверхности КЭМ при синусоидальных изменениях внешней температуры. Условные обозначения: 1 – температура в точке с нулевой радиальной координатой; 2 – температура в точках поверхности КЭМ; 3 – внешняя температура.

Figure 3. Graphs of cooling of the center and surface of the finite element model with sinusoidal changes in external temperature. Symbols: 1 – temperature at a point with zero radial coordinate; 2 – temperature at points on the surface of the finite element model; 3 – ambient temperature.

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

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

ОБСУЖДЕНИЕ

Поскольку возможности феноменологических уравнений охлаждения трупа в плане термометрической диагностики ДНС практически исчерпаны, основным методом решения данной научно-практической задачи остается конечно-элементный анализ посмертного теплообмена мертвого тела в целом или его отдельных частей. Однако существенным препятствием на пути внедрения МКЭ остаются его сложность, необходимость приобретения коммерческих пакетов программного обеспечения и наличия у пользователя умений по их использованию, наиболее важными из которых следует считать построение геометрии и сетки КЭМ, задание ее теплофизических свойств и краевых условий. В ходе проведенного исследования впервые разработано открытое судебно-медицинское онлайн-приложение конечно-элементного анализа посмертного теплообмена, позволяющее в режиме real-time рассчитать кривую охлаждения для любой точки мозговой части головы трупа и не требующее от пользователя компетенций по заданию геометрии, теплофизических свойств и краевых условий КЭМ.

Геометрия КЭМ созданного онлайн-приложения является масштабируемой и аппроксимирована многослойным шаром, поскольку при условии адекватной фиксации общих объемов тканей геометрическое разрешение модели практически не влияет на точность определения ДНС [5]. Правомерность использованной аппроксимации уже была подтверждена ранее на двумерной КЭМ с грубой сеткой всего из 254 конечных элементов однотипной треугольной формы [8]. Названное упрощение геометрии, выполненное с повышением дискретизации расчетной области до 9277 конечных элементов, позволило значительно уменьшить вычислительные затраты и автоматизировать процесс создания КЭМ, а также сделало возможным выполнение расчетов в онлайн-режиме, при этом повысив их точность по сравнению с двумерными КЭМ.

При наличии соответствующих данных пользователь может задавать в расчетах толщину каждого из 5 актуальных анатомических слоев: кожно-апоневротического лоскута, свода черепа, лептоменинкса с ликвором, поверхностного и глубокого блоков мозга, а также их основные теплофизические свойства и начальную температуру их внешних границ. Это позволяет учитывать реальную геометрию и размеры расчетной области и составляющих ее тканевых структур. Кроме того, доступно задание величины конвективного теплообмена, значение коэффициента теплоотдачи которого зависит от типа и скорости перемещения охлаждающей труп среды [15-18].

Важным преимуществом разработанного онлайн-приложения, как и всех аналогичных программных пакетов, основанных на МКЭ, является возможность расчета кривой охлаждения для любой точки мозгового отдела головы. Интерфейс онлайн-приложения позволяет указать декартовы координаты сечения головы диаметральной плоскостью. Однако ввиду осевой симметрии КЭМ достаточно, приняв одну из декартовых координат равной нулю, вместо второй задать радиальную координату точки, соответствующую расстоянию от центра КЭМ, расположенного в толще межталамического сращения головного мозга, до конца раневого канала от игольчатого термощупа. Возможность учета в расчетах координат диагностической точки позволяет устранить такой присущий феноменологическим уравнениям тип погрешностей определения ДНС, как диагностические промахи, связанные с отклонениями конца термощупа от точки с наивысшей температурой при проведении краниоэнцефальной термометрии [19].

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

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

В целом перечисленные ограничения могут быть устранены, но ценой необходимости возложения на пользователя обязанности самостоятельного задания геометрии и краевых условий КЭМ. По этой причине выбранный разработчиками приложения диапазон учитываемых при вычислениях условий охлаждения и геометрических характеристик расчетной области является компромиссом между возможностями моделирования последних и простотой интерфейса приложения.

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

ЗАКЛЮЧЕНИЕ

Одним из наиболее перспективных современных направлений в области термометрической диагностики ДНС является метод конечно-элементного анализа посмертного теплообмена, позволяющий преодолеть ограничения, присущие феноменологическим уравнениям охлаждения трупа. Однако доступные пакеты программного обеспечения МКЭ характеризуются высокой стоимостью, а также предполагают наличие у пользователя навыков самостоятельного задания условий решаемых задач и анализа результатов выполненных вычислений. В ходе проведенного исследования разработано открытое онлайн-приложение конечно-элементного анализа посмертного конвективного теплообмена головы человека, предназначенное для определения ДНС путем краниоэнцефальной термометрии трупа. Для достижения возможности выполнения вычислений в онлайн-режиме геометрия головы при расчетах аппроксимирована многослойным шаром. В отличие от феноменологических уравнений разработанное онлайн-приложение при определении ДНС позволяет учитывать особенности начального температурного поля, размеры и теплофизические свойства основных анатомических слоев головы, интенсивность конвективного теплообмена, координаты диагностической точки и изменения внешней температуры в процессе охлаждения трупа. Результатом работы онлайн-приложения являются визуализация кривых охлаждения в диагностической точке и на поверхности головы в первые сутки посмертного периода с генерацией соответствующих числовых значений, а также вывод информации о геометрии и свойствах сетки КЭМ и распределении температуры вдоль контура интегрирования от центра КЭМ до диагностической точки. История и результаты расчетов сохраняются в личном кабинете пользователя. В отличие от других программных пакетов разработанный решатель не предъявляет высоких системных требований к компьютеру, поскольку реализован в виде веб-приложения с выполнением всех вычислений на сервере, и не требует наличия у пользователя специальной математико-инженерной подготовки. Изложенное позволяет использовать МКЭ в формате разработанного онлайн-инструмента в судебно-медицинской экспертной практике при установлении ДНС путем краниоэнцефальной термометрии трупа.

ДОПОЛНИТЕЛЬНАЯ ИНФОРМАЦИЯ

Источник финансирования. Авторы заявляют об отсутствии внешнего финансирования при проведении исследования.

Конфликт интересов. Авторы декларируют отсутствие явных и потенциальных конфликтов интересов, связанных с публикацией настоящей статьи.

Участие авторов. Авторы подтверждают соответствие своего авторства международным критериям ICMJE (все авторы внесли существенный вклад в разработку концепции, проведение исследования и подготовку статьи, прочли и одобрили финальную версию перед публикацией). Наибольший вклад распределён следующим образом: В.Г. Недугов – разработка вычислительного алгоритма МКЭ, написание исходного кода онлайн-приложения, научное редактирование рукописи, рассмотрение и одобрение окончательного варианта рукописи; Г.В. Недугов - концепция и дизайн исследования, разработка геометрии конечно-элементной модели, сбор данных по теплофизическим свойствам тканей головы, тестирование и валидация онлайн приложения, написание текста рукописи.

About the authors

German Nedugov

Author for correspondence.
Email: nedugovh@mail.ru
ORCID iD: 0000-0002-7380-3766
SPIN-code: 3828-8091
Scopus Author ID: 25947646500

PhD, docent of the Department of forensic medicine of Samara State Medical University

Vladimir Nedugiv

Email: nedugovvg@gmail.com
ORCID iD: 0009-0007-7542-7235
SPIN-code: 2407-7937
Scopus Author ID: 58092580600
Россия, 443086, Russia, Samara, Moskovskoye shosse, 34

References

Supplementary files