Численное дифференцирование в задачах исследования двумерных скалярных полей океанаграфических характеристик

1. Введение

В последние десятилетия в математике вопросам численного дифференци¬рования уде¬ляет¬ся особенно пристальное внимание. Это объясняется тем обстоятель¬ством, что большинство актуальных инженерно-технических задач, подлежащих ре¬шению на разных уровнях проектирования, которые описываются системами нели¬нейных интегро-дифференциальных уравнений, не могут быть решены посредством существующего аналитического аппарата. А это означает, что решение требуемых за¬дач оказывается возможным лишь на основе приближенных, или, как их чаще при¬нято сейчас называть, численных методов математики.
Первоначально, объектом исследования вычислительной математики были разностные схемы, с помощью которых в том или ином приближении осуществля¬лось решение любого порядка дифференциальных уравнений. В первую очередь речь идет о конечных и разделенных разностях [1, 2], которые с успехом применялись в качестве т.н. метода сеток в решении многих задач математической физики [3, 4] и теории упругости [5]. Впрочем, не следует забывать, что при моделировании дву¬мер¬ных скалярных полей разностными схемами допускалась гораздо большая ошибка, чем при приближении теми же схемами одномерных одноименных полей.
Во второй половине XX века был разработан метод конечных элементов применительно к решению плоской задачи теории упругости [6]. Эта модификация позволила многим успешнее решать плоские задачи теории упругости, чем метод сеток, который с этим нововведением обнаружил всю свою несостоятельность в свя¬зи с исследованием двумерных скалярных полей [6, 7]. В дальнейшем были сфор¬му¬ли¬ро¬ваны пространственная и осесимметричная задачи метода конечных элемен¬тов [6, 7], однако этих модификаций не было достаточно для того, чтобы утверждать о том, что развитый метод конечных элементов суть естественное обобщение метода конечных (разделен¬ных) разностей на пространства бóльших размерностей.
После появления на свет публикации [8], в которой была сформулирована и решена методом конечных элементов антиплоская задача теории упругости, ока¬за¬лось вполне очевидным, что весь арсенал сформированного метода конечных эле¬ментов представляет собой естественное обобщение метода разделенных разнос¬тей. Более того, решение плоской и антиплоской задач теории упругости на основе МКЭ дало возможность численно решать практически любые задачи двумерных скаляр¬ных полей на плоскости. Эти две задачи образуют в совокупности применительно к изгибу тонких гибких пластин комплексную задачу МКЭ [9]. Из этого следует, что антиплоская задача МКЭ вытекает из комплексной задачи МКЭ как частный случай, адаптированный для случая изгиба тонких жестких пластин [9].
Теперь выясняется, что для полноценного численного дифференцирования дифференциальных уравнений, описывающих задачи двумерных скалярных полей, недостает только инструмента для исследо¬ва¬ния напряженно-деформированного со¬стояния произвольной конфигурации двумерных скалярных полей. Эта модифика¬ция предложена в публикации [10], где приводится распростране¬ние вышеуказанной ан¬типлоской задачи МКЭ на случай тонких жестких и гибких оболочек.
Таким образом, мы имеем в порядке возрастающей конкретности три вло¬жен¬ных в известном смысле задач, подлежащих детальному освещению:

• плоская задача метода конечных элементов;
• антиплоская задача МКЭ для плоских скалярных полей;
• комплексная задача МКЭ для плоских скалярных полей;
• комплексная задача МКЭ для произвольных скалярных полей.

 
2. Плоская задача метода конечных элементов

В обеих задачах теории упругости – о плоском на¬пря¬жённом и плоском деформированном состояниях [6] – поле перемещений од¬но¬значно определяется перемеще¬ния¬ми u и v в направлениях осей x и y прямоугольной сис¬те¬мы координат. В обоих случаях рассмат¬ри¬вают¬ся только по три компоненты напряжения и деформации в плос¬кос¬ти xOy. В случае плоского напряжённого состояния все ос¬таль¬ные компоненты напряжения равны нулю, а в случае плоского деформи¬ро¬ванного состояния напряжение в на¬прав¬лении, перпендикулярном плос¬кости xOy, не равно нулю, и может быть определено при необ¬хо¬ди¬мости через зна¬чения главных компонент напряжения [6].
Рассмотрим тонкую пластину, находящуюся в усло¬в謬ях об¬об¬щённого плос¬кого напряжённого состояния [6]. Пластина мыслен¬но разбивается на треугольные конеч¬ные элементы, после чего выде¬ляется один из них с узлами l,m,n [1] (рис. 1).
Перемещения каждого узла конечного элемента lmn, напри¬мер l, имеют две компоненты:

(2.1)

откуда следует, что вектор узловых перемещений эле¬мен¬та опреде¬ляет¬ся следующим шестимерным вектор-столб¬цом [6]:

(2.2)


Рис. 1. Формирование функций форм в пределах конечного элемента

Пусть перемещения произвольной точки внут¬ри эле¬мента однознач¬но определяются через узловые пере¬м嬬ще¬ния сле¬дую¬щим образом [2]:

(2.3)

где элементы прямоугольной матрицы размерности яв¬ляют¬ся функциями координат рассматрива¬е¬мой точки. Эти функ¬ции должны быть выбраны так, что¬бы при подстановке координат уз¬лов элемента в зависимость (2.3) они обра¬ща¬лись бы в соот¬вет¬ст¬вую¬щие узловые пере¬ме¬щения.
Для рассматриваемой плоской задачи МКЭ пере¬ме¬ще¬ния при¬ни-маются линейными относительно [2]:

(2.4)

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

(2.5)

Постоянные коэффициенты системы уравнений (2.5) мо¬гут быть оп¬ре¬делены, напри¬мер, по правилу Крамера:

(2.6)

где – площадь треугольника ijk на рис. 6, т.е.

(2.7)

Подставляем значения коэффициентов, выража¬емые ра¬вен¬ст¬вами (2.6), в за¬ви¬симости (2.4); полученные фор¬му¬лы [6] для ком¬по¬нент перемещения произ-вольной точки кнеч¬но¬го элемента будем на¬зы¬вать функциями перемеще¬ний [6]:

(2.8)

где

(2.9)

а остальные коэффициенты находятся путём круговой пе¬ре¬становки индексов , т.е.




Зависимость (2.8) можно также представить в мат¬рич¬ной форме [2]:

(2.10)

где

(2.11)

причём интерполяционные функции в (2.10) час¬то назы¬вают¬ся функциями формы [6].
Функциям формы можно дать любопытное геометри¬ческое истолкование. Нетрудно удостовериться, что для функ¬ции , напри¬мер, справедливо следующее представ¬л嬬ние [6] (рис. 1):

(2.12)

откуда следует, что функция формы принимает значе¬ние 1 в узле l и обращается в нули в узлах m и n. Свой¬ст¬вами, аналогичными (2.12), обладают также функции фор¬мы и ; как видно из рис. 1, они связаны между со¬бой очевидным тождеством:

(2.13)

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

(2.14)

В силу определённости сформулированных равенств (2.10) и (2.14) будем иметь [6]:

(2.15)

где

(2.16)

При принятых выражениях (2.10) для функций пере¬ме¬щений , как и сле¬довало ожидать, матрица , уста¬нов¬ленная равенством (2.15) с учётом (2.16), не зависит от координат точек внутри конеч¬но¬го элемента, т.е. дефор¬ма¬ции в его точках по¬стоянны [6].
Компоненты напряжения, как известно, связаны с ком¬понен¬тами деформа¬ции посредством закона Гука; эта зависимость запи¬сы¬вает¬ся в матричном виде для плоской задачи следующим образом:

(2.17)

В формуле (2.17) – матрица упругих постоянных ма¬те¬риала, кото¬рая в случае однородного и изотропного тела формулируется [6]:

(2.18)

Для задачи плоской деформации матрица упругих пос¬тоян¬ных (2.18) принимает вид [6]:

(2.19)

В случае обобщённого плоского напряжённого со¬сто¬яния будем иметь [6]:

(2.20)

Что касается вектора в зависимости (2.17), то он выра¬жает т.н. началь¬ную деформацию, т.е. деформацию, не зависящую от на¬пряжений. Чаще всего на практике к на¬чальным деформациям при¬во¬дят колебания темпе¬ра¬ту¬ры: тогда под на¬чальной деформацией бу¬дет по¬ни¬маться тем¬пературная деформация. В случае плоского на¬пря¬¬жён¬ного состояния изотропного материала для на¬гре¬того до тем¬пе¬ратуры T элемента при коэффициенте ли¬ней¬ного рас-ширения бу¬дем иметь:

(2.21)

В случае плоского деформированного состояния в изо¬троп¬ном ма¬те¬риа¬ле величина температурной деформации за¬висит от упругих по¬стоян¬ных, т.е.

(2.22)

При отсутствии начальной деформации, т.е. при , за¬ви¬симость (2.17) принимает вид

(2.23)

Путём исключения из равенства (2.23) с помощью соотно¬ше¬ния (2.15) вектор-столбца находим:

(2.24)

где – матрица, которую называют матрицей напря¬же¬ний плоской задачи [6].
Итак, напряжённо-деформированное состояние конечного элемен¬та l,m,n (рис. 1), описываемое зависимостями (2.15) и (2.17), рас¬смат¬ри¬вается в МКЭ как результат действия узловых сил, которые дол¬жны быть статически экви¬ва¬лент¬ны на¬пряжениям на границе ко¬неч¬ного эле¬мен¬та.