Программные продукты и системы

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

Математическая модель задачи Постановка плоской задачи теории упругости. Для плоской кусочно-однородной области W, на внешней границе которой заданы некоторые нагрузки или перемещения, рассматривается статическая задача теории упругости, заключающаяся в определении в каждой из составляющих область W зон W(1), W(2), …, W(n) вектора перемещений ui, тензора деформаций eij и тензора напряжений sij, которые в рассматриваемой зоне удовлетворяют системе уравнений

sij,i=0,                                                                      (1)

,                                                 (2)

                                     (3)

и заданным на внешней границе S=S1ÈS2 области W граничным условиям

, .                        (4)

Здесь по повторяющемуся индексу производится суммирование от 1 до 2; ui,j=¶ui/¶xj; m – модуль упругости при сдвиге; n – коэффициент Пуассона; dij – символ Кронекера; ni – вектор нормали к границе зоны. Значения модуля упругости и коэффициента Пуассона для разных зон различные. На внутренних границах, являющихся границами между зонами, принимается условие непрерывности перемещений и напряжений.

Метод граничных элементов для неоднородной области. Согласно данному методу [1], внутри каждой зоны W(i) выполняется тождество Сомильяны:

, (5)

где S(i) – граница зоны W(i); x – внутренняя точка W(i); функции влияния  и  для двухмерной задачи теории упругости имеют следующий вид:

,                          (6)

. (7)

Здесь ; c2=(3–4n); ; c4=(1–2n); r=r(x, x) – расстояние между точками x и x; .

Выражение (5) с учетом уравнений (2) и (3) дает непрерывное решение задачи теории упругости в области W(i), если найдены компоненты векторов поверхностного напряжения fi(x) и перемещения ui(x) на границе S(i). Частично эти функции определяются из граничных условий (4).

Предельный переход x®x0, x0ÎS(i), приводит к граничному интегральному уравнению для каждой зоны:

   (8)

Здесь ,

,

,

,                         (9)

где w - величина угла внутри деформируемой области, образованного в граничной точке x0. Когда точка x0 лежит на гладкой поверхности, .

Граничное интегральное уравнение позволяет определить недостающие значения поверхностных перемещений и напряжений. В однородной области (одной зоне) решение уравнений (8) сводится к разбиению границы S на конечное число элементов, построению интерполяционных функций для неизвестных значений ui(x), fi(x) и решению получающейся из (8) системы алгебраических уравнений относительно параметров интерполяционных функций. В простейшем случае можно считать компоненты векторов поверхностного напряжения и перемещения постоянными на каждом прямолинейном элементе. Тогда уравнение (8) составляется для единственного узла на каждом элементе, в результате чего получается следующая система уравнений:

,

i=1, 2; p=1, …, N,                                                        (10)

где N – число элементов; e(k) – граничные элементы; ,  – значения компонентов векторов перемещения и поверхностного напряжения на элементе e(k). В неоднородной области, где граничные условия (4) заданы только на внешней границе всей области W, систему уравнений, объединяющую системы вида (10), построенные для каждой из зон, необходимо дополнить условиями непрерывности перемещений и напряжений на общих границах зон. Решение полученной системы линейных алгебраических уравнений даст значения поверхностных напряжений и перемещений на всех границах S(i), что позволит определить из соотношений (5), (2), (3) НДС внутри зон W(i).

Для реализации описанного алгоритма необходимо вычислять входящие в уравнения (5) и (10) криволинейные интегралы по граничным элементам от компонент функций влияния. Стандартный подход к решению предполагает численное вычисление этих интегралов в каждой конкретной задаче. Более универсальным подходом является получение аналитических формул для точного вычисления всех необходимых при решении задачи интегралов. Получение таких формул позволяет включить их в алгоритм, тем самым повысив его универсальность. Кроме того, аналитическое интегрирование по сравнению с численным имеет очевидные преимущества по точности и скорости вычислений. В данной работе использованы аналитические формулы, полученные в [2]. Их применение ускоряет счет в 5–10 раз, а выигрыш в точности соответствует погрешности интерполяции, применяемой при численном интегрировании. Сопоставление полученных решений различных задач деформирования для плоских областей с известными аналитическими решениями, а также c решениями, полученными в ANSYS, показало высокую эффективность разработанных алгоритмов [2].

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

Пакет программ для анализа НДС плоской неоднородной области

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

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

-      модуль обмена данными с внешними файлами;

-      графический редактор – построение геометрической модели плоского сечения неоднородной конструкции;

-      модуль численного решения плоской задачи упругости – вычисление НДС плоского сечения конструкции;

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

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

На рисунке 2 представлена общая блок-схема пакета.

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

Графический редактор предназначен для построения новой и редактирования существующей геометрической модели плоской неоднородной области, определения внешней статической нагрузки (точек приложения, направлений и значений внешних сил), указания закрепленной части внешней границы области, формирования дополнительных структур данных (вариабельные части контактных границ образующих область зон, анализируемая зона области). Инструментальные средства этого модуля позволяют пользователю в различных режимах строить замкнутые ломаные линии, определяющие границы соответствующих зон (рис. 1). В частности, при построении новой ломаной автоматически активизируются режимы захвата вершин и ребер уже существующих ломаных, позволяющие обеспечить слипание ломаных на контактирующих участках (отсутствие «дырок»). Для редактирования существующей ломаной реализованы операции удаления ее вершин, добавления вершин на ребра ломаной, перемещения вершин (рис. 1).

Закрепленная граница области задается набором ребер ломаных, принадлежащих ее внешней границе (на рис. 1 отображена жирной черной линией).

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

В рамках графического редактора реализована подсистема управления параметрами и режимами визуализации: масштаб визуализации, маркеры (кружки) вершин и их размер, заливка зон и маркеров вершин ломаных, направление обхода по ломаным, размеры рабочего поля, величина отступов на рабочем поле и т.д.

Модуль вычисления НДС плоской неоднородной области предназначен для численного решения плоской задачи упругости (1)−(4) на основе текущей геометрической модели плоского сечения конструкции и исходных данных – заданной внешней нагрузки, закрепленной границы, значений упругих параметров зон (коэффициентов Пуассона и модулей Юнга). В этом модуле реализован описанный выше алгоритм решения задачи (1)−(4) на основе метода граничных элементов. На рисунке 3 показан вывод результатов вычисления НДС плоской неоднородной области.

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

Модуль оптимизации формы границы служит для численного решения одного класса задач оптимизации формы области.

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

Анализ НДС восстановленного при помощи цельнокерамического микропротеза депульпированного зуба

Интересным примером задачи анализа свойств неоднородной упругой конструкции может служить возникающая в стоматологии задача, которая заключается в определении оптимальной формы цельнокерамической реставрации – микропротеза, используемого для восстановления депульпированных зубов [3]. Как показывает клиническая практика, при восстановлении зуба с помощью такого микропротеза геометрические свойства его формы являются одним из важных факторов, определяющих успех протезирования. Оказывается, за счет оптимального выбора формы микропротеза во многих случаях значительно увеличивается жизненный цикл реставрации. Поскольку возможности варьирования формы жевательной (внешней) поверхности микропротеза очень ограниченные, изменения формы этой части его поверхности могут быть весьма незначительными. В этих условиях увеличение жизненного цикла реставрации может достигаться путем формирования рациональной формы лишь внутренней поверхности микропротеза – контактной границы «микропротез–ткани зуба». Таким образом, возникает целесообразность (а в ряде клинических случаев иПодпись: необходимость) в решении задачи оптимизации формы упругонапряженного тела – микропротеза в депульпированном зубе, подверженном воздействию заданной внешней нагрузки. При этом на допустимую форму микропротеза накладываются как габаритные ограничения – форма его внешней границы фиксирована, так и изопериметрические – объем материала микропротеза не должен превышать заданной величины, что при протезировании позволяет в максимальной степени сохранить здоровые ткани зуба. В такой задаче вид оптимизируемого показателя, необходимые и достаточные условия оптимальности контактной границы «микропротез–ткани зуба», а также численные методы ее построения [4] основаны на свойствах распределения возникающих на контактной границе напряжений. Перечисленные обстоятельства обусловливают необходимость достаточно точного вычисления значений компонент тензоров напряжений и деформаций на контактных границах, образующих конструкцию областей.

Геометрическая модель конструкции. Рассматривается неоднородная упругая область, подверженная воздействию заданной внешней статической нагрузки. Эта область представляет собой объединение конечного числа упругоконтактирующих друг с другом однородных подобластей (зон), обладающих различными механическими свойствами. Частным примером такой неоднородной упругой конструкции является реставрированный при помощи цельнокерамического микропротеза депульпированный зуб, состоящий из трех материалов – дентина, зубной эмали и керамики (рис. 4). Все материалы жестко связаны друг          с другом – при упругой деформации конструкции смещения на контактных границах зон, соответ-ствующих различным материалам, ничтожно    малы.

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

В данной работе рассматривается упрощенная постановка задачи упругости, основанная на общепризнанном экспертами в области сопротивления материалов методе плоских сечений трехмерной конструкции [5]. Для построения двухмерной геометрической модели плоского сечения неоднородной конструкции, например зуба, может использоваться рентгенограмма (рис. 4, 5).

Геометрическая модель плоского сечения конструкции (рис. 5) описывается плоской неоднородной областью, состоящей из конечного числа контактирующих зон. Каждая зона соответствует некоторому материалу (рис. 4) и задается аппроксимирующей ее границу замкнутой ломаной, определяемой последовательностью пар координат ее вершин в декартовой системе координат.

Внешняя статическая нагрузка задается набором внешних сил, приложенных к участкам внешней границы области. Как уже отмечалось, каждая из сил характеризуется направлением (вектором), длиной площадки приложения и величиной силы (рис. 5). Центр площадки приложения силы является серединной точкой ребра ломаной, принадлежащего внешней границе области. Часть внешней границы области может быть жестко закреплена – в точках этой части границы перемещения полагаются равными нулю (на рис. 5 отображена жирной черной линией).

Анализ НДС плоского сечения восстановленного зуба. Результаты численного моделирования. Описанный выше алгоритм численного решения плоской задачи упругости используется для анализа НДС плоского сечения восстановленного с помощью цельнокерамического микропротеза зуба. Эта задача решается при следующих предположениях:

1)    сечение зуба жестко закреплено на его нижней части границы (смещения равны нулю в точках этой границы) – игнорируется отдача корней зуба (рис. 5);

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

Указанные допущения адекватно отражают лабораторные условия проведения эксперимента.

Геометрическая модель плоского сечения зуба (рис. 5) описывается плоской неоднородной областью, состоящей в общем случае из трех типов (по типам материалов) контактирующих зон. На рисунке 5 зона дентина отображается светло-серым цветом, зона зубной эмали – белым цветом, зона керамики – темно-серым цветом, закрепленная граница указана жирной черной линией. Замкнутые ломаные, аппроксимирующие границы соответствующих зон, строятся по изображению рентгенограммы (рис. 4, 5). При этом на контактных участках вершины ломаных, соответствующих различным контактирующим зонам, совпада-      ют (рис. 5).

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

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

Целью решения задачи теории упругости являлась оценка распределения напряжений на границах зон. Особый интерес представляет контактная граница «микропротез–ткани зуба» (на рис. 5 изображена пунктирной линией белого цвета), так как распределение напряжений на этой границе может использоваться для оценки качества формы плоского сечения микропротеза и формирования прогноза длительности жизненного цикла реставрации в целом.

Свойства распределения напряжений на контактной границе «микропротез–ткани зуба» положены в основу численного метода решения задачи оптимизации формы этой границы. Как уже отмечалось, подобная задача заключается в определении такой ее геометрической формы, при которой реализовывалось бы оптимальное распределение напряжений на этой границе. При этом на допустимую форму сечения микропротеза накладывается изопериметрическое ограничение – площадь сечения должна быть равна заданной величине. Оптимальной формой контактной границы «микропротез–ткани зуба» считается такая, которая при заданной внешней нагрузке на область обеспечивает равномерное распределение напряжений на этой контактной границе.

Содержательно выполнение настоящего условия приводит к тому, что при заданной внешней нагрузке материал восстановленного зуба работает одинаково во всех точках вблизи контактной границы «микропротез–ткани зуба». Это обстоятельство позволяет, в свою очередь, избежать возникновения пиковых нагрузок на контактной границе, приводящих к уменьшению жизненного цикла реставрации. Поэтому в рассматриваемой задаче оптимизируемой характеристикой является величина (величины), отражающая степень равномерности распределения напряжений на варьируемой границе сечения зуба.

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

Пакет прикладных программ прошел опытную апробацию на кафедре ортопедической стоматологии Уральской государственной медицинской академии (г. Екатеринбург). Реализованные в пакете функциональные возможности использовались для прогнозирования результатов реставрации депульпированных зубов реальных пациентов цельнокерамическими микропротезами. С помощью этого пакета проводился численный анализ НДС плоских сечений восстанавливаемых зубов с целью выработки рекомендаций по формированию рациональной формы контактных границ «микропротез–ткани зуба». В ряде случаев полученные результаты такого численного анализа (рис. 3) учитывались в клинической практике.

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

Литература

1.     Бенерджи П., Баттерфилд Р. Методы граничных элементов в прикладных науках. М.: Мир, 1984. 494 с.

2.     Федотов В.П., Спевак Л.Ф. Модифицированный метод граничных элементов в задачах механики, теплопроводности и диффузии. Екатеринбург: УрО РАН, 2009. 164 с.

3.     Кандоба И.Н., Коледа П.А., Жолудев С.Е. Оптимизация формы накладки (pinlay). Моделирование упругонапряженного состояния зуба // Проблемы стоматологии. 2007. № 5. С. 19–23.

4.     Кандоба И.Н. Метод решения одного класса задач оптимизации формы неоднородной области. Забабахинские научные чтения: сб. тез. докл. X Междунар. конф. (15–19 марта 2010 г.). Снежинск: Изд-во РФЯЦ-ВНИИТФ, 2010. С. 282–283.

5.     Беляев Н.М. Сопротивление материалов. М.: Физматлит, 1962. 856 с.

width=