Аналитические и численные методы решения задач математической физики

Федеральное агентство по образованию Государственное образовательное учреждение высшего профессионального образования Тульский государственный университет КАФЕДРА "СТРОИТЕЛЬСТВО, СТРОИТЕЛЬНЫЕ МАТЕРИАЛЫ И КОНСТРУКЦИИ" Теличко Виктор Григорьевич ассистент АНАЛИТИЧЕСКИЕ И ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ МАТЕМАТИЧЕСКОЙ ФИЗИКИ КОНСПЕКТ ЛЕКЦИЙ Строительные конструкции для студентов, обучающихся по направлению 550100 – «Строительство» очной формы обучения (третий уровень высшего образования – уровень магистра) Количество часов 154 Тула 2005 СОДЕРЖАНИЕ Лекция № 1 Лекция № 2 Лекция № 3 Лекция № 4 Лекция № 5 Лекция № 6 Лекция № 7 Лекция № 8 Лекция № 9 Лекция № 10 Лекция № 11 Лекция № 12 Лекция № 13 Лекция № 14 Лекция № 15 Лекция № 16 Лекция № 17 Лекция № 18 ЛЕКЦИЯ № 1 ПЛАН 1.1. Положительные операторы и энергия. Оператор краевой задачи. Положительные и положительно определенные операторы 1.2. Энергетическое пространство 1.3. Обобщенные производные. Теоремы вложения. Неравенство Пуанкаре Оператор краевой задачи. Всякая краевая задача математической физики может быть сведена к уравнению вида Au = f , (1) где u – искомый элемент некоторого функционального пространства, A – данный оператор и f данный элемент того же или иного функционального пространства. Оператор A будем называть оператором данной краевой задачи. Положительные и положительно определенные операторы. Симметричный оператор A , действующий в некотором гильбертовым пространстве, называется положительным, если для любого элемента из области определения оператора, справедливо неравенство ( Au, u ) ≥ 0, причем знак равенства имеет место только тогда, когда u = 0 , т.е. когда u – нулевой элемент пространства. Если A положительный оператор, то скалярное произведение ( Au , u ) называется энергией элемента u по отношению к оператору A . Эта терминология оправдана тем, что во всех случаях, когда элемент u можно трактовать как смещение некоторой системы, величина ( Au, u ) совпадает, при подходящем выборе единиц измерения, с потенциальной энергией деформации этой системы. Теорема 1. Если оператор A положителен, то уравнение (1) имеет не более одного решения. Симметричный оператор A называется положительно определенным, если существует такая положительная постоянная γ 2 , что для любого элемента u из области определения оператора A справедливо неравенство: ( Au, u ) ≥ γ 2 Энергетическое пространство. u 2 Со (2) всяким положительным (в частности, со всяким положительно определенным) оператором можно связать особое гильбертово пространство, которое мы будем называть энергетическим пространством. Пусть A – положительный оператор, действующий в некотором гильбертовом пространстве H , и пусть M = D ( A ) – область определения это оператора. Введем на M новое скалярное произведение (которое будем обозначать квадратными скобками): если u и v элементы M , то положим [u, v ] = ( Au, v ) . Величину [u , v ] Доказывается, что назовем энергетическим энергетическое (3) произведением произведение u и удовлетворяет v. всем аксиомам скалярного произведения. Раз на множестве M введено новое скалярное произведение, это множество стало гильбертовым пространством. В общем случае неполное – пополним его. Построенное таким образом новое гильбертово пространство назовем энергетическим пространством и обозначим его через H A . Норму в энергетическом пространстве назовем энергетической нормой и будем обозначать символом u . Для элементов области M определения оператора A энергетическая норма определяется формулой u = ( Au, u ). (4) Сходимостью в энергетическом пространстве называется сходимость по энергии. Важную роль играет вопрос о природе тех элементов, которые служат для пополнения энергетического пространства. Если оператор A положительно определенный, то имеет место теорема, в силу которой все элементы пространства H A принадлежат также исходному гильбертову пространству H ; если u – элемент пространства H A , то имеет место неравенство u ≤ где символ 1 γ (5) u, означает норму в исходном пространстве H , а γ – постоянная входящая в неравенство (2). Если A – положительно определенный оператор, то их сходимости некоторой последовательности по энергии вытекает также ее сходимость в норме пространства: если u n ∈ H A , u ∈ H A и u n − u → 0 , то и u n − u → 0 . Обобщенные производные. Пусть Ω – некоторая конечная область m -мерного пространства и S – ее граница. Пограничной полосой области Ω называется совокупность тех ее точек, расстояния которых до границы S не превосходят некоторого заданного числа δ ; это число называется шириной пограничной полосы. Область Ω′ называется подобластью Ω , если все точки Ω′ принадлежат Ω , и внутренней подобластью, если Ω содержит не только все точки области Ω′ , но и все точки ее границы, иначе говоря, если существует пограничная полоса области Ω , не содержащая точек подобласти Ω′ . Обозначим через Φk множество функций, k раз непрерывно дифференцируемых в Ω и обращающихся в нуль в пограничной полосе этой области. Допустим сперва, что функция u ( x ) имеет в замкнутой области Ω = Ω + S непрерывную производную ∂ ku ; ∂ x i1 ∂ x i 2 ...∂ x i k одновременно она имеет и все предшествующие этой производные, также непрерывные в Ω . Пусть ϕ ( x ) – любая функция множества Φ k . На поверхности S функция ϕ ( x ) и все ее производные обращаются в нуль, поэтому интегрирование по частям приводит к формуле ∂ kϕ ∂ ku k ∫ u ∂ xi ∂ xi ...∂ xi dx = ( −1) Ω∫ ϕ ∂ xi ∂ xi ...∂ xi dx. Ω 1 2 k 1 2 k Пусть теперь u ( x ) – некоторая функция, суммируемая в любой внутренней подобласти Ω , и пусть существует такая функция w ( x ) , суммируемая в любой внутренней подобласти Ω , что для любой функции ϕ ( x ) ∈Φ k справедливо тождество ∂ kϕ k u dx 1 = − ( ) ∫ ∂ x i ∂ x i ...∂ x i ∫ ϕ w dx. Ω Ω 1 2 k Тогда w ( x ) называется обобщенной производной порядка k по переменным x i1 , x i 2 ,..., x i k от функции u ( x ) в области Ω . Обозначается обобщенная производная обычным символом ∂ ku w( x ) = . ∂ x i1 ∂ x i 2 ...∂ x i k Можно дать другое определение (6) обобщенной производной, равносильное приведенному выше. Пусть u ( x ) и w ( x ) суммируемы в любой внутренней подобласти обобщенной Ω′ производной последовательность k области вида (6) Ω . Функция от u ( x) , w( x) называется если существует раз непрерывно дифференцируемых внутри Ω функций u n ( x ) таких, что u n → u и ∂ k un ∂ x i1 ∂x i 2 ... ∂x i k →w в смысле метрики пространства L ( Ω′ ) , иначе говоря, если lim ∫ u n ( x ) − u ( x ) dx = 0 n →∞ Ω′ и lim ∫ n →∞ Ω′ ∂ k un ∂ x i1 ∂x i 2 ... ∂x i k − w ( x ) dx = 0; через Ω′ обозначена любая внутренняя подобласть области Ω . Теоремы вложения. Область называется звездной относительно данной точки, если любой луч, исходящий из этой точки, только один раз пересекается с границей области. Так, например, всякая выпуклая область звездна относительно любой своей точки. Здесь и ниже будут рассматриваться только конечные области, которые можно представить в виде суммы конечного числа взаимно налегающих подобластей, каждая из которых – звездная относительно любой точки некоторого шара. Теорема 2. Пусть функция u ( x ) суммируема в области Ω и имеет всевозможные обобщенные производные некоторого порядка k > 1 , также суммируемы в Ω . Тогда u ( x ) имеет суммируемые в Ω всевозможные обобщенные производные порядков, меньших k . Пространство W p( ) ( Ω ) С. Л. Соболева, где p > 1 , состоит из функций, l которые имеют в Ω всевозможные обобщенные производные порядка l , причем как сами функции, так и их производные порядка l суммируемы в Ω со степенью p ; норма в этом пространстве задается формулой p u p = ∫ u dx + ∫ ∑ p W p( ) ( Ω ) l Ω Ω ∂lu dx, ∂x i1 ∂x i 2 ...∂x i l (7) где сумма распространена на всевозможные выборы равных или не равных между собой чисел i1 , i 2 , ..., i l , каждое из которых не превосходит числа m – размерности области Ω . Теорема 3. Если u ( x ) ∈W p( ) ( Ω ) и pl > m , то u ( x ) эквивалентна l функции, непрерывной в замкнутой области Ω . Имеет место неравенство u C(Ω) ≤M u W p( ) ( Ω ) l (8) , где C ( Ω ) обозначает пространство функций, непрерывных в замкнутой области Ω , а M – постоянная, которая не зависит от функции u ( x ) . Всякое множество, ограниченное в W p( ) ( Ω ) , компактно в C ( Ω ) . l Следствие. Если u ( x ) ∈W p( ) ( Ω ) и натуральное число k таково, что l p ( l − k ) > m , то u ( x ) ∈ C ( k) (Ω) (то есть производные порядка k от u ( x ) непрерывны в Ω ), при этом u k C( ) (Ω) ≤M u W p( ) ( Ω ) l и всякое множество, ограниченное в W p( ) ( Ω ) , компактно в C ( l k) (Ω) . Теорема 4. Пусть Γ s , s ≤ m , – достаточно гладкое многообразие, все точки которого принадлежат Ω . Если pl ≤ m и s > m − pl , то всякая функция u ( x ) ∈W p( ) ( Ω ) эквивалентна функции, которая определена почти всюду на l Γ s с любой степенью q , удовлетворяющей неравенству q≤ ps ; m − pl (9) при этом имеет место неравенство u L q( Γ ) ≤ M 1 u s W p( ) ( Ω ) l , (10) где постоянная M 1 не зависит от функции u ( x ) . Если q удовлетворяет строгому неравенству q< ps , m − pl то множество, ограниченное в W p( ) ( Ω ) , компактно в L q ( Γ s ) . l Следствие. Если u ∈W p( ) ( Ω ) и k < l , то l (11) u ∈Wq( l −k ) (Ω), где q≤ mp , m − p (l − k ) (12) и справедливо неравенство u Wq( l −k ) ≤ M2 u W p( ) l (13) , где постоянная M 2 не зависит от функции u ( x ) . Теорема 5. Пусть N – число линейно независимых полиномов степени ≤ l − 1 от m переменных x1 , x 2 , ..., x m . Пусть линейные ограниченные в W p( ) ( Ω ) функционалы l j u ( j = 1, 2,..., n ) l таковы, что они не обращаются одновременно в нуль ни на одном полиноме степени ≤ l − 1 . Тогда норма, определяемая формулой u p N = ∑ l ju + ∫ ∑ j =1 p Ω ∂lu dx, ∂ x i1 ∂ x i 2 ... ∂ x i l (14) эквивалентна норме (7). Отметим два наиболее простых и важных неравенства, вытекающих из теорем вложения. 1. Неравенство Пуанкаре. Если u ( x ) ∈W2( ) ( Ω ) , то 1 2 ⎛ ⎞ ′ ≤ + u dx A grad u dx A u dx ( ) ⎜ ∫ ∫ ∫ ⎟⎠ . ⎝Ω Ω Ω 2 2 (15) 2. Если u ( x ) ∈W2( ) ( Ω ) и u = 0 на части границы S , то 1 2 2 ∫ u dx ≤ B ∫ ( grad u ) dx. Ω (16) Ω Если u = 0 на всей границе S , то неравенство (16) называют неравенством Фридрихса. В неравенствах (15) и (16) A , A′ , B – постоянные, которые не зависят от выбора функций u ( x ) , но зависят от выбора области Ω . ЛЕКЦИЯ № 2 ПЛАН 1.1. Функционал энергетического метода 1.2. Построение решения вариационной задачи 1.3. Метод Ритца. Методы решения системы Ритца Функционал энергетического метода. Если оператор A положителен, то решение уравнения (1) можно свести к решению некоторой вариационной задачи, как это вытекает из следующей теоремы. Теорема 6. Пусть A – положительный оператор. Если уравнение Au = f имеет решение, то это решение сообщает функционалу F ( u ) = ( Au, u ) − 2 ( u, f ) (17) наименьшее значение. Обратно, если существует элемент, реализующий минимум функционала (17), носит в литературе название энергетического метода. Функционал (17) будем называть функционалом энергетического метода. Теорема 5 не дает указаний ни на условия существования решения вариационной задачи, ни за то, как такое решение можно строить. Такие указания могут быть даны, если оператор краевой задачи положительно определенный. В этом случае введем в рассмотрение энергетическое пространство H A . По формуле (4) ( Au , u ) = u . Далее, по неравенству Коши2 Буняковского и неравенству (5) ( u, f ) ≤ f u ≤ f γ 2 u . Это означает, что линейный функционал ( u , f ) ограничен в H A ; по теореме Риса существует элемент u 0 ∈ H A такой, что ( u, f ) = ⎡⎣u, u 0 ⎤⎦ , если только u ∈ H A . Теперь функционал (17) приводится к виду 2 2 F ( u ) = u − 2 ( u , f ) = u − 2 ⎡⎣u , u 0 ⎤⎦ = u − u 0 − u 0 . 2 2 (18) Из формулы (18) вытекают два простых и важных следствия: 1) эта формула позволяет определить функционал F ( u ) не только на всех элементах области определения оператора A , но и на всех элементах энергетического пространства H A ; 2) в пространстве H A функционал F ( u ) достигает минимума при u = u 0 . Если u 0 ∈ D ( A ) , то по теореме 5 u 0 есть решение уравнения Au = f ; однако энергетическое пространство H A , вообще говоря, шире, чем D ( A ) , и может случиться, что элемент u 0 , построенный по теореме Риса и реализующий в энергетическом пространстве минимум функционала F ( u ) , не попадет в D ( A ) . В этом случае можно рассматривать u 0 как обобщенное решение уравнения Au = f . Построение решения вариационной задачи. Здесь мы укажем два общих приема построения элемента u 0 , реализующего минимум функционала F ( u ) ; весьма важный метод Ритца, являющийся конкретизацией каждого из упомянутых выше методов. 1) Допустим, что в пространстве HA существует полная ортонормированная (в смысле метрики этого пространства) счетная система ω n ( n = 1, 2, ...) , так что ⎧1, n = k , ⎡⎣ω n , ω k ⎤⎦ = ⎨ ⎩0, n ≠ k . Тогда ∞ u 0 = ∑ ( f , ω n )ω n . (19) n =1 Ряд (19) представляет собой ряд Фурье (в энергетическом пространстве) элемента u 0 по ортонормированной системе {ω n } ; этот ряд сходится как в метрике пространства H A (т.е. по энергии), так и в метрике пространства H . 2) Пусть u n – какая угодно минимизирующая последовательность для функционала F ( u ) . Тогда u n → u 0 по энергии, а также, следовательно, и в метрике исходного пространства. Любой достаточно далекий член минимизирующей последовательности можно рассматривать как приближенное решение задачи о минимуме функционала (18). Метод Ритца В пространстве H A выберем последовательность элементов ϕ 1 ,ϕ 2 ,...,ϕ n , ..., (20) удовлетворяющих двум условиям: 1) при любом n элементы ϕ 1 ,ϕ 2 ,...,ϕ n линейно независимы, 2) последовательность (20) полна по энергии; под этим мы понимаем следующее: каковы бы ни были элемент u ∈ H A и число ε > 0 , можно найти такое натуральное число N и такие постоянные α 1 , α 2 , …, α N , чтобы выполнялось неравенство N u − ∑ α kϕ k < ε . k =1 Элементы (20) называются координатными. По методу Ритца задаются натуральным числом n и строят приближенное решение u n вариационной задачи в виде n u n = ∑ a kϕ k , (21) k =1 где a k – постоянные, которые выбираются так, чтобы величина F ( u n ) была минимальной; для определения этих постоянных получается линейная алгебраическая система ⎡⎣ϕ 1 ,ϕ 1 ⎤⎦ a1 + ⎡⎣ϕ 2 ,ϕ 1 ⎤⎦ a 2 + ... + ⎡⎣ϕ n ,ϕ 1 ⎤⎦ a n = ( f ,ϕ 1 ) , ⎫ ⎪ ⎡⎣ϕ 1 ,ϕ 2 ⎤⎦ a1 + ⎡⎣ϕ 2 ,ϕ 2 ⎤⎦ a 2 + ... + ⎡⎣ϕ n ,ϕ 2 ⎤⎦ a n = ( f ,ϕ 2 ) , ⎪⎪ ⎬ ....................................................................................⎪ ⎪ ⎡⎣ϕ 1 ,ϕ n ⎤⎦ a1 + ⎡⎣ϕ 2 ,ϕ n ⎤⎦ a 2 + ... + ⎡⎣ϕ n ,ϕ n ⎤⎦ a n = ( f ,ϕ n ) , ⎪ ⎭ (22) называется системой Ритца. В силу линейной независимости элементов ϕ 1 ,ϕ 2 ,...,ϕ n определитель системы Ритца отличен от нуля и эта система единственным образом разрешима. Метод Ритца находится в тесной связи с представлением решения в виде ряда (19). Последовательность (20) координатных элементов подвергнем процессу ортогонализации в метрике пространства H A , что приведет нас к некоторой полной ортонормированной в H A системе {ω n } . Теперь искомый элемент u 0 можно представить в виде ряда (19). Оказывается, что приближенное решение u n , построенное по методу Ритца с помощью формул (21) и (22), есть n − я частичная система ряда (19). Последовательность приближенных решений, построенных по методу Ритца, есть также минимизирующая последовательность для функционала (18). С возрастанием n энергетические нормы приближенных решений по Ритцу возрастают (точнее, не убывают) и стремятся к энергетической норме точного решения u k ≤ u n , k < n, lim u n = u 0 . n→∞ Далее, имеет место тождество 2 2 2 u0 − un = u0 − un . (23) Условимся погрешность приближенного решения величиной u0 − un . Из формулы (23) следует, что un оценивать погрешность приближенного решения по Ритцу убывает (точнее, не возрастает) при добавлении новых координатных элементов. Если координатные элементы принадлежат не только пространству H A , но также и области определения данного оператора A , то ϕ j , ϕ k = ( Aϕ j ,ϕ k ) и систему Ритца можно представить в виде ( Aϕ ,ϕ ) a + ( Aϕ ,ϕ ) a + ... + ( Aϕ ,ϕ ) a = ( f ,ϕ ) , ( Aϕ ,ϕ ) a + ( Aϕ ,ϕ ) a + ... + ( Aϕ ,ϕ ) a = ( f ,ϕ ) , ⎫ ⎪ ⎪⎪ 1 2 1 2 2 2 n 2 n 2 ⎬ ...........................................................................................⎪ ( Aϕ1,ϕ n ) a1 + ( Aϕ 2 ,ϕ n ) a 2 + ... + ( Aϕ n ,ϕ n ) a n = ( f ,ϕ n ) , ⎪⎪⎭ 1 1 1 2 1 2 n 1 n 1 (24) Методы решения системы Ритца. Матрица системы Ритца (22) положительно определенная, и для решения этой системы пригодны все методы, применяемые для решения систем линейных алгебраических уравнений матрицей, например метод итераций. с положительно определенной ЛЕКЦИЯ № 3 ПЛАН 1.1. Основы метода. Достаточный принцип сходимости 1.2. Применение к задачам математической физики 1.3. Обобщение метода Бубнова-Галеркина Пусть дано уравнение Au − f = 0 , где A – линейный оператор, действующий в некотором гильбертовом пространстве H; мы не предполагаем об операторе A ничего кроме его линейности. По методу Бубнова-Галеркина выбираем последовательность элементов ϕ n ∈ D ( A ) и ищем приближенное решение в виде n u n = ∑ a kϕ k (25) k =1 ak коэффициенты определяем из условия, чтобы левая часть данного уравнения после подстановки в нее u n вместо u оказалась ортогональной к элементам ϕ 1 ,ϕ 2 ,...,ϕ n . Это приводит к системе уравнений ( Aϕ ,ϕ ) a + ( Aϕ ,ϕ ) a + ... + ( Aϕ ,ϕ ) a = ( f ,ϕ ) , ( Aϕ ,ϕ ) a + ( Aϕ ,ϕ ) a + ... + ( Aϕ ,ϕ ) a = ( f ,ϕ ) , ⎫ ⎪ ⎪⎪ 1 2 1 2 2 2 n 2 n 2 ⎬ ...........................................................................................⎪ ( Aϕ1,ϕ n ) a1 + ( Aϕ 2 ,ϕ n ) a 2 + ... + ( Aϕ n ,ϕ n ) a n = ( f ,ϕ n ) , ⎪⎪⎭ 1 1 1 2 1 2 n 1 n 1 (26) по форме совпадающей с системой (24). Отсюда сразу следует, что метод Бубнова-Галеркина совпадает с методом Ритца, если A – положительно определенный оператор. Метод Бубнова-Галеркина можно применять и в задаче о собственных числах. Так, если потребуется найти собственные числа уравнения Au − λu = 0 , то по методу Бубнова-Галеркина их приближенные значения находятся как корни уравнения ( Aϕ ,ϕ ) − λ (ϕ ,ϕ ) ( Aϕ ,ϕ ) − λ (ϕ ,ϕ )...( Aϕ ,ϕ ) − λ (ϕ ,ϕ ) ( Aϕ ,ϕ ) − λ (ϕ ,ϕ ) ( Aϕ ,ϕ ) − λ (ϕ ,ϕ )...( Aϕ ,ϕ ) − λ (ϕ ,ϕ ) = 0. 1 1 1 1 2 1 2 1 2 2 1 2 2 2 1 n 2 1 n n 2 1 n 2 .......................................................................................... ( Aϕ ,ϕ ) − λ (ϕ ,ϕ ) ( Aϕ ,ϕ ) − λ (ϕ ,ϕ )...( Aϕ 1 n 1 n 2 n 2 n n (27) , ϕ n ) − λ (ϕ n , ϕ n ) Приближенные значения собственных чисел более общего уравнения Au − λ Bu = 0 находятся как корни уравнения ( Aϕ ,ϕ ) − λ ( Bϕ ,ϕ ) ( Aϕ ,ϕ ) − λ ( Bϕ ,ϕ )...( Aϕ ,ϕ ) − λ ( Bϕ ,ϕ ) ( Aϕ ,ϕ ) − λ ( Bϕ ,ϕ ) ( Aϕ ,ϕ ) − λ ( Bϕ ,ϕ )...( Aϕ ,ϕ ) − λ ( Bϕ ,ϕ ) = 0. (28) 1 1 1 1 2 1 2 1 2 2 1 2 2 2 1 n 2 1 n n 2 1 n 2 .......................................................................................... ( Aϕ ,ϕ ) − λ ( Bϕ ,ϕ ) ( Aϕ ,ϕ ) − λ ( Bϕ ,ϕ )...( Aϕ 1 n 1 n 2 n 2 n n ,ϕ n ) − λ ( Bϕ n ,ϕ n ) Достаточный признак сходимости. Предположим, что оператор A имеет вид A = A0 + B , где A0 – положительно определенный в данном гильбертовом пространстве H оператор, и D ( B ) ⊃ D ( A ) . Допустим еще, что оператор A0−1B вполне непрерывен в энергетическом пространстве H A0 . Координатные элементы ϕ n выберем так чтобы: 1) они входили в D ( A0 ) ; 2) взятые в любом конечном числе, они были линейно независимы; 3) их совокупность была полна в H A0 . В перечисленных условиях верны следующие теоремы. Теорема 7. Приближенные решения уравнения, построенные по методу Бубнова-Галеркина, существуют при достаточно больших n и сходятся в норме пространства H A0 к точному решению данного уравнения, если это последнее разрешимо и имеет только одно решение. Теорема 8. Приближенные собственные числа уравнения A0u − λ Bu = 0 , построенные по методу Бубнова-Галеркина, сходятся к соответствующим точным значениям этих чисел. Применение к задачам математической физики. Мы перечислим здесь некоторые задачи, к которым применимы теоремы 7 и 8. Уравнение ds ⎛ d su ⎞ ( −1) s ⎜ p ( x ) s ⎟ + Bu = f ( x ) dx ⎝ dx ⎠ s (29) При краевых условиях u ( a ) = u′ ( a ) = ... = u ( ( a ) = 0,⎫⎪ ⎬ s −1 u ( b ) = u′ ( b ) = ... = u ( ) ( b ) = 0 ⎪⎭ s −1) (30) Допускает применение метода Бубнова-Галеркина, если p ( x ) ≥ p 0 , где p 0 – положительная постоянная и B – дифференциальный оператор порядка ≤ 2 s − 1 с ограниченными коэффициентами. Метод Бубнова-Галеркина применим к невырождающемуся уравнению эллиптическому типа m ∂ 2u ∂u − ∑ A jk + ∑Bj + Cu = f ( x ) ∂x j ∂x k j =1 ∂x j j , k =1 m (31) при краевых условиях задачи Дирихле u S =0 (32) ⎡ m ⎤ ∂u cos (ν , x k ) + σ u ⎥ = 0; ⎢ ∑ A jk ∂x j ⎢⎣ j ,k =1 ⎥⎦ S (33) или смешанной задача Неймана (σ ≡ 0 ) не представляет исключения. Приближенные собственные числа, построенные по методу Галеркина, для уравнения ds ⎛ d su ⎞ ( −1) s ⎜ p ( x ) s ⎟ − λ Bu = 0 dx ⎝ dx ⎠ s При краевых условиях (30) сходятся к соответствующим точным собственным числам, если p ( x ) и B удовлетворяют тем же условиям, что и в уравнении (29). Тоже верно для собственных чисел уравнения ⎛ m ⎞ ∂ ⎛ ∂u ⎞ ∂u −∑ + Cu ⎟ = 0 ⎜⎜ A jk ⎟⎟ − λ ⎜⎜ ∑ B j ⎟ ∂x k ⎠ j , k =1 ∂x j ⎝ ⎝ j =1 ∂x j ⎠ m При краевых условиях (32) и (33). Обобщение метода Бубнова-Галеркина (проекционный метод) Пусть дано линейное уравнение Au = f , (34) где u ∈ E1 , u ∈ E 2 , E1 и E 2 – некоторые банаховы пространства и A – оператор, действующий из E1 в E2 . Зададим последовательность подпространств L n размерности n , где n = 1, 2,... , причем Ln ⊂ D ( A ) ⊂ E1 при любом n . Зададим, далее, последовательность n -мерных подпространств M n ⊂ E 2 , и пусть Pn – оператор проектирования из пространства E 2 в подпространство M n . Проекционный метод состоит в том, что точное уравнение (34) заменяется приближенным уравнением Pn Au n = Pn f , u n ∈ L n , (35) Очевидно, проектор Pn имеет вид n Pnu = ∑ l k ( u ) ψ k , k =1 где ψ 1 ,ψ 2 ,...,ψ n – базис подпространства M n , а l k ( u ) – линейные ограниченные в E 2 функционалы. Если ϕ 1 , ϕ 2 ,…, ϕ n есть базис для подпространства Ln , то u n имеет вид (25) и уравнение (35) сводится к линейной алгебраической системе n ∑ l ( Aϕ ) a k =1 j k k =lj( f ) ( j = 1, 2,..., n ) ; (36) при такой записи нет необходимости явно указывать подпространство M n достаточно указать функционалы l j . ЛЕКЦИЯ № 4 ПЛАН 1.1. Сведение решения краевой задачи для дифференциальных уравнений в частных производных к краевым задачам для обыкновенных дифференциальных уравнений 1.2. Метод Власова-Канторовича как вариант метода БубноваГалеркина Этот метод позволяет свести решение краевой задачи для дифференциальных уравнений в частных производных к краевым задачам для обыкновенных дифференциальных уравнений, т.е. понизить размерность задачи. Это оказывается очень удобным при решении задач динамики, нестационарных задач теплопроводности и в тех случаях, когда по одной из координат заданы сложные краевые условия. Если в методе Бубнова-Галеркина в решение входил неизвестный числовой параметр a k , то в методе Власова-Канторовича в решение входит неизвестная функция одной из переменных. Эта функция находится из условий, аналогичных тем, что и в методе Бубнова-Галеркина, только теперь получим не систему линейных алгебраических уравнений, а систему обыкновенных дифференциальных уравнений. Итак, рассмотрим следующую краевую задачу. Найти функцию u ( x, y , t ) , удовлетворяющую в области D : {0 ≤ x ≤ a; 0 ≤ y ≤ b; 0 ≤ t ≤ T } дифференциальному уравнению Lu = f ( x, y, t ) , На границе области ( Γ ) по пространственным координатам x, y – краевым условиям R [ u ] Γ = ϕ ( x, y ) в общем случае может быть ϕ = ϕ ( x, y, t ) , а по переменной t – начальным условиям. Приближенное решение возьмем в виде n u n ( x, y, t ) = ψ ( x, y ) + ∑ C i ( t ) ϕ i ( x, y ), (37) i =1 где C i ( t ) – неизвестные функции переменной t ; ϕ i ( x, y ) – известные (аппроксимирующие) функции, удовлетворяющие однородным краевым условиям, т.е. R ⎡⎣ϕ i ( x, y ) ⎤⎦ = 0 , удовлетворяющая однородным R ⎡⎣ψ ( x, y ) ⎤⎦ = ϕ ( x, y ) . Γ функциям а Γ ϕ i ( x, y ) Условие дает ψ ( x, y ) – краевым ортогональности систему обыкновенных известная функция, условиям, невязки Lu n − f т.е. к дифференциальных уравнений относительно функций C i ( t ) a b ∫ ∫ ( Lu n − f )ϕ j ( x, y ) dx dy = 0 ( j = 1,2,..., n ) (38) 0 0 Решив эту систему при заданных начальных условиях по переменной t , найдем искомые функции C i ( t ) , а подставив их в (37), найдем приближенное решение. ЛЕКЦИЯ № 5 ПЛАН 1.1. Варианты метода последовательных приближений. 1-й вариант метода последовательных приближений 1.2. 2-й вариант метода последовательных приближений 1.3. 3-й вариант метода последовательных приближений Метод возмущений или метод малого параметра, можно рассматривать как некоторую модификацию метода последовательных приближений, при применении которого для вычисления n -го приближения пренебрегают членами, имеющими порядок выше n -го. Изложим варианты метода последовательных приближений и метода малого параметра. 1. Пусть формула Bx = P символизирует совокупность системы дифференциальных уравнений и ее начальных или граничных условий, а x – совокупность подлежащих определению неизвестных функций. Представим предыдущие уравнения в виде Ax + ( B − A ) x = P, Ax + Cx = P, где A – некоторый линейный оператор, C = B − A . В качестве первого приближения примем решение уравнения Ax1 = P . Погрешность решения будет Δ 2 = P − Ax1 − Cx1 . Для устранения неуравновешенности Δ 2 определим поправку δ 2 к первому приближению из уравнения Aδ 2 = Δ 2 и в качестве второго приближения к истинному решению принимается x 2 = x1 + δ 2 . Аналогично определяется неуравновешенность второго приближения, новая поправка к решению и т.д. 2. Пусть дана система уравнений Ax + Cx = pP, где P – заданная функция координат, а p – число, характеризующее величину нагрузки. Зададимся каким-либо значением параметра p = p1 и определим x1 (например, прогиб) в первом приближении из уравнения Ax1 ( p1 ) = p1P. Далее внесем полученное решение в уравнения Aδ 2 ( p1 ) = p 2 P − Ax1 ( p1 ) − Cx1 ( p1 ) и определим поправку δ 2 ( p1 ) к первому приближению, а также величину второго приближения к решению x 2 ( p 1 ) = x1 ( p 1 ) + δ 2 ( p 1 ) . Число p 2 определятся из условия равенства прогибов центра оболочки в первом и втором приближениях. Неуравновешенность – неравенство левой и правой частей уравнений – устраняется при этом путем изменения приближенной формы прогиба и величины нагрузки, которая подбирается так, чтобы «неуравновешенность» как можно меньше сказывалась на прогибе оболочки. Таким образом, в качестве второго приближения к истинной зависимости между перемещением оболочки и ее нагрузкой следует принять ее параметрическое представление x ≈ x 2 ( p1 ) ; p = p 2 ( p1 ) Через параметр p1 . Для определения следующих поправок к параметрическому представлению искомой зависимости решаем систему уравнений Aδ 2 ( p1 ) = p 3 P − Ax 2 ( p1 ) − Cx 2 ( p1 ) и полагаем x 3 ( p1 ) = x 2 ( p1 ) + δ 2 ( p1 ) , причем число p 3 определяем из условия равенства поперечных прогибов центра оболочки для решений x 3 ( p1 ) и x 2 ( p1 ) . Задаваясь рядом последовательных значений p 2 , можно найти соответствующие друг другу значения x 3 и p 3 и построить графики их зависимости. 3. Для описания третьего метода последовательных приближений предварительно дадим некоторые общие соображения. Пусть Bω = P – нелинейное уравнение, связывающее между собой нагрузку и прогиб оболочки. Увеличив давление на бесконечно малую величину δ P , получим изменение прогиба, определяемого из уравнения B′ (ω ) = δω = δ P. Здесь B′ (ω ) – линейный относительно δω оператор, зависящий линейно от ω . Обозначим его собственные числа через λ1 , λ 2 , ...; λ1 ≤ λ 2 ≤ ..., а через ϕ 1 , ϕ 2 , … – собственные функции оператора B′ (ω ) , т.е. функции, удовлетворяющие уравнениям B′ ( ω ) ϕ i = λ i ϕ i ( i = 1, 2,...) , причем предполагаем, что эти функции образуют полную систему, нормированы и ортогональны, т.е. (ϕ (ϕ i ⋅ ϕ k ) = 1, если i = k , i ⋅ ϕ k ) = 0, если i ≠ k , где (ϕ i ⋅ ϕ k ) – скалярное произведение функций ϕ i и ϕ k . Разложим δ P в ряд по собственным функциям ∞ dP = ∑ d c i ϕ i , i =1 где коэффициенты dc i определяются непосредственно, если умножить скалярно уравнение на ϕ i и воспользоваться условием ортогональности dc i = ( d P, ϕ i ) ( i = 1, 2, ...) . Будем искать решение уравнения B′ (ω ) δ ω = δ P в виде ряда ∞ δω = ∑ δ y iϕ i , где d y i = (δ ω , ϕ i ) . i Внося δ ω и d P в наше уравнение и учитывая, что B′ (ω )ϕ i = λ i ϕ i , получим ∞ ∞ i =1 i =1 ∞ B′ ( ω ) ∑ d y i ϕ i = ∑ d c i ϕ i ; ∞ ∑ λ dy ϕ = ∑ d c ϕ i =1 i i i i =1 i Сравнивая коэффициенты при одинаковых собственных функциях обеих частей последнего уравнения, имеем dy i = dc i λi . Внося это в выражение для δ ω , получим ∞ dc i i =1 λi δω = ∑ ϕi. Рассматривая это выражение и считая, что ϕ i ограничены по совокупности, замечаем, что наибольшее влияние на изменение прогиба δ ω имеет величина dc1 , которая делится на наименьший по абсолютной величине делитель λ i . Заметим, что для многих задач теории пластин и оболочек величины λ i быстро возрастают по абсолютной величине при увеличении их номера, поэтому часто можно ограничиться только первым членом, полагая приближенно δω ≈ dc1 λ1 ϕi, например, при решении задачи о прогибах свободно опертой пластины под действием равномерно распределенной нагрузки методом Навье двойных тригонометрических радов первый член ряда дает значение прогиба центра пластины с погрешностью, не превышающей 2,6%. Этот метод отличается от второго тем, что при его применении на каждом этапе величина p корректируется так, чтобы в двух последовательных приближениях оставались одинаковыми не поперечные прогибы центра оболочки, а «обобщенные перемещения»: ⎛ δω ⎞ ⎛ δω ⎞ , ωn ⎟ = ⎜ , ω n+1 ⎟ , ⎜⎜ ⎟ ⎜ ∂ p1 ⎟ ⎝ dp1 ⎠ ⎝ ⎠ δω при этом вместо величины приближенное значение δ ωn ∂ p1 подставить вместо величины в эту формулу подставляется ее ∂ p1 . δ ωn ∂ p1 Для упрощения величины δ ω1 ∂ p1 вычислений можно . При этом сходимость может ухудшиться, но вычисления упрощаются. Особенностью этого метода является требование, чтобы скалярное произведение приближенное неизменным. решение уравнений, при изменении (ω ,ψ ) , n p где ψ – оставалось ЛЕКЦИЯ № 6 ПЛАН 1.1. Модификации метода последовательных приближений 1.2. Разложение по степеням малого параметра 1.3. Использование метода возмущений Методы возмущений или методы разложения по степеням малого параметра для решения уравнений типа Ax + μ Cx = pP, где μ – малый параметр, представляют собой видоизменение приведенных выше методов последовательных приближений и отличаются лишь тем, что при их применении во время выполнения n -ого шага вычислений пренебрегаются члены, содержащие μ в степени выше n -й. Этот метод применим лишь в тех случаях, когда оператор C аналитичен, т.е. когда он может представлен в виде степенного ряда ∞ Cx = ∑ C s ( x1 ... x s ), s =0 где C s ( x1 ... x s ) – полилинейные операторы n -ой степени. При этом подлежащие решению уравнения имеют вид ∞ Ax + μ ∑ C s ( x1...x s ) = pP, s =0 где C s ( x1 , x 2 ... x s ) аргументов – операторы, линейные относительно каждого из x1 , x 2 ,..., x s . Рассмотренная в следующем пункте система уравнений пологой оболочки принадлежит к этому классу уравнений. При этом в последней системе будут участвовать лишь линейные и билинейные операторы. Рассмотрим применение метода малого параметра, соответствующему третьему методу последовательных приближений. В качестве первого приближения выбирается решение уравнения Ax 0 = p 0 P, где p 0 – произвольно выбранное число. Решение задачи разыскивается в виде рядов x = x 0 + μ x1 + μ 2 x 2 + ... , p = p 0 + μ p1 + μ 2 p 2 + ... . Внося эти выражения в подлежащие решению уравнения, разлагая правую и левую части в ряды по степеням малого характера и сравнивая коэффициенты при одинаковых степенях μ , получим систему уравнений Ax1 + C 0 = p1 P, Ax 2 + C1 ( x1 ) = p 2 P, Ax 2 + ⎡⎣C 2 ( x1 x1 ) + C1 ( x 2 ) ⎤⎦ = p 2 P, ........................................................ ........................................................ Коэффициенты p1 ,..., p n определяются из условия (x + ... + μ ( n −1) ) x n−1 + μ n x n , x 0 = ( x 0 + ... + μ n−1 x n−1 , x 0 ) . Отсюда следует, что при n > 0 ( x n , x 0 ) = 0 . Итак, если определены x 0 ,..., x n−1 и p 0 ,..., p n−1 , то следующие коэффициенты x n и p n определяются путем решения системы уравнений Ax n + γ n ( x 0 ... x n−1 ) = p n P, ( x , x ) = 0. n Если применять вариант метода возмущений, соответствующий второму варианту метода последовательных приближений, то последнее уравнение заменяется условием равенства нулю функций x n ЛЕКЦИЯ № 7 ПЛАН 1.1. Поиск приближенного аналитического решения 1.2. Выбор точек коллокации 1.3. Рациональный выбор аппроксимирующих функций В этом параграфе излагается метод, дающий возможность найти приближенное решение краевой задачи в виде аналитического выражения. Пусть требуется найти методом коллокаций приближенное решение, удовлетворяющее внутри области D уравнению Lw = f (39) И граничным условиям на контуре Γ l i w = f i , i = 1, 2, 3, ... , k , (40) где L и l i – дифференциальные операторы. Для решения выберем некоторую функцию w , зависящую от n параметров a n : w = w ( x, y, a1 , a 2 , ..., a n ) . (41) Подставив (41) в (39) и (40), получим Δ = Lw − f 1 , (42) δ i = li w − f i . (43) Параметры a n определяем из условия обращения в нуль (43) и (42) на некоторой системе точек, называемых точками коллокаций, расположенных внутри области D и на ее границе Γ . Это условие приводит к алгебраической системе n уравнений относительно a n Δ M j = 0, j = 1, 2, ..., n1 , (44) δ i M α = 0, ∑ iα = n − n , (45) 1 где M j и M α – точки коллокаций внутри области D и на ее границе соответственно. Однако в таком общем виде метод коллокаций применяется редко, так как приходится вводить большое число параметров a n . Обычно w подбирают так, чтобы она удовлетворяла (40), а параметры определяются из (44). В математической литературе методу коллокаций уделено значительно меньше внимания, рационального чем выбора другим приближенным аппроксимирующих методам. функций, Вопросы выбор точек коллокации, сходимость и оценка погрешностей практически не разработан. Недостатком метода коллокаций является произвол в выборе точек коллокаций, следствие чего при одной и той же аппроксимирующей функции методом коллокации может быть получено решение, вообще говоря, менее точное, чем вариационным методом. Однако этот недостаток сглаживается при большом числе неизвестных параметров. Другой важный вопрос метода коллокаций – рациональный выбор аппроксимирующих функций. Например, М.С. Корнишин в качестве аппроксимирующих функций рекомендует брать точные решения родственных задач. Например, для решения двумерных задач выгодно брать прогиб в виде w = w1 ( x ) w2 ( y ) wk ( x, y, a1 , a 2 , ..., a n ) , где w1 ( x ) и w2 ( y ) – точные решения одномерных (46) задач, wk ( x, y, a1 , a 2 ,..., a n ) – корректирующая функция. При решении нелинейных задач функцию w можно брать в виде w = w0 ( x, y ) wk ( x, y, a1 , a 2 ,..., a n ) , (47) где w0 ( x, y ) – точное решение линейной задачи, wk ( x, y, a1 ,..., a n ) – корректирующая функция. Функцию w можно задавать и в виде суммы n w = w0 ( x, y ) + ∑ a i wi ( x, y ), (48) i =1 причем если решение родственной задачи неизвестно или имеет сложную структуру, то за w0 ( x, y ) принимают простую функцию, близкую к искомой. ЛЕКЦИЯ № 8 ПЛАН 1.1. Значение расчета за пределом упругости 1.2. Особенности нелинейной работы материала 1.3. Нелинейно-упругие системы Большинство строительных материалов или совсем не подчиняется закону Гука, или подчиняется ему лишь при относительно малых напряжениях, не превосходящих предела упругости материала. Между тем в работе инженерных сооружений довольно часто напряжения в отдельных точках или даже в целой области существенно превосходят предел упругости и приближаются к пределу прочности материала. Поэтому очень важно уметь хотя бы приближенно рассчитать конструкции и сооружения в неупругой стадии их работы. Особенно большое значение такой расчет имеет для оценки действительной несущей способности конструкций и для выяснения условий, при которых может произойти их разрушение. Однако расчет конструкций за пределом упругости является значительно более сложным, чем расчет в предположении линейной зависимости между напряжениями и деформациями, а многие задачи такого рода до сих пор еще не могут считаться вполне решенными. Это является одной из причин, почему в практике проектирования пока еще используют главным образом простые и хорошо разработанные методы линейной строительной механики, основанной на законе Гука. Другая, весьма существенная причина широкого применения расчета по упругой стадии работы материала заключается в том, что такой расчет полностью гарантирует безопасность конструкций, как при однократных, так и при многократных и переменных силовых воздействиях, тогда как разработанная к настоящему времени методика расчета неупругих систем пока еще не дает точного ответа, в каких случаях повторные и длительные нагружения не будут вызывать накопление незаметных на глаз повреждений, создающих повышенную опасность разрушения всей конструкции. Таким образом, часто приходится мириться с излишними запасами прочности, создаваемыми расчетом по упругой стадии, и идти на заведомый перерасход материалов. Расчет за пределом упругой работы конструкции имеет целью, насколько это возможно, гарантированно уменьшить этот перерасход, получив, как предсказывает теория, существенную экономию. Особенности нелинейной работы материала Основное, что отличает неупругую работу материала от упругой – это отсутствие потенциала внутренних сил, обусловленное частичным рассеиванием механической энергии, которая переходит в другие виды энергии, главным образом в теплоту, и полностью не возвращается обратно при разгрузке элементов конструкции. Это наглядно видно на обычной диаграмме работы образца при одноосном растяжении или сжатии (рис. 1), т. е. зависимости между деформациями и растяжениями. Форма кривой этой зависимости отвечает определенному режиму нагружения, принятому при испытаниях. При иных скоростях нагружения, остановках на определенных этапах и т. п. вид диаграммы работы несколько меняется. При переходе же от нагружения к разгрузке зависимость между деформациями и напряжениями приобретает совершенно иной вид, показанный на рис. 2 штриховой линией. Легко видеть, что при одной и той же величине деформации могут быть различные напряжения в образце. Следовательно, между напряжениями и деформациями за пределом упругости нет функциональной связи. Еще в большей степени это присуще двух- и трехосному деформированию, в частности, когда компоненты сложного напряженного состояния меняются во времени пропорционально друг другу – так называемое с л о ж н о е н а г р у ж е н и е. Таким образом, нельзя говорить о потенциальной энергии внутренних сил, как функции компонентов деформаций за пределом упругости. Точно также для системы, составленной из неупругих элементов, нельзя говорить о потенциальной энергии внутренних усилий, как функции компонентов вектора перемещений этой системы. Это обстоятельство, как мы увидим ниже, обходится в расчетах за пределом упругости при определенных ограничениях в режиме работы конструкции. σ σ ε Рис. 1 σ O ε Рис. 2 σ A C B ε Рис. 3 ε Рис. 4 Определим работу, совершаемую продольной силой в стержне при деформировании за пределом упругости с определенным режимом нагружения, который принят для получения стандартной диаграммы работы материала. При бесконечно малой деформации эта работа dA = − N d λ = − Flσ d ε , где N – продольная сила; λ – абсолютная деформация растяжения; F – площадь сечения стержня; l – его длина; σ – продольное напряжение; ε – относительная деформация растяжения. При монотонном увеличении деформации от нуля до некоторого значения ε работа ε A = − Fl ∫ σ d ε . Поскольку диаграмма работы материала при заданном монотонном увеличении нагрузки определяет напряжение как функцию деформаций σ = σ (ε ), то работа A в данном случае тоже будет представлять собой однозначную функцию деформации A = A(ε ). Геометрически эта работа может быть выражена умноженной на − Fl площадью OAB ,ограниченной осью абсцисс и кривой диаграммы работы материала (рис. 3). При разгрузке продольная сила выполняет работу на отрицательных приращениях деформаций же согласно кривой разгрузки материала. Эта работа может быть выражена площадью ABC , ограниченной осью абсцисс и кривой разгрузки. Как видим, возвращенная при разгрузке работа меньше, чем работа, произведенная во время нагружения. Разность затраченной и возвращенной работы соответствует площади OAC (рис. 3). Она представляет собой поглощенную механическую энергию, перешедшую в теплоту или затраченную на структурные преобразования материала. При повторных нагрузках кривые нагружений и разгрузок приближаются все ближе друг к другу, однако полного совпадения не достигают, образуя так называемые п е т л и г и с т е р е з и с а, выражающие потери механической энергии при каждом цикле нагружения и разгрузки (рис. 4). Нелинейно упругие системы Гипотетически можно представить себе материал с нелинейной зависимостью напряжений от деформаций, который при разгрузке полностью возвращает накопленную механическую энергию. Диаграмма разгрузки для такого материала совпадает с диаграммой нагружения, и поглощения механической энергии не происходит. Если известно, что напряжения в конструкции в процессе ее нагружения нигде не будут уменьшаться, то можно условно считать, что конструкция выполнена из нелинейно упругого материала с функциональной зависимостью σ = f (ε ) , совпадающей с диаграммой работы материала при нагружении. При этом становится возможным использование понятия потенциальной энергии системы. Понятия нелинейно упругого материала и нелинейно упругой системы очень удобны для расчета неупругих систем при первичном их нагружении. Необходимо только следить, чтобы нигде не возникало уменьшения напряжений (об этом можно судить на основании расчета). При появлении уменьшения напряжений схема нелинейно упругой системы будет уже неприменима. В нелинейно упругих системах внутренние силы N j ( j = 1, 2, ..., m) , где m – число элементов системы, выражаются через деформации λ j в общем случае системой нелинейных уравнений: N1 = N1 (λ1 , λ2 , ..., λm ); ⎫ N 2 = N 2 (λ1 , λ2 , ..., λm ); ⎪⎪ ⎬ """"""""" ⎪ N m = N m (λ1 , λ2 , ..., λm ).⎪⎭ (49) Уже составление этой системы уравнений, а тем более решение ее совместно с условиями равновесия и условиями неразрывности деформаций является делом довольно сложным, однако при наличии современной электронно-вычислительной техники вполне разрешимым. Во многих случаях, например при расчете ферм, система уравнений (49) распадается на отдельные уравнения. Это бывает тогда, когда усилие в одном элементе зависит от деформации только этого элемента и не зависит от деформаций других элементов. Так, для ферм мы получаем N j = N j (λ j ) = Fj f (λ j l j ) ( j = 1, 2, ..., m), (50) где Fj и l j – площади и длины стержней фермы; f – функция, выражающая зависимость напряжений σ от деформаций ε материала фермы; σ = f (ε ) (ε = λ l ). Зависимости (50) легко обращаются при помощи обратной функции g , выражающей деформации через напряжения ε = g (σ ), и имеют вид λ j = l j g (σ j ) = l j g ( N j Fj ). При помощи этих зависимостей легко решается задача расчета статически определимых ферм из нелинейно упругого материала, но при расчете статически неопределимых ферм здесь остаются еще значительные трудности. ЛЕКЦИЯ № 9 ПЛАН 1.1. Нелинейные дифференциальные уравнения и их системы 1.2. Метод простой итерации Пусть требуется найти простой вещественный корень уравнения f ( x) = 0, (51) где f ( x) – непрерывная функция. Считается, что известен интервал изоляции [a, b] , в котором находится только один корень уравнения (51), и на концах этого интервала функция f ( x) имеет разные знаки. Уравнение (51) представляют в виде x = ϕ ( x), а последовательные приближения к искомому корню находят по итерационной формуле xk +1 = ϕ ( xk ) (k = 1, 2, ...). (52) В качестве первого приближения x1 к искомому корню можно принять любое значение x из отрезка [a, b] . Если в окрестности искомого корня модуль производной ϕ ′( x) удовлетворяют неравенству ϕ ′( x) ≤ K < 1 , итерационный процесс (52) сходится, причем тем быстрее, чем меньше K . ПРИМЕР 1 Найти корень алгебраического уравнения x5 − x − 0, 2 = 0, лежащий в интервале [1;1,1] . Запишем данное уравнение в виде x = ϕ ( x) = 5 x + 0,2. Условие сходимости метода в интервале [1;1,1] выполнено. 1 1 −4 Действительно, ϕ ′( x) = ( x + 0,2) 5 < . 5 5 За начальное приближение примем x1 = 1 . Последовательно по формуле (52) находим: x2 = 5 1 + 0, 2 = 1,0371; x3 = 5 1,0371 + 0,2 = 1,0440; x4 = 5 1,0440 + 0,2 = 1,0446; ..., x7 = 5 1,0448 + 0, 2 = 1,0448. за приближенное значение искомого корня можно принять x = 1,0448 . Рассмотрим теперь метод простой итерации применительно к системе нелинейных алгебраических уравнений. Для нахождения решения системы нелинейных уравнений f i ( x1 , x2 , ..., xn ) = 0 (i = 1, 2, ..., n) (53) методом простой итерации ее записывают в виде xi = ϕ i ( x1 , x2 , ..., xn ) (i = 1, 2, ..., n), где ϕ i – непрерывные функции, имеющие непрерывные первые производные в некоторой области G , содержащей вещественное решение (a1 , a2 , ..., an ) системы уравнений (53). Если решение единственное, а xi(1) – числа, близкие к ai , эти числа принимают за начальное приближение решения, то последовательные приближения к искомому решению находят по итерационной формуле xi( k +1) = ϕ i ( x1( k ) , x2( k ) , ..., xn( k ) ) (i = 1, 2, ..., n) (k = 1, 2, ...). (54) Условия сходимости итерационного процесса (54) n ∑M j =1 ij <1 (i = 1, 2, ..., n), где M ij = max ∂ϕ i . ∂x j Если обозначить X = ( x1 , x2 , ..., xn )Τ , F = (ϕ1 , ϕ 2 , ..., ϕ n )Τ , то формулу (54) можно записать в виде X k +1 = F ( X k ) (k = 1, 2, ...). Рассмотрим решение дифференциальных уравнений методом простой итерации. Пусть нужно найти решение дифференциального уравнения y ′ = f ( x, y ) в некотором промежутке x0 ≤ x ≤ x 0 + a , удовлетворяющее начальному условию y ( x0 ) = y0 . Если в области G , содержащей точку ( x0 , y0 ) , функция непрерывна и имеет непрерывную частную производную f ( x, y ) ∂f , ∂y то рассматриваемая задача имеет единственное решение. Последовательные приближения находят по итерационной формуле x yk +1 ( x) = y0 + ∫ f (t , yk (t ))dt (k = 1, 2, ...). (55) x0 Причем начальное приближение y1 можно выбрать произвольно. ПРИМЕР 2 Найти решение уравнения (3 − y ) y′ + 4 = 0 на отрезке [0,1], удовлетворяющее начальному условию y (0) = 0. Начальное приближение возьмем в виде y1 = 0 и по формуле (55), принимающей вид x 4dt 3 − yk (t ) yk +1 ( x) = − ∫ (k = 1, 2, ...). Последовательно находим 4 4 y2 = − x, y3 = −3ln(1 + x). 3 9 За приближенное решение поставленной задачи можно взять y3 ( x) = −3ln(1 + 4 9 x) , так как y3 (1) = −1,102 , а значение точного решения y3 (1) = −1,123. Рассмотрим теперь решение системы дифференциальных уравнений yi′ = f i ( x, y1 ( x), y2 ( x), ..., yn ( x)) (i = 1, 2, ..., n), удовлетворяющее начальным условиям yi ( x0 ) = yi 0 . Последовательные приближения будем искать по итерационной формуле x yi k +1 ( x) = yi 0 + ∫ f (t , y1 k (t ), y2 k (t ), ..., yn k (t ))dt x0 (i = 1, 2, ..., n), (k = 1, 2, ...). При этом за начальные приближения yi 1 можно взять начальные условия yi 0 . ЛЕКЦИЯ № 10 ПЛАН 1.1. Границы применимости метода упругих решений 1.2. Алгоритм метода упругих решений 1.3. Улучшение сходимости метода Идея метода последовательных приближений – метода простой итерации – легла в основу метода решения физически нелинейных задач механики деформируемого твердого тела, разработанного А. А. Ильюшиным и получившего название метода упругих решений. Для решения задач упругопластического деформирования трехмерного тела была получена система нелинейных дифференциальных уравнений 1 ∂Δ X + = Rx , 1 − 2 μ ∂x G 1 ∂Δ Y ∇ 2V + + = Ry , 1 − 2 μ ∂y G 1 ∂Δ Z ∇ 2W + + = Rz , 1 − 2 μ ∂z G ∇ 2U + (56) левая часть которой представляет собой уравнения, соответствующие решению упругих задач, а правая – содержит члены, учитывающие пластические деформации. Для решения системы (56) методом упругих решений процесс строится так, что на каждом этапе решается упругая задача. В качестве первого приближения принимается w( x, y, z ) = 0 , при этом Rx = Rx(1) = 0 , Ry = Ry(1) = 0 , Rz = Rz(1) = 0 , и отыскивается решение системы (56) , удовлетворяющее условиям закрепления на границе тела. Полученные решения U (1) , V (1) , W (1) используются для отыскания ε i ( x, y, z ) , σ i ( x, y, z ) , ω ( x, y, z ) и вычисления Rx(2) , Ry(2) , Rz(2) . Затем отыскивается второе приближение решения системы (56), при Rx = Rx(2) , Ry = Ry(2) , Ry = Ry(2) . После нахождения U (2) , V (2) , W (2) снова определяются ε i ( x, y, z ) , σ i ( x, y, z ) , ω ( x, y, z ) и вычисляются правые части Rx(3) , Ry(3) , Rz(3) и т. д. Метод упругих решений можно применять и для решения упругопластических задач для пластин и оболочек. Так, если в уравнении ⎛ ∂ 2 Fxx ∂ 2 Fxy ∂ 2 Fyy ⎞ 3 + + = q D1∇ W − D2 ⎜ 2 ⎟ ⎜ ∂x 2 ⎟ 4 ∂ ∂ ∂ x y y ⎝ ⎠ 4 (57) член в скобках перенести вправо, то слева останутся выражения, соответствующие решению упругой задачи ⎛ ∂ 2 Fxx (W ) ∂ 2 Fxy (W ) ∂ 2 Fyy (W ) ⎞ 3 + + D1∇ W = q + D2 ⎜ ⎟⎟ . 2 ⎜ ∂x 2 ∂ ∂ ∂ 4 x y y ⎝ ⎠ 4 Итерационный процесс выражается зависимостью ⎛ ∂ 2 Fxx (Wk ) ∂ 2 Fxy (Wk ) ∂ 2 Fyy (Wk ) ⎞ 3 D1∇ Wk +1 = q + D2 ⎜ + + ⎟⎟ ⎜ ∂x 2 4 x y ∂ ∂ ∂y 2 ⎝ ⎠ 4 (k = 1, 2, ...), где W1 = 0. ЛЕКЦИЯ № 11 ПЛАН 1.1. Метод продолжения решения по параметру 1.2. Сведение нелинейной задачи к последовательному решению линейных задач 1.3. Применение метода продолжения решения по параметру Этот метод позволяет свести нелинейную задачу к последовательному решению линейных задач. Рассмотрим систему m нелинейных уравнений относительно m неизвестных x1 , x2 , ..., xm , содержащих параметр P , f i ( x1 , x2 , ..., xm , P) (i = 1, 2, ..., m). Если ввести вектор-функцию F = ( f1 , f 2 , ..., f m )Τ (58) и вектор X = ( x1 , x2 , ..., xm )Τ , то систему (58) можно записать в виде F ( X , P) = 0. (59) Пусть известно некоторое решение системы (59), т. е. F ( X 0 , P0 ) = 0. Уравнение (59) определяет некоторую неявную функцию X = X ( P ) . Рассмотрим некоторую окрестность A точки ( X 0 , P0 ) . Свойства решений системы (59) в этой окрестности устанавливает теорема о неявных функциях. В ней доказывается, что если: 1) вектор-функция F (т. е. все fi при i = 1, 2, ..., m ) определена и непрерывна в A ; 2) в A существуют и непрерывны частные производные от F по всем аргументам x j ( j = 1, 2, ..., m) и по параметру P ; 3) в точке ( X 0 , P0 ) отличен от нуля якобиан вектор-функции F ( det J ≠ 0 ), то в некоторой окрестности точки ( X 0 , P0 ) решением системы (59) являются однозначные непрерывные функции параметра P xi = xi ( P) (i = 1, 2, ..., m), (60) такие, что xi ( P0 ) = xi 0 (i = 1, 2, ..., m), и производные dxi (i = 1, 2, ..., m) также непрерывны в этой окрестности. dP Якобиан вектор-функции F – это определитель ее матрицы Якоби J ∂F ⎡ ∂fi ⎤ =⎢ ⎥ ∂X ⎢⎣ ∂ j ⎥⎦ (i, j = 1, 2, ..., m). J= Таким образом, теорема о неявных функциях, что при выполнении условий 1–3 решение системы (59) в некоторой окрестности B точки ( X 0 , P0 ) образует единственную кривую K , которая имеет параметрическое представление (60) и проходит через точку ( X 0 , P0 ) . Чтобы получить теперь решение X 1 системы (59) при близком к P0 значение P1 , нужно продвинуться вдоль кривой K . При этом точка ( X 1 , P1 ) должна оставаться внутри окрестности B . Если условия 1–3 выполняются в точке ( X 1 , P1 ) , то решение снова можно продолжить, и т. д. Рассмотрим, как практически перейти от точки ( X 0 , P0 ) к точке ( X 1 , P1 ) . Продифференцируем уравнение (59) по параметру P , учитывая, что X = X ( P ) . Тогда получим систему дифференциальных уравнений ∂F dX ∂F + = 0. ∂X dP ∂P Система (61) линейна относительно (61) dX . Ее решение при условии dP det J ≠ 0 приводит к системе обыкновенных дифференциальных уравнений dX ∂F = J −1 . dP ∂P (62) Так как при P = P0 решение системы (59) известно, то задача решения системы (61) или (62) является задачей Коши при начальном условии X ( P0 ) = X 0 . (63) Для решения начальной задачи (61), (63) можно применять известные методы – Эйлера, Рунге-Кутта и другие. Если для решения начальной задачи (61), (63) применить метод Эйлера, то получим расчетную схему X k +1 = X k + Δ X k , Pk +1 = Pk + Δ Pk , где Δ Pk задается, а Δ X k находится из линейного относительно Δ X k уравнения FX′ ( X k , Pk )ΔX k + FP′ ( X k , Pk )ΔPk = 0. (64) Метод продолжения по параметру может быть применен для решения нелинейных алгебраических и дифференциальных уравнений. Кроме того, идея метода продолжения по параметру послужила развитию целого ряда методов решения сложных задач механики деформируемого твердого тела. Рассмотрим несколько примеров решения нелинейных задач методом продолжения решения по параметру. ПРИМЕР 1 Найти действительное решение нелинейного алгебраического уравнения x5 − x = 0,2. (65) Если правую часть этого уравнения взять за параметр p , то получим F ( x, p ) ≡ x5 − x − p = 0. (66) При p = 0 решение уравнения (66) легко найти (по такому принципу и вводится в уравнение параметр). Это решение будет x = 1 . Необходимо найти решение уравнения (66) при p = 0, 2 . Уравнение (64) для уравнения (66) принимает вид (5 x 4 − 1)Δ xk = Δ pk Причем x1 = 1 . (k = 1, 2, ...). (67) Если Δ pk = 0,2 , то для получения искомого решения следует решить только одно уравнение 4Δ x1 = 0,2. Откуда Δ x1 = 0,05 и x2 = x1 + Δ x1 = 1,05 . Если Δ pk = 0,1 , то для получения искомого решения следует дважды решить уравнение (67) при k = 1 и k = 2 . Имеем 4Δ x1 = 0,1; Δ x1 = 0,025; x2 = x1 + Δ x1 = 1,025; 4,519Δ x2 = 0,1; Δ x2 = 0,022; x3 = x2 + Δ x2 = 1,047. Если Δ pk = 0,05 , то, решив уравнение (67) при k = 1 , k = 2 , k = 3 и k = 4 , получим x5 = 1,0458 . При Δ pk = 0,025 решение уравнения (66) при p = 0, 2 x9 = 1,0453 . Это значение будет x9 и можно принять за приближенное решения уравнения (65). ПРИМЕР 2 Найти действительное решение нелинейного дифференциального уравнения (3 − y ) y′ + 4 = 0, (68) удовлетворяющее начальному условию y (0) = 0 . Нас будет интересовать решение начальной задачи на отрезке [0,1]. Параметр p введем в уравнение (68) так, чтобы при p = 0 легко можно было найти его решение. Итак, получим F ( y, p ) ≡ (3 − y ) y′ + p = 0. (69) Уравнение (64) будет иметь вид (3 − yk )Δ yk′ − yk′ Δ yk + Δ pk = 0 (k = 1, 2, ...). (70) Начальное условие для Δ yk остается тем же, что и для y : Δ yk (0) = 0 . При p = 0 уравнение (69) имеет решение y1 = 0 ( y1′ = 0 ). Если Δ pk = 4 , то уравнение (70) решаем один раз. 4 4 3Δ y1′ + 4 = 0, Δ y1 = − x, y2 = y1 + Δ y1 = − x. 3 3 Если Δ pk = 4 , то уравнение (70) нужно решать при k = 1 и k = 2 : 2 2 3Δ y1′ + 2 = 0, Δ y1 = − x, y2 = − x; 3 3 (3 + 2 2 6x x)Δ y2′ + Δ y2 + 2 = 0, Δ y2 = − , 3 3 9 + 2x 2 6x y3 = − x − . 3 9 + 2x 2 6x . За приближенное решение берем y3 = − x − 3 9 + 2x ЛЕКЦИЯ № 12 ПЛАН 1.1. Определение напряженно-деформированного состояния пологих оболочек 1.2. Сведение к последовательному решению линейных задач 1.3. История нагружения 1.4. Физически нелинейные задачи Нелинейные уравнения равновесия для пологих оболочек ∂ ∂ ∂ ∂ ( Δ1 ) + ( Δ 2 ) = 0, ( Δ 2 ) + ( Δ 3 ) = 0, ∂x ∂y ∂x ∂y ⎡⎛ ∂ 2W ∇ 4W − 12h −2 ⎢⎜ k x + 2 ∂x ⎣⎝ = ⎞ ⎛ ∂ 2W ∂ 2W 2 k Δ + Δ + + 2 ⎟ 1 ⎜ y ∂x∂y ∂y 2 ⎠ ⎝ 2 (1 − μ 2 ) q Eh3 ⎞ ⎤ ⎟ Δ3 ⎥ = ⎠ ⎦ (71) , 2 ⎛ ∂V ∂U 1 ⎛ ∂W ⎞ ⎜ Δ1 = − k xW + ⎜ + − k yW + μ ⎟ ⎜ ∂y ∂x 2 ⎝ ∂x ⎠ ⎝ ⎛ ∂U ∂V ∂W ∂W ⎞ Δ 2 = μ1 ⎜ + + ⎟, ∂ ∂ ∂ ∂ y x x y ⎝ ⎠ 2 1 ⎛ ∂W ⎞ ⎞ ⎜ ⎟ ⎟, 2 ⎝ ∂y ⎠ ⎠⎟ 2 2 ⎛ ∂U 1 ⎛ ∂W ⎞ ⎞ ∂V 1 ⎛ ∂W ⎞ Δ3 = μ ⎜ − k xW + ⎜ + − k yW + ⎜ ⎟ ⎟ ⎟ , ⎜ ∂x ⎟ ∂y ∂ ∂ 2 2 x y ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ∂4 ∂4 ∂4 ∇ 4 = 4 + 2 2 2 + 4 , μ1 = 0,5(1 − μ ), ∂x ∂x ∂y ∂y и уравнения в смешанной форме ⎛ ∂ 2W ⎞ ∂ 2 Φ ⎛ ∂ 2W ⎞ ∂ 2 Φ ∂ 2W ∂ 2 Φ − ⎜ kx + 2 ⎟ 2 + ⎜ k y + 2 ⎟ 2 − 2 ∂x ⎠ ∂y ∂y ⎠ ∂x ∂x∂y ∂x∂y ⎝ ⎝ Eh 2 q 4 − ∇ + = 0, W 12(1 − μ 2 ) h 2 ⎛ ∂ 2W ⎞ ∂ 2W ∂ 2W 1 4 ∂ 2W ∂ 2W ∇ Φ=⎜ − kx 2 − k y 2 ⎟ − 2 2 E x y x y ∂ ∂ ∂ ∂ ∂y ∂x ⎝ ⎠ кратко можно записать в виде (72) F (r , P ) = 0, (73) где r = (U , V , W )Τ , если рассматриваются уравнения (71) и r = (W , Φ )Τ , если рассматриваются уравнения (72); P = (0, 0, q)Τ – параметр нагрузки. Причем известно, что в ненагруженном состоянии перемещения и напряжения равны нулю, т. е. при P = 0 , r = 0 , так что r = r ( P) . Если за параметр взята нагрузка, то расчетная схема метода продолжения решения по параметру (64), соответствующая методу последовательных нагружений, имеет вид rk +1 = rk + Δ rk , Pk +1 = Pk + Δ Pk , где Δ Pk задается, а Δ rk находится из линейного относительно Δ rk уравнения Fr′(rk , Pk )Δ rk + FP′ (rk , Pk )Δ Pk = 0. (74) Причем r1 = 0, P1 = 0. Этот метод был разработан В. В. Петровым в 1959 году и получил большое применение для решения нелинейных задач для пластин и оболочек, так как нелинейную задачу (73) этот метод сводит к последовательному решению линейных задач (74). На каждом этапе нагружения коэффициенты уравнения (74) зависят от «истории» нагружения. Рассмотрим уравнения в смешанной форме (72) для пластин и пологих оболочек, допускающих прогибы, соизмеримые с толщиной. Для этого обозначим ΔWk = wk , ΔΦ k = ϕk , ΔPk = pk . Причем k −1 k −1 k −1 i =1 i =1 i =1 Wk = ∑ wi , Φ k = ∑ϕi , Pk = ∑ pi . Теперь уравнения (74) примут вид (индекс k опущен) ⎛ ∂ 2W ⎞ ∂ 2ϕ ∂ 2 Φ ∂ 2 w ⎛ ∂ 2W ⎞ ∂ 2ϕ ∂ 2 Φ ∂ 2 w + ⎜ ky + 2 ⎟ 2 + 2 − ⎜ kx + 2 ⎟ 2 + 2 ∂x ⎠ ∂y ∂y ∂x 2 ⎝ ∂y ⎠ ∂x ∂x ∂y 2 ⎝ ⎛ ∂ 2W ∂ 2ϕ ∂2Φ ∂2 w ⎞ Eh 2 p + ∇ 4W + = 0, 2⎜ ⎟− 2 h ⎝ ∂x∂y ∂x∂y ∂x∂y ∂x∂y ⎠ 12(1 − μ ) (75) ∂ 2W ∂ 2 w ∂ 2W ∂ 2 w ∂ 2W ∂ 2 w ∂2 w ∂2 w 1 4 ∇ ϕ =2 − − 2 − kx 2 − k y 2 . ∂x∂y ∂x∂y ∂x 2 ∂y 2 ∂y ∂x 2 ∂y ∂x E Линейные системы дифференциальных уравнений в частных производных (75) при заданных краевых условиях можно решить методом Бубнова-Галеркина, МКЭ и другими методами решения линейных краевых задач. Методом последовательных нагружений с успехом можно применять и для решения физически нелинейных задач. Рассмотрим уравнение ⎛ ∂ 2 Fxx ∂ 2 Fxy ∂ 2 Fyy D1∇ W − D2 ⎜ + + ⎜ ∂x 2 x y ∂ ∂ ∂y 2 ⎝ 2 ⎞ 3 ⎟⎟ = q, ⎠ 4 (76) где Fxx , Fxy , Fyy – нелинейные дифференциальные выражения функции W . 3 3 2 3 3 2 ⎛ ∂ 2W ⎞ 3 ∂ 2W ∂ 2W 2 1 ⎛ ∂ 2W ⎞ ⎛ ∂ 2W ⎞ 2 Fxx = ⎜ 2 ⎟ + ∇ W + ⎜ 2 ⎟ +⎜ ⎟ ∇1 W , 2 2 x 2 x y 2 y x y ∂ ∂ ∂ ∂ ∂ ∂ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎛ ∂ 2W ⎞ 3 ∂ 2W ∂ 2W 2 1 ⎛ ∂ 2W ⎞ ⎛ ∂ 2W ⎞ 2 Fyy = ⎜ 2 ⎟ + ∇ W + ⎜ 2 ⎟ +⎜ ⎟ ∇ 2W , 2 2 2 ⎝ ∂x ⎠ ⎝ ∂x∂y ⎠ ⎝ ∂y ⎠ 2 ∂x ∂y 3 2 2 ⎛ ∂ 2W ⎞ ∂ 2W ⎡⎛ ∂ 2W ⎞ ⎛ ∂ 2W ⎞ ⎤ ∂ 2W ∂ 2W ∂ 2W , Fxy = ⎜ ⎢⎜ 2 ⎟ + ⎜ 2 ⎟ ⎥ + 2 ⎟ + 2 x y x y x y x y x y ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎢⎣⎝ ⎝ ⎠ ⎠ ⎝ ⎠ ⎥⎦ 1 ∂2 1 ∂2 ∂2 ∂2 ∂2 1 ∂2 2 2 , ∇2 = . ∇ = 2 + 2 , ∇1 = 2 + + 2 ∂x 2 2 ∂x 2 2 ∂x 2 ∂x ∂y ∂x 2 Обозначим Δ Wk = wk , Δ Pk = pk , тогда линеаризованные уравнения (74) для уравнения (76) будут иметь вид (индекс k опущен) D1∇ 2 w − D2 L(W , w) = 3 p, 4 (77) где L(W , w) – линеаризованный дифференциальный оператор для стоящего в скобках члена в (76), содержащий W и w (ввиду его громоздкости здесь не приводится). Для решения уравнения (77) при заданных краевых условиях для W (и соответственно для w ) можно применять метод Бубнова-Галеркина. ЛЕКЦИЯ № 13 ПЛАН 1.1. Уравнения в конечных разностях 1.2. Разностная схема для уравнений в частных производных 1.3. Аппроксимация точной сеточной задачи 1.4. Выбор решетки и сходимость сеточных решений Рассмотрим этот метод на примере решения уравнения 2-го порядка y ′′ + p( x) y ′ + q( x) y = f ( x) (78) при следующих краевых условиях: α 0 y (a ) + α1 y′(a ) = A, β 0 y (b) + β1 y′(b) = B. (79) Нужно найти функцию y ( x) , удовлетворяющую на отрезке [a, b] дифференциальному уравнению (78), а на концах отрезка – краевым условиям (79). Метод конечных разностей позволяет свести краевую задачу к решению системы алгебраических уравнений. Разобьем отрезок [a, b] на n равных частей точками x1 , x2 , ..., xn−1 (рис. 5) с шагом h = xi = a + ih (i = 0,1, ..., n) . y y0 y1 y2 a = x0 x1 x2 yi xi yn b = xn Рис. 5 Введем обозначения: y ( xi ) = yi , y ′( xi ) = yi′, y ′′( xi ) = yi′′, p ( xi ) = pi , q ( xi ) = qi , f ( xi ) = f i . x b−a . Тогда n Считается, что на отрезке [a, b] функции p ( x) , q ( x) , f ( x) непрерывны; α 0 , α1 , β 0 , β1 – заданные числа, причем α 0 + α1 ≠ 0, β 0 + β1 ≠ 0. Идея метода конечных разностей заключается в замене производных конечными разностями, в результате чего дифференциальное уравнение превращается в систему алгебраических уравнений относительно значений искомой функции y ( x) в точках xi . Для представления производных через конечные разности нужно использовать разложение функции в ряд Тейлора y ( xi + h) = y ( xi ) + y ′( xi ) y ′′( xi ) 2 h+ h + ... 1! 2! Если отбросить третий член ряда, то получим y ( xi + h) = y ( xi + h) − y ( xi ) + O(h 2 ). h O (h 2 ) говорит о порядке погрешности. Аналогично выразим y ′′( xi ) : y ′′( xi ) = y ( xi + 2h) − y ( xi + h) − ( y ( xi + h) − y ( xi ) ) h 2 + O(h 2 ). Итак, используя формулу Тейлора, получим следующие формулы: yi +1 − yi ; h y − 2 yi +1 + yi . yi′′ ≈ i + 2 h2 yi′ ≈ (80) Это формулы первого порядка точности (O(h 2 )) . Можно получить и более точные формулы: yi +1 − yi −1 ; 2h y − 2 yi + yi −1 yi′′ ≈ i +1 . h2 yi′ ≈ Это формулы второго порядка точности (O(h3 )) . (81) Для граничных точек имеем y0′ ≈ y1 − y1 ; h yn′ ≈ yn − yn−1 . h (82) Так как в дальнейшем мы будем пользоваться приближенными равенствами (80) – (82), то и решение соответственно будет приближенным. Заменяя производные в точке xi в уравнении (78) приближенными равенствами (81), мы сведем дифференциальное уравнение к системе алгебраических уравнений. Понятно, что точность замены будет тем выше, чем меньше выбран шаг h . yi +1 − 2 yi + yi −1 y − yi −1 + pi i +1 + qi yi = f i , 2 2h h (i = 1, 2, ..., n − 1). (83) Мы получили систему (n − 1) уравнений относительно неизвестных y0 , y1 , ..., yn . Недостающие два уравнения получим из краевых условий (79): α 0 y0 + α1 y1 − y0 y − yn−1 = A, β 0 yn + β1 n = B. h h (84) Преобразуем уравнения (83) pi ⎞ pi ⎞ ⎛ 1 ⎛ 1 ⎛ 2 ⎞ ⎜ h 2 − 2h ⎟ yi −1 + ⎜ − h 2 + qi ⎟ yi + ⎜ h 2 + 2h ⎟ yi +1 = fi (i = 1, 2, ... n − 1). ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ После введения обозначений Ai = 1 − pi h ph , Ci = 2 − qi h 2 , Bi = 1 + i , Fi = f i h 2 2 2 (85) получится простая система Ai yi −1 − Ci yi + Bi yi +1 = Fi (i = 1, 2, ... n − 1). (86) Это – система с трехдиагональной матрицей, то есть в каждом уравнении всего три неизвестных. Следовательно, для ее решения можно использовать метод прогонки. ПРИМЕР Используя метод конечных разностей, найти решение дифференциального уравнения y ′′ + 2 x 2 y ′ − xy = 2 при краевых условиях y (1) = 2, y (2) = 3. Примем h = 1 , тогда x0 = 1 ; x1 = 1, 25 ; x2 = 1,5 ; x3 = 1,75 ; x4 = 2 . Здесь 4 pi = 2 xi2 , qi = − xi , fi = 2 . Система линейных алгебраических уравнений будет иметь вид (86) Ai yi −1 − Ci yi + Bi yi +1 = Fi , (i = 1, 2, 3), где Ai = 1 − xi2 x2 x 1 , Bi = 1 + i , Ci = 2 + i , Fi = . 4 4 16 8 Краевые условия дадут еще два уравнения: y0 = 2, y4 = 3. Процесс вычисления коэффициентов Ai , Bi , Ci сведем в табл. 1. Таблица 1 i xi xi2 Ai Ci Bi Fi 1 2 3 1,250 1,500 1,750 1,562 2,250 3,063 0,610 0,437 0,234 2,078 2,094 2,109 1,390 1,563 1,766 0,125 0,125 0,125 Вычислим прогоночные коэффициенты по формулам αk = Bk Aβ − f , β k = k k −1 k , Ck − Ak α k −1 Ck − Ak α k −1 учитывая, что α 0 = 0, β 0 = 2, α 4 = 0, β 4 = 3, α1 = 0,669; β1 = 0,527; α 2 = 0,868; β 2 = 0,058; α3 = 0,927; β3 = −0,058. Так как y4 известно, то по формуле xk = α k xk +1 + β k последовательно находим: y3 = 2,723; y2 = 2,419; y1 = 2,144 . Разностная схема для уравнений в частных производных Метод конечных разностей для уравнений в частных производных называют методом сеток. Рассмотрим решение дифференциального уравнения ⎛ ∂ 2u ∂ 2u ⎞ − ⎜ 2 + 2 ⎟ = f ( x, y ), ∂y ⎠ ⎝ ∂x (87) или кратко −Δ u = f ( x, y ). Итак, нужно найти решение этого уравнения – функцию u ( x, y ) в некоторой области D , которая бы на границе этой области Γ удовлетворяла краевому условию u Γ = 0. (88) Считаем, что функция f ( x, y ) в области D является достаточно гладкой функцией (чтобы задача имела решение). Для простоты область D возьмем в виде D : {0 ≤ x ≤ 1; 0 ≤ y ≤ 1}. Для решения краевой задачи (87), (88) применяем метод сеток. При этом производные заменяем конечными разностями и сводим задачу к решению системы линейных алгебраических уравнений. Разобьем область D линиями, параллельными осям координат, на части: по оси x – на n частей, по оси y – на m часей. Тогда 1 n 1 y j = jτ, τ = m xi = ih, h = (i = 0,1, ..., n), ( j = 0,1, ..., m). Введем обозначения: u ( xi , y j ) = uij ; f ( xi , y j ) = fi j . Для замены производных конечными разностями воспользуемся формулами (81), т. е. ∂ 2u ∂x 2 ∂u ∂y 2 x = xi y= y j 2 x = xi y= y j uij+1 − 2uij + uij−1 ≈ , h2 ≈ ui j +1 (89) j −1 − 2ui + ui . τ2 j Подставив (89) в выражение (87), получим ⎛ uij+1 − 2uij + uij−1 uij +1 − 2uij + uij −1 ⎞ j −⎜ + ⎟ = fi 2 2 h τ ⎝ ⎠ (90) (i = 1, 2, ..., n − 1), ( j = 1, 2, ..., m − 1). Распишем краевые условия, т. к. в (90) не хватает 4 уравнений. Всего должно быть (n + 1)(m + 1) уравнений. На границе u0j = 0, unj = 0 ( j = 0,1, ..., m), (91) ui0 = 0, uim = 0 (i = 0,1, ..., n). Систему (90), (91) можно решить методом Гаусса. Точка ( xi , y j ) называется узлом. Узлы разделяют на внутренние и граничные. Обозначим через W внутренние узлы { } W : ( xi , y j ) ; i = 1, 2, ..., n − 1; j = 1, 2, ..., m − 1 , а через W ′ – граничные узлы сетки { } W ′ : ( xi , y j ) ; i = 0, n; j = 0, m . Выражение (90) называют разностной схемой краевой задачи (87) и кратко записывают так: Auij W = fi j , (92) где A – некоторый разностный оператор. Краевое условие (88) можно кратко записать: l uij W′ = 0. (93) Разностная схема (92), (93) аппроксимирует краевую задачу (87), (88). Для решения разностной схемы (92), (93) разработан специальный прием, позволяющий свести задачу к последовательному решению систем с трехдиагональной матрицей (чтобы воспользоваться методом прогонки), который получил название метода переменных направлений. ЛЕКЦИЯ № 14 ПЛАН 1.1. Явная и неявная схемы 1.2. Устойчивость разностной схемы Рассмотрим простейший случай, когда область по пространственным координатам одномерная. Требуется найти непрерывную в замкнутой области D функцию u ( x, t ) , которая внутри области удовлетворяет уравнению теплопроводности ∂u ∂ 2u − = f ( x, t ), ∂t ∂x 2 (94) при t = 0 начальному условию u ( x,0) = s ( x) , при x = 0 , x = 1 удовлетворяет краевым условиям u (0, t ) = p (t ), u (1, t ) = q (t ). (95) Функции s ( x) , p (t ) , q (t ) считаются гладкими функциями, причем s (0) = p (0), s (1) = q (0). Область D возьмем прямоугольную D : {0 ≤ x ≤ 1; 0 ≤ T ≤ T1}. Разобьем область D на части линиями, параллельными осям координат: по x шаг возьмем равным h = 1 , n по t шаг возьмем равным τ = 1 . m Тогда xi = ih (i = 0,1, ..., n), t j = jτ ( j = 0,1, ..., m). Введем обозначения: u ( xi , t j ) = uij ; f ( xi , t j ) = fi j . Рассмотрим сетки { } Wh : ( xi , t j ) ; i = 0,1, ..., n; j = 0,1, ..., m , (96) т. е. Wh – множество внутренних и граничных точек. { } Wh′ : ( xi , t j ) ; i = 1, 2, ..., n − 1; j = 1, 2, ..., m − 1 , (97) т. е. Wh′ – множество внутренних точек. W * : {u0j , unj , ui0 , uim } , т. е. W * – сетка, состоящая из граничных узлов. Производные аппроксимируем выражениями uij − uij −1 ⎛ ∂u ⎞ , = ⎜ ⎟ τ t ∂ ⎝ ⎠(i , j ) ⎛ ∂ 2u ⎞ uij+1 − 2uij + uij−1 . ⎜ 2⎟ = h2 ⎝ ∂x ⎠(i , j ) Введем разностный оператор uij+1 − 2uij + uij−1 Aui = − . h2 j На сетке W * зададим оператор luij W* = g. (В сеточную функцию g входят начальные и краевые условия.) Рассмотрим две разностные схемы для уравнения (94): uij − uij −1 + Auij −1 = f i j −1 , L1ui ≡ τ lu W * = g ; j uij − uij −1 + Auij = f i j , τ = g. L2uij ≡ lu W * (98) (99) Обе разностные схемы (98) и (99) называются двухслойными, т. к. рассматривается слой j − 1 и слой j . Разрешив разностную схему (98) относительно uij , получим uij = uij −1 + τ j −1 u − 2uij −1 + uij+−11 ) + τfi j −1. 2 ( i −1 h (100) Поскольку ui0 , u0j , unj (при i = 1, 2, ..., n − 1; j = 0,1, ..., m ) известны, решения разностной схемы (98) находятся по формуле (100) явно, слой за слоем. Поэтому разностная схема (98) называется явной. Разностное уравнение (99) можно представить в виде: τ j ⎛ 2τ ⎞ j τ j u − ⎜1 + 2 ⎟ ui + 2 ui +1 = −uij −1 − τfi j . 2 i −1 h h ⎝ h ⎠ (101) Неизвестных в каждом уравнении три. Это уже другая схема. Исходя из краевых и начальных условий u0j = p, unj = q, разностная схема (99) представляет собой неявную разностную схему, решение которой может быть найдено с помощью уравнений (101) методом прогонки. Разностная схема (98) (явная схема) проста для решения задач, но она устойчива лишь при соотношении τ и h : τ 1 ≤ . h2 2 Неявная разностная схема (99) устойчива при любом соотношении шагов τ и h . ЛЕКЦИЯ № 15 ПЛАН 1.1. Идеи метода конечных элементов 1.2. Вариационная постановка задачи Метод конечных элементов (МКЭ) получил большое применение в различных областях техники при расчете конструкций и целых сооружений и аппаратов. Его можно рассматривать как один из вариантов вариационных методов. Если в методе Ритца и Бубнова-Галеркина аппроксимирующие функции задаются на всей области, занимаемой конструкцией, то в МКЭ эти функции задаются свои на каждом элементе, на которые разбивается вся область. При увеличении конечных элементов приходится изменять все аппроксимирующие функции. Кроме того по сравнению с методом Ритца или Бубнова-Галеркина МКЭ требует гораздо больше вычислительных затрат при формировании системы алгебраических уравнений. Поэтому его целесообразно применять для расчетов сложных объектов: целых сооружений, элементов конструкций сложного очертания и не стандартного закрепления на границе, составных конструкций. Рассмотрим этот метод на примере определения прогиба жестких прямоугольных плит (пластин). X x a B ek A b C D y Y Рис. 6. Пусть прямоугольная пластинка (рис. 6) испытывает изгиб под действием произвольной поперечной нагрузки q ( x, y ) . Разобьем пластинку на n прямоугольных элементов ek со сторонами a и b . Связь конечных элементов между собой осуществляется в узлах. В каждом узле задаем по три ∂w ∂w и ). Потребуем ∂x ∂y перемещения (прогиб w и два угла поворота совместности вертикальных перемещений и углов поворота относительно местных осей x и y в узловых точках для примыкающих к узлу конечных элементов. Обобщенные перемещения в узлах конечного элемента обозначим через qi . Прямоугольный конечный элемент должен иметь, как это показано на рис. 7, двенадцать обобщенных перемещений, или двенадцать степеней свободы [11]. x q6 (a,0) q4 q3 (0,0) q9 ( a, b) q5 q7 q12 (0, b) q2 q1 q8 q10 y q11 z, w Рис. 7. Перемещения между узловыми точками на конечном элементе ek зададим многочленом с таким числом произвольных коэффициентов, какое число степеней свободы имеет конечный элемент (в рассматриваемом примере 12). На рис. 3 показаны положительные направления векторов перемещений узловых точек q1 - q12 . При принятом направлении осей координат положительным направлениям углов поворота q2 , q5 , q8 , q11 соответствуют отрицательные производные ∂w , ∂x а положительным направления углов поворота q3 , q6 , q9 , q12 соответствуют положительные производные ∂w . ∂y Если, например, элемент пластинки (0,0) прогибается в положительном направлении (вниз), то при этом касательная поворачивается по часовой стрелке, т. е. q3 = время положительному значению от узла ∂w >0 и ∂y ∂w (0,0) . В то же ∂y ∂w будет соответствовать вращение ∂x касательной в направлении, обратном q2 , т. е. q2 = − ∂w (0,0) . ∂x Таким образом, для узловых перемещений можно записать следующие выражения: ∂w ∂w (0,0); q3 = (0,0); q4 = w(a,0); ∂x ∂y ∂w ∂w ∂w (a,0); q7 = w(a, b); q8 = − (a, b); q5 = − (a,0); q6 = ∂x ∂y ∂x ∂w ∂w ∂w q9 = (a, b); q10 = w(0, b); q11 = − (0, b); q12 = (0, b). ∂y ∂x ∂x q1 = w(0,0); q2 = − (102) Представим функцию wkn ( x, y ) на конечном элементе ek в виде суммы многочленов Эрмита (эти многочлены третьей степени благодаря своим свойствам очень удобны для МКЭ) 2 2 wkn = ∑∑ aij Э0i ( x)Э0 j ( y ) + bij Э1i ( x)Э0 j ( y ) + cij Э0i ( x)Э1 j ( y ), i =1 j =1 где многочлены Эрмита Э0i ( x) , Э1i ( x) имеют вид a 3 − 3ax 2 + 2 x3 3ax 2 − 2 x3 Э01 ( x) = ; Э02 = ; a3 a3 − ax 2 + x3 a 2 x − 2ax 2 + x 3 Э11 ( x) = ; Э12 = . a2 a2 Заменой a на b и x на y получаются многочлены Э0 j ( y ) , Э1 j ( y ) . Многочлены Эрмита обладают следующим свойством: (103) Э01 (0) = 1; Э02 (0) = Э11 (0) = Э12 (0) = 0; ∂Э01 ∂Э ∂Э ∂Э (0) = 02 (0) = 12 (0) = 0; 11 (0) = 1; ∂x ∂x ∂x ∂x Э01 (a ) = Э11 (a) = Э12 (a ) = 0; Э02 (a ) = 1; ∂Э01 ∂Э ∂Э ∂Э (a ) = 02 (a ) = 11 (0) = 0; 12 (a ) = 1. ∂x ∂x ∂x ∂x Аналогичные выражения имеют место для Э0 j ( y ) , Э1 j ( y ) и их производных по y при y = 0 и y = b . С учетом этих свойств легко выразить коэффициенты разложения функции wkn ( x, y ) в ряд (103) через обобщенные перемещения qi на конечном элементе ek w(0,0) = q1 = a11; − ∂w ∂w (0,0) = q2 = −b11; (0,0) = q3 = c11; ∂x ∂y w(a,0) = q4 = a21; − ∂w ∂w (a,0) = q5 = −b21; (a,0) = q6 = c21; ∂x ∂y w(a, b) = q7 = a22 ; − ∂w ∂w (a, b) = q8 = −b22 ; ( a, b) = q9 = c22 ; ∂x ∂y w(0, b) = q10 = a12 ; − ∂w ∂w (0, b) = q11 = −b12 ; (0, b) = q12 = c12 . ∂x ∂y Если теперь ввести обозначения ϕ1 = Э01 ( x)Э01 ( y ); ϕ2 = −Э11 ( x)Э01 ( y ); ϕ3 = Э01 ( x)Э11 ( y ); ϕ4 = Э02 ( x)Э01 ( y ); ϕ5 = − Э12 ( x)Э01 ( y ); ϕ6 = Э02 ( x)Э11 ( y ); ϕ7 = Э02 ( x)Э02 ( y ); ϕ8 = − Э12 ( x)Э02 ( y ); ϕ9 = Э02 ( x)Э12 ( y ); ϕ10 = Э01 ( x)Э02 ( y ); ϕ11 = − Э11 ( x)Э02 ( y ); ϕ12 = Э01 ( x)Э12 ( y ), то функцию wkn ( x, y ) на конечном элементе ek можно представить в виде 12 wkn = ∑ qiϕi , i =1 ( x, y ) ∈ ek , а на остальных элементах она берется равной нулю. (104) Полученное выражение для функции прогиба обеспечивает непрерывность прогибов w и их первых производных по x и y между узлами по линии контакта конечных элементов. Так как аппроксимация wkn ( x, y ) на каждом элементе ek (104) записана в местной системе координат x0 y , то необходимо осуществить переход к общей системе координат X 0Y (в данном случае нужно осуществить только параллельный перенос осей координат, но в некоторых случаях может понадобиться и поворот осей координат). На всей области, занимаемой пластинкой, функция прогиба w( X , Y ) представляется в виде (область разбита на n конечных элементов) n wn ( X , Y ) = ∑ wkn ( X , Y ). k =1 Теперь к функционалу полной энергии деформации жесткой пластинки ⎡⎛ ∂ 2 w ⎞ 2 ∂ 2 w ∂ 2 w ⎤ D ⎧⎪ 2 2 q ⎫⎪ n n n 2 − − Э = ∑ ∫ ⎨(∇ wn ) + 2(1 − μ ) ⎢⎜ wn ⎬ dXdY ⎥ ⎟ 2 2 2 X Y X Y D ∂ ∂ ∂ ∂ k =1 ⎢⎣⎝ ⎥⎦ ⎠ Vk ⎪ ⎪⎭ ⎩ n ( Vk – область, занимаемая конечным элементом ek ) применим метод Ритца и получим систему алгебраических уравнений (в данном случае линейную) относительно неизвестных узловых перемещений каждого конечного элемента. Эта система будет разрежена (содержит много нулей), так как функция wkn на всех элементах кроме ek равна нулю. Покажем формирование системы алгебраических уравнений на примере одного конечного элемента ek . В данном случае метод Ритца применяем к функционалу a b ⎫⎪ ⎡⎛ ∂ 2 w ⎞ 2 ∂ 2 w ∂ 2 w ⎤ D ⎧⎪ 2 q 2 kn kn kn ⎥ − 2 wkn ⎬ dxdy. Э = ∫∫ ⎨(∇ wkn ) + 2(1 − μ ) ⎢⎜ ⎟ − 2 2 x y x y D 2 0 0⎪ ∂ ∂ ∂ ∂ ⎢ ⎥⎦ ⎠ ⎪⎭ ⎣⎝ ⎩ (105) Заготовим выражения подынтегральной функции через узловые перемещения (считаем q ( x, y ) = q0 ) 12 12 12 ∇ 2 wkn = ∑ qi ∇ 2ϕi ; (∇ 2 wkn ) 2 = ∑∑ qi q j ∇ 2ϕi ∇ 2ϕ j ; i =1 i =1 j =1 2 2 12 12 ∂ 2 wkn 12 ∂ 2ϕi ⎛ ∂ 2 wkn ⎞ ∂ 2ϕi ∂ ϕ j = ∑ qi ;⎜ ⋅ ; ⎟ = ∑∑ qi q j ∂x∂y i =1 ∂x∂y ⎝ ∂x∂y ⎠ ∂x∂y ∂x∂y i =1 j =1 2 12 ∂ 2 wkn ∂ 2 wkn 12 12 ∂ 2ϕi ∂ ϕ j ⋅ = ∑∑ qi q j 2 ⋅ 2 ; qwkn = ∑ q0 qiϕi . ∂x 2 ∂y 2 ∂x ∂y i =1 j =1 i =0 Теперь функционал (105) принимает вид 12 1 12 12 Э = ∑∑ qi q j kij − ∑ qi Pi , 2 i =1 j =1 i =1 (106) ⎡ 2 ⎛ ∂ 2ϕi ∂ 2ϕ j ∂ 2ϕi ∂ 2ϕ j ⎞ ⎤ 2 kij = D ∫∫ ⎢∇ ϕi ∇ ϕ j + 2(1 − μ ) ⎜ dxdy; ⋅ − ⋅ ⎜ ∂x∂y ∂x∂y ∂x 2 ∂y 2 ⎟⎟ ⎥ 0 0⎢ ⎝ ⎠ ⎥⎦ ⎣ (107) где обозначено a b a b Pi = q0 ∫∫ ϕi dxdy. (108) 0 0 Найдем ∂Э ( p = 1, 2, ...,12) и приравняем их к нулю, в результате ∂q p получим систему алгебраических уравнений относительно q1 , q2 , ..., q12 ( p заменим на i и считаем матрицу коэффициентов kij симметричной) 12 ∑q k j =1 j ij = Pi (i = 1, 2, ...,12). (109) ЛЕКЦИЯ № 16 ПЛАН 1.1. Решение дифференциальных уравнений методом конечных элементов 1.2. Типы конечных элементов 1.3. Расчет НДС пластинки В литературе, относящейся к МКЭ, вводятся понятия узловых усилий Ri , которые имеют вид (на примере одного конечного элемента ek ) 12 Ri = ∑ q j kij , (110) j =1 коэффициентов матрицы жесткости kij и узловых эквивалентных сил Pi (внешних сил), а для построения системы (109) пользуются принципом возможных перемещений, который гласит, что сумма работ узловых усилий на возможных приращениях перемещений δ qi (работа внутренних сил в узловых точках) 12 ∑Rδq i =1 i i равна сумме работ узловых эквивалентных сил Pi на тех же перемещениях δ qi (работа внешних сил в узловых точках) 12 ∑ Pδ q , i i =1 i т. е. 12 12 ∑ R δ q = ∑ Pδ q . i =1 i i i =1 i Так как δ qi произвольны, то получим Ri = Pi , что и есть система (109). В матричной форме R = K ⋅ q, i где R – вектор-столбец узловых усилий; K – матрица жесткости конечного элемента пластины, размерность которой 12 × 12 ; q – векторстолбец узловых перемещений. Кроме того, P – вектор-столбец узловых эквивалентных сил. Следовательно, систему (109) можно записать в матричной форме так: K ⋅ q = P, а ее решение q = K −1 ⋅ P. Система (109) – это система для четырех узлов одного конечного элемента (12 неизвестных по три в каждом из четырех узлов). Аналогичную систему составляют для каждого конечного элемента. Для граничных узлов конечного элемента используются краевые условия. Выражения для внутренних усилий Ri (110) и узловых сил Pi (108) записаны в местной системе координат, относящейся к выделенному конечному элементу. Если для разных конечных элементов выражения для Ri и Pi записаны в разных местных системах координат, то для определения матрицы внутренних усилий всей конструкции необходимо построить матрицу жесткости всей конструкции, а для этого необходимо обеспечить переход от местной к общей системе координат. На рис. 8 представлены простейшие типы конечных элементов, используемые в МКЭ: x y x z y x x z y y x z x y z Рис. 8. Простейшие типы конечных элементов 1. Элемент стержня, работающий в условиях напряженного состояния на растяжение-сжатие (рис 4, а). Число степеней свободы конечного элемента 2. 2. Элемент стержня, работающий на растяжение-сжатие в осевом направлении, на скручивание и изгиб в двух взаимно перпендикулярных плоскостях (рис 4, б). Число степеней свободы конечного элемента 12. 3. Плоский треугольный элемент пластины (рис. 4, в). Число степеней конечного элемента 6. 4. Прямоугольный плоский элемент пластины (рис. 4, г). Число степеней конечного элемента 8. 5. Тетраэдр (рис. 4, д). Конечный элемент такой формы используется при решении объемной задачи теории упругости. Число степеней конечного элемента 12. 6. Плоский прямоугольный элемент в задачах изгиба жестких пластин (рис. 4, е). Число степеней конечного элемента 12. 7. Треугольный элемент произвольной оболочки (рис. 4, ж). Число степеней конечного элемента 18. Расчет НДС пластинки МКЭ Найдем МКЭ прогиб квадратной в плане со стороной a жесткой пластинки, защемленной по контуру и находящейся под действием распределенной поперечной нагрузки q . Разобьем пластинку на четыре равные части, т. е. на четыре конечных элемента (рис. 5). x a a B 2 A C D a 2 a y Рис. 5 В силу симметрии можем рассматривать только один четырехугольный элемент ABCD . Обозначим перемещения в точке A через q1 , q2 , q3 , в точке B – q4 , q5 , q6 , в точке C – q7 , q8 , q9 , в точке D – q10 , q11 , q12 . Учитывая граничные условия жесткого закрепления (первые производные от прогиба равны нулю) и симметричность деформации пластинки относительно центральных осей, заключаем, что из двенадцати перемещений только q7 не равно нулю. В результате получим из (109) только одно уравнение q7 ⋅ k77 = P7 . Используя (107) и (108), найдем (в выражениях многочленов Эрмита нужно a и b заменить на a 2 ) a 2a 2 k7,7 = D ∫ ∫ ⎧⎪ ⎡⎛ ∂ 2ϕ ⎞2 ∂ 2ϕ ∂ 2ϕ ⎤ ⎫⎪ 79,44 D 2 2 7 7 ⋅ 27 ⎥ ⎬dxdy = ; ⎨(∇ ϕ7 ) + 2(1 − μ ) ⎢⎜ ⎟ − 2 2 x y x y a ∂ ∂ ∂ ∂ ⎢ ⎥ ⎝ ⎠ ⎪⎩ ⎣ ⎦ ⎪⎭ a 2a 2 P7 = q ∫ ∫ ϕ dxdy = 0,0625qa . 2 7 Следовательно, qa 4 q7 = 0,000787 . D Полученное значение прогиба в центре пластинки при разбиении ее на 2 × 2 конечных элемента существенно отличается от точного решения wmax qa 4 = 0,00126 . Для получения более точного решения необходимо D пластинку разбить на большее число элементов. В работе Постнова В. А. и Хархурима И. Я. [11] приводятся значения прогиба в центре квадратной со стороной a пластины, шарнирно неподвижно закрепленной по контуру и нагруженной равномерно распределенной поперечной нагрузкой q (табл. 2). Таблица 2 Размер сетки 2× 2 4× 4 8×8 16 × 16 Точное решение Перемещение wmax qa 4 D 0,003444 0,003945 0,004040 0,004058 0,004062 Из вышеизложенного видно, что МКЭ не имеет смысла применять для таких простых конструкций как прямоугольные пластинки (и оболочки), т. к. для них первое приближение метода Бубнова-Галеркина дает практически точный результат при минимальных вычислительных затратах. Однако для сложных конструкций (с вырезами, подкреплениями, сложного очертания контура, с изломами поверхности) МКЭ является эффективным. Многие стандартные пакеты прикладных программ расчета строительных конструкций и сооружений основаны на МКЭ. Численные методы решения начальных задач для дифференциальных уравнений и их систем Математическими моделями задачи расчета элементов строительных конструкций являются не только краевые задачи для дифференциальных уравнений, но и начальные задачи. Кроме того, после применения метода Власова-Канторовича к смешанным задачам для дифференциальных уравнений в частных производных получается начальная задача для обыкновенного обыкновенных дифференциального дифференциальных уравнения уравнений. или Это для систем происходит при исследовании колебательных процессов для стержней, пластин и оболочек. Аналогичные задачи возникают при определении температурных полей для нестационарных задач теплопроводности. И в том, и в другом случае получается начальная задача для обыкновенных дифференциальных уравнений, где искомые функции зависят от временной координаты. Выше был рассмотрен вариант метода продолжения решения по параметру – метод последовательных нагружений. Расчетная схема этого метода является результатом применения метода Эйлера к начальной задаче для функционального уравнения. ЛЕКЦИЯ № 17 ПЛАН 1.1. Метод Эйлера 1.2. Метод Рунге-Кутта Метод Эйлера Пусть требуется найти решение обыкновенного дифференциального уравнения y′ = f ( x, y ), (111) удовлетворяющее начальному условию y ( x0 ) = y0 . (112) Предполагается, что функция f ( x, y ) в рассматриваемой области D , содержащей точку ( x0 , y0 ) , имеет непрерывные частные производные ∂f и ∂x ∂f до некоторого порядка m и решение задачи (111), (112) существует на ∂y всем заданном отрезке [ x0 , a ] . Отрезок интегрирования [ x0 , a ] разбивается на n одинаковых участков длиной h , и вводятся обозначения: xi = x0 + ih, (i = 0,1, ..., n). Метод Эйлера является простейшим численным методом решения начальной задачи (111), (112), которую называют задачей Коши. Суть метода заключается в том, что приближенные значения искомой функции y ( x) находят последовательно по формулам yi +1 = yi + Δyi , Δyi = hf ( xi , yi ) (i = 0,1, ..., n). (113) При этом искомая интегральная кривая y = y ( x) заменяется ломаной. Метод Эйлера прост, но имеет малую точность. Погрешность метода на каждом шаге имеет порядок Ο(h 2 ) , то есть метод первого порядка точности, поэтому при его использовании необходимо принимать малый шаг h. Более точным является усовершенствованный метод Эйлера, в котором на каждом шаге вначале вычисляется k1 = h f ( xi , yi ) , а затем Δyi = h f ( xi + h 2, yi + k1 2). (114) Этот метод дает на каждом шаге погрешность порядка Ο(h3 ) , то есть имеет второй порядок точности. ПРИМЕР 1 Найти приближенное решение на отрезке [0; 0,6] уравнения y′ = x 2 + y 2 при условии y (0) = 0 . Пусть h = 0,1 тогда в соответствии со схемой (113) имеем x0 = 0; y0 = 0; Δy0 = 0,1(02 + 02 ) = 0; y1 = y0 + Δy0 = 0; x1 = 0,1; Δy1 = 0,1(0,12 + 02 ) = 0,001; y2 = y1 + Δy1 = 0,001; x2 = 0, 2; Δy2 = 0,1(0,22 + 10−6 ) = 0,004; y3 = y2 + Δy2 = 0,005; x3 = 0,3; Δy3 = 0,1(0,32 + 25 ⋅ 10−6 ) = 0,009; y4 = y3 + Δy3 = 0,014; x4 = 0, 4; Δy4 = 0,1(0,42 + 1,96 ⋅ 10−4 ) = 0,0162; y5 = y4 + Δy4 = 0,0302; x5 = 0,3; Δy5 = 0,1(0,52 + 9,12 ⋅ 10−4 ) = 0,0251; y6 = y5 + Δy5 = 0,0553. Значение точного решения поставленной задачи при x = 0,6 дает результат y ( x) = 0,07244 . Рассмотрим применение метода Эйлера для решения начальных задач для систем обыкновенных дифференциальных уравнений. Пусть даны система обыкновенных дифференциальных уравнений dY = F ( x, Y ) dx (115) Y ( x0 ) = Y0 , (116) и начальные условия где Y = ( y1 , y2 , ..., ym )Τ , F = ( f1 , f 2 , ..., f m )Τ . Приближенное решение начальной задачи (115), (116) находят в виде Yi +1 = Yi + ΔYi , ΔYi = hF ( xi , Yi ) (i = 0,1, ..., n), (117) где Yi = Y ( xi ) . ПРИМЕР 2 Найти приближенное решение уравнения продольных колебаний прямоугольного стержня длиной l 2 ∂2w 2 ∂ w =a ∂x 2 ∂t 2 (118) при условиях w(t , 0) = w(t , l ) = 0, w(0, x) = w0 , ∂w = 0, ∂t t =0 (119) где w(t , x) – продольное смещение стержня. Так как заданы и краевые и начальные условия, то такую задачу называют смешанной. Сведем ее к начальной задаче. Искомую функцию w(t , x) представим в виде w(t , x) = T (t )sin kπ x . При l этом краевые условия в (119) будут выполнены. Из уравнения (118) получим обыкновенное дифференциальное уравнение T ′′ + bT = 0, (120) 2 ⎛ kπ ⎞ где b = ⎜ ⎟ , а штрихами обозначены производные по t . ⎝ al ⎠ Решение этого уравнения должно удовлетворять начальным условиям T (0) = w0 , T ′(0) = 0 . Уравнение (120) можно свести к системе уравнений T ′ = Q, Q′ = −bT . (121) Теперь для решения системы (121) на отрезке [0; 0,3] применим метод Эйлера, выбрав шаг h = 0,1 . В результате, используя (117), получим T0 = w0 ; Q0 = 0; ΔT0 = h ⋅ 0 = 0; ΔQ0 = h(−bw0 ) = −bw0 h; T1 = T0 + ΔT0 = w0 ; Q1 = Q0 + ΔQ0 = −bw0 h; ΔT1 = h(−bw0 h) = −bw0 h 2 ; ΔQ1 = h(−bw0 ) = −bw0 h; T2 = T1 + ΔT1 = w0 (1 − bh 2 ); Q2 = Q1 + ΔQ1 = −2bw0 h; ΔT2 = h(−2bw0 h) = −2bw0 h 2 ; ΔQ2 = h(−bw0 (1 − bh 2 )) = −bw0 h(1 − bh 2 ); T3 = T2 + ΔT2 = w0 (1 − 3bh 2 ); Q3 = Q2 + ΔQ2 = −3bw0 h + b 2 w0 h 2 . Таким образом, w(t , x) при t = 0,3 имеет вид w(0,3, x) = w0 (1 − 3bh 2 )sin kπ x . l Метод Рунге-Кутта Наиболее распространенным методом решения начальных задач типа (111), (112) является метод Рунге-Кутта. Основная схема этого метода имеет вид 1 yi +1 = yi + Δyi , Δyi = (k1( i ) + 2k2(i ) + 2k3(i ) + k4( i ) ), 6 (122) где (i ) 1 k = hf ( xi , yi ); k (i ) 2 ⎛ h k1( i ) ⎞ = hf ⎜ xi + , yi + ⎟; 2 2 ⎝ ⎠ ⎛ h k (i ) ⎞ k3(i ) = hf ⎜ xi + , yi + 2 ⎟ ; k4(i ) = hf ( xi + h, yi + k3( i ) ) 2 2 ⎠ ⎝ (i = 0,1, ..., n). (123) Погрешность приведенной схемы метода Рунге-Кутта на каждом шаге имеет Ο(h5 ) . ПРИМЕР Рассмотрим решение приведенной в примере 1 п. 1.1.1 начальной задачи методом Рунге-Кутта. В соответствии со схемой (122), (123) покажем начало процесса (i = 0) при шаге h = 0,1 : k1(0) = 0,1(0) = 0; k2(0) = 0,1(0,052 + 0) = 0,00025; k3(0) = 0,1(0,052 + 0,0001252 ) = 0,00025; k4(0) = 0,1(0,12 + 0,000252 ) = 0,001; 1 Δ y0 = (0 + 2 ⋅ 0,00025 + 2 ⋅ 0,00025 + 0,001) = 0,000333; 6 y1 = y0 + Δ y0 = 0,000333. Продолжая процесс до i = 5 , получаем y5 = 0,04179; x5 = 0,5; k1(5) = 0,1(0,52 + 0,041792 ) = 0,025175; k2(5) = 0,1(0,0552 + 0,0543782 ) = 0,030546; k3(5) = 0,1(0,0052 + 0,0570642 ) = 0,030576; k4(5) = 0,1(0,62 + 0,072367 2 ) = 0,036524; 1 Δ y5 = (0,025175 + 2 ⋅ 0,030546 + 2 ⋅ 0,030576 + 6 + 0,036524) = 0,030657; y6 = y5 + Δ y5 = 0,072448. Это значение практически не отличается от точного значения решения задачи при x = 0,6 , которое равно 0,07244. Если применить метод Рунге-Кутта для решения начальной задачи для системы уравнений (115), (116), то получим расчетную схему 1 Yi +1 = Yi + ΔYi , ΔYi = ( K1( i ) + 2 K 2( i ) + 2 K 3( i ) + K 4( i ) ), 6 где (i ) 1 K K (i ) 3 = hF ( xi , Yi ); K (i ) 2 ⎛ h K1(i ) ⎞ = hF ⎜ xi + , Yi + ⎟; 2 2 ⎠ ⎝ ⎛ h K 2( i ) ⎞ ( i ) (i ) = hF ⎜ xi + , Yi + ⎟ ; K 4 = hF ( xi + h, Yi + K 3 ) 2 2 ⎠ ⎝ (i = 0,1, ..., n). Например, для системы двух дифференциальных уравнений y′( x) = ϕ ( x, y, z ), z′( x) = ψ( x, y, z ) при начальных условиях y ( x0 ) = y0 , z ( x0 ) = z0 процесс нахождения решения y ( x) и z ( x) на i -м шаге сводится к следующим операциям. При заданных или найденных значениях h , xi , yi , zi вычисляют последовательно: k1(i ) = hϕ ( xi , yi , zi ); m1( i ) = hψ( xi yi , zi ); h k1( i ) m1(i ) k = hϕ ( xi + , yi + , zi + ); 2 2 2 h k1( i ) m1( i ) (i ) m2 = hψ( xi + , yi + , zi + ); 2 2 2 h k (i ) m(i ) k3(i ) = hϕ ( xi + , yi + 2 , zi + 2 ); 2 2 2 (i ) h k2 m2( i ) (i ) m3 = hψ( xi + , yi + , zi + ); 2 2 2 k4(i ) = hϕ ( xi + h, yi + k3( i ) , zi + m3( i ) ); (i ) 2 m4( i ) = hψ( xi + h, yi + k3(i ) , zi + m3( i ) ). Затем находят Δyi и Δzi по формулам 1 Δ yi = (k1(i ) + 2k2( i ) + 2k3(i ) + k4( i ) ), 6 1 Δ zi = (m1(i ) + 2m2( i ) + 2m3( i ) + m4(i ) ). 6 Значения y ( x) и z ( x) на шаге i + 1 определяются значениями yi +1 = yi + Δyi ; zi +1 = zi + Δzi . ЛЕКЦИЯ № 18 ПЛАН 1.1. Метод последовательных приближений 1.2. Методы решения жестких систем обыкновенных дифференциальных уравнений Метод последовательных приближений Решение начальной задачи (111), (112) можно получить в виде предела последовательных приближений, построенных на рекуррентной формуле x yn ( x) = y0 + ∫ f (t , yn−1 (t ))dt. (124) x0 Начальное приближение здесь может быть любым. Проще всего за начальное приближение принять y0 . Последовательность функций yn , вычисленных по формуле (124), равномерно сходится на заданном отрезке к искомому решению y ( x) , если функция f ( x, y ) в уравнении (111) удовлетворяет условиям теоремы существования и единственности решения. Практически процесс последовательных приближений следует прекратить, когда yn−1 и yn совпадают в пределах допустимой точности. Метод легко распространяется и на решение начальных задач для систем уравнений типа (115), (116). ПРИМЕР Рассмотрим решение задачи примера 2 п. 1.1.1 методом последовательных приближений. Как видим из примера, решение начальной задачи (118), (119) свелось к решению системы (121). В соответствии с формулой (124) последовательно получаем (помня, что T0 = w0 , Q0 = 0 ) t t T1 = T0 + ∫ Q0 dz = w0 ; Q1 = Q0 + ∫ (−bT0 )dz = −bw0t ; t ⎛ t ⎞ T2 = T0 + ∫ Q1dz = w0 ⎜1 − b ⎟ ; Q2 = Q0 + ∫ (−bT1 )dz = −bw0t ; 2⎠ ⎝ t 2 t ⎛ ⎛ t2 ⎞ t2 ⎞ T3 = T0 + ∫ Q2 dz = w0 ⎜1 − b ⎟ ; Q3 = Q0 + ∫ (−bT2 )dz = −bw0t ⎜ 1 − b ⎟ ; 2⎠ 6⎠ ⎝ ⎝ t ⎛ ⎛ t2 t4 ⎞ t2 ⎞ T4 = w0 ⎜1 − b + b 2 ⎟ ; Q4 = −bw0t ⎜1 − b ⎟ . 2 24 ⎠ 6⎠ ⎝ ⎝ Таким образом, решение поставленной задачи имеет вид 4 ⎛ t2 kπ x 2 t ⎞ . w(t , x) = w0 ⎜1 − b + b ⎟ sin 2 24 l ⎝ ⎠ Методы решения жестких систем обыкновенных дифференциальных уравнений При решении задач динамики для пластин и оболочек, когда учитывается, что система имеет несколько степеней свободы, после применения метода Власова-Канторовича получается система обыкновенных дифференциальных уравнений, которая относится к жестким системам. Кроме того, при решении нестационарных задач теплопроводности также получается система обыкновенных дифференциальных уравнений, которая относится к жестким системам. Решение жесткой системы содержит как быстро убывающие, так и медленно убывающие составляющие. Начиная с некоторого t > 0 ( t – параметр времени) решение системы почти полностью определяется медленно убывающей составляющей. Однако при использовании явных схем метода Рунге-Кутта быстро убывающая составляющая отрицательно влияет на устойчивость, что вынуждает брать шаг интегрирования τ = Δti слишком малым. Выход из этого положения был найден в применении неявных абсолютно устойчивых схем метода Рунге-Кутта. Так, если рассматривать систему уравнений dY = F (t , Y ) dt с начальными условиями, то неявная схема первого порядка точности имеет вид (Δti = τ) Yi +1 − Yi = F (ti +1 , Yi +1 ), τ а неявная схема второго порядка – Yi +1 − Yi = 0,5( F (ti , Yi ) + F (ti +1 , Yi +1 )), τ (125) или 3 1 Yi +1 − 2Yi + Yi −1 = τ F (ti +1 , Yi +1 ). 2 2 Можно записать и неявную схему четвертого порядка точности 25Yi +1 − 48Yi + 36Yi −1 − 16Yi −2 + 3Yi −3 = 12τF (ti +1 , Yi +1 ). При применении неявных схем Рунге-Кутта явно Yi +1 нельзя выразить и для нахождения Yi +1 приходится решать алгебраические уравнения, в общем случае нелинейные, одним из итерационных методов. В качестве примера рассмотрим решение уравнения теплопроводности для пластинки, которое после применения метода Власова-Канторовича сводится к системе уравнений dT0 a = (−1,5T0 + 0,5θ1 + 80); dt h 2 dθ1 12a = 2 (0,5T0 − 1,375θ1 − 35). dt h Здесь T0 (126) – температура срединной плоскости пластинки; θ – температурный коэффициент по толщине пластинки; θ1 = hθ ; h – толщина пластинки; a – коэффициент температуропроводности. Начальные условия для системы (126) имеют вид T0 (0) = 20°; θ1 (0) = 0°. Если использовать явную схему Рунге-Кутта первого порядка точности (метод Эйлера), то τ = Δti нужно брать равным 0,1 с, так как при τ = 1 с процесс неустойчив. При этом температура стабилизируется на значении T0 = 50,01° , θ1 = −7,05° по истечении 15 с, т. е. необходимо произвести 150 шагов (здесь принималось a = 0,222 см 2 с , h = 1 см , линейный размер пластинки 60h ). При использовании неявной схемы (125) шаг можно взять τ = 1 с и решение уже при t = 6 с полностью совпадает с решением, полученным с шагом τ = 0,1 с по явной схеме Эйлера. Процесс нахождения решения по схеме (125) принимает вид T0i +1 = 0,726T0i + 0,034θ1i + 13,71; θ1i +1 = 0,2345T0i + 0, 2345T0i +1 − 0, 2895θ1i − 32,824. и, как видим, не намного усложняется по сравнению с явной схемой метода Эйлера.