1 МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ Белорусский национальный технический университет Кафедра «Теоретическая механика» П. И. Ширвель ОСНОВЫ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ В МЕХАТРОНИКЕ Учебно-методическое пособие для студентов технических специальностей высших учебных заведений В 2 частях Часть 1 Рекомендовано учебно-методическим объединением по образованию в области машиностроительного оборудования и технологий Минск БНТУ 2015 2 УДК 372.853+372.862+519.3/6 ББК 22.2+22.3 Ш64 Рецензенты: Н. И. Юрчук, М. Г. Ботогова Ширвель, П. И. Основы метода конечных элементов в мехатронике : учебно- методическое пособие для студентов технических специальностей высших учебных заведений : в 2 ч. / П. И. Ширвель. – Минск : БНТУ, 2015– . – Ч. 1. – 2015. – 89 с. ISBN 978-985-550-506-9 (Ч. 1). Включено краткое изложение основных теоретических сведений, знание которых необходимо для сознательного решения задач механики и мехатроники численными и численно-аналитическими методами. При этом основное внимание уделено обоб- щению и систематизации данных по численным методам решения дифференциаль- ных уравнений, что в дальнейшем должно помочь правильно ориентироваться в на- правлении построения новых механико-математических моделей и эффективных расчетных методов решения конкретных прикладных задач. Учебно-методическое пособие составлено в соответствии с программой дисциплины «Метод конечных элементов в мехатронике» и содержит методические указания к выпол- нению практических и лабораторных работ, определения, формулы и индивидуальные задания, необходимые для самостоятельного выполнения лабораторных работ, а также основные термины и понятия, представленные с пояснениями и примерами. Издание предназначено для студентов IV курса машиностроительного факультета заочного и дневного отделений специальности 1-55 01 03 «Компьютерная мехатроника». УДК 372.853+372.862+519.3/6 ББК 22.2+22.3 ISBN 978-985-550-506-9 (Ч. 1) © Ширвель П. И., 2015 ISBN 978-985-550-507-6 © Белорусский национальный технический университет, 2015 Ш64 3 ПРЕДИСЛОВИЕ В настоящее время методы виртуального моделирования позволяют значительно ускорить процесс разработки элементов конструкций и ком- понентов оборудования различных мехатронных устройств, в том числе и снизить затраты, необходимые для изготовления опытных образцов и прототипов. Конструкции таких устройств определяются физическими принципами, на которых построена их работа, требованиями, предъявляе- мыми к надежности, чувствительности, быстродействию и т. д. Поэтому первостепенной целью моделирования является математический расчет и прогнозирование характеристик разрабатываемой конструкции. Наибо- лее эффективным и широко используемым средством достижения постав- ленной цели является применение метода конечных элементов (МКЭ), который позволяет с помощью компьютера в короткие сроки оценить па- раметры различных вариантов конструкции и выбрать наилучший. В связи с этим настоящее пособие имеет своей целью обучить компьютерному моделированию с помощью современных численных и численно-аналити- ческих подходов. В первую очередь курс посвящен изучению основ МКЭ, который является математической базой таких хорошо зарекомендовав- ших себя программных комплексов, как ANSYS, ABAQUS, LS-DYNA, Temper-3D, COMSOL Multiphysics и др. Также в издании представлены классические методы построения математических моделей. Этапы моде- лирования содержат подробное описание выполняемых операций, что поз- воляет студентам выполнять работу самостоятельно и, главное, осознанно. В учебно-методическом пособии основное внимание уделено матема- тическим аспектам и численно-алгоритмическим вопросам реализации МКЭ. Для решения задач применяются пакеты вычислительных программ (MathCAD, MatLAB, Maple, Mathematics), а также средства программирова- ния (Microsoft Visual Studio 2012). Такой подход позволяет избавиться от абстрактных вычислений, принятых в справочной литературе, и способ- ствует развитию у студентов здорового интереса именно к математиче- ским аспектам решения задач. Необходимость получить от машины ответ на поставленный вопрос требует более глубокого проникновения в суть изучаемой проблемы. Полученный ответ порождает новый вопрос, и в ре- зультате происходит закрепление теоретического материала. Выполнение предложенных лабораторных работ в дальнейшем даст возможность лучше разобраться в принципах конструирования и работы мехатронных устройств и освоить основные возможности современных методов их проектирования. В свою очередь, совокупность приобретенных навыков позволит студентам выполнять последующие проекты, в том чис- ле курсовые и дипломные работы, на более высоком научном уровне. 4 Лабораторная работа № 1 ПОСТРОЕНИЕ ФИЗИЧЕСКИХ, РАСЧЕТНЫХ И МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ФИЗИЧЕСКИХ ЯВЛЕНИЙ И ТЕХНИЧЕСКИХ ПРОЦЕССОВ. МЕТОДЫ РЕШЕНИЯ ОПРЕДЕЛЯЮЩИХ УРАВНЕНИЙ МАТЕМАТИЧЕСКОЙ МОДЕЛИ. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ. ОСНОВНЫЕ ЭТАПЫ ПРАКТИЧЕСКОЙ РЕАЛИЗАЦИИ Цель работы 1. Изучить основные термины и базовые понятия численного моде- лирования. 2. Провести анализ общих методов решения определяющих урав- нений математической модели. 3. Рассмотреть основные этапы практической реализации метода конечных элементов. Основные теоретические сведения Под математической моделью в механике и мехатронике в даль- нейшем будем подразумевать некую замкнутую систему математи- ческих соотношений, позволяющую с приемлемой точностью изу- чать интересующие особенности поведения рассматриваемого объ- екта. Такая замкнутая система уравнений получается с помощью методов дискретизации, при помощи которых непрерывная матема- тическая модель преобразуется в дискретную, состоящую из конеч- ного числа степеней свободы. Суть такой математизации – построе- ние моделей для описания различных явлений и изучения этих моде- лей с целью объяснения старых или предсказания новых эффектов. Как известно, исторически математика возникла и развивается как часть естествознания. Долгое время ее развитие существенным образом определялось потребностями прежде всего именно меха- ники и физики. По выражению выдающегося ученого эпохи Воз- рождения Леонардо да Винчи: «Наука механика потому столь бла- городна и полезна более всех прочих наук, что все живые тела, имеющие способность к движению, действуют по ее законам». Естественно, что прикладные исследования имеют непосредст- венную отдачу, а это усиливает доверие общества к механике, рас- 5 ширяет понимание ее проблем и, как следствие, способствует уве- личению вложения средств (с целью ее развития). Таким образом, уметь формулировать на языке математики конкретные задачи ме- ханики и мехатроники и означает в первую очередь способность построить математическую модель рассматриваемых процессов. Вычислительный (математический) эксперимент – метод изуче- ния конструкций или физических процессов с помощью виртуального (численного или математического) моделирования. Он предполагает, что за построением модели проводится ее численное исследование, которое позволяет проверить ее при различных условиях (рис. 1.1). Сфера применения – промежуточное положение между натурным экс- периментом и теоретическими исследованиями. Рис. 1.1. Схема организации вычислительного эксперимента Этапы вычислительного эксперимента:  построение математической модели (составление определяющих уравнений, описывающих модель);  выбор численных методов расчета (построение дискретной мо- дели, аппроксимирующей математическую задачу);  разработка вычислительного алгоритма; 6  создание или выбор программы, реализующей выбранный ал- горитм;  проведение расчетов и получение результатов вычисления;  анализ результатов вычисления, сравнение с физическими и на- турными испытаниями. Заметим, что в некоторых исследованиях доверие к результатам расчетов так велико, что при расхождении между данными числен- ных и натурных испытаний погрешность ищут именно в результа- тах натурных экспериментов. Следует подчеркнуть, что сейчас в технической литературе термин «модель» несколько перегружен. Как правило, в него вкладываются совершенно разные понятия. Как следствие, существуют значительные противоречия в определении таких понятий, как физическая, расчетная или математическая модель. Вместе с тем проблема создания моделей имеет фундаментальное значение для всех технических и естествен- ных дисциплин, поэтому далее попытаемся сформулировать по воз- можности развернутые определения употребляемых в этой области понятий, которые потом проиллюстрируем простейшими примерами. Под физической моделью какого-либо процесса, например про- цесса работы мехатронной конструкции при внешнем нагружении, следует понимать по возможности полное (в соответствии с достиг- нутым уровнем знаний) описание этого процесса в физически со- держательных терминах. В физическую модель должны входить без всяких упрощений все известные функциональные, дифференциаль- ные и прочие соотношения и связи между параметрами процесса. Эти связи могут носить как детерминированный, так и стохастиче- ский характер. Также физическая модель должна содержать имею- щиеся экспериментальные данные, относящиеся к рассматриваемо- му процессу, изложение гипотез, которые могут быть сформули- рованы по поводу еще неизученных связей и соотношений между параметрами системы. Другими словами, физическая модель пред- ставляет собой содержательное отражение реальных явлений или процессов на уровне современных знаний. Подчеркнем, что физи- ческая модель не может быть создана путем чисто эмпирического наблюдения данного класса явлений или процессов. Конечно, для создания любой физической модели (которая может рассматривать- ся как физическая теория) необходимы эксперименты и наблюде- 7 ния. Однако понимание самих экспериментов (и это следует особо подчеркнуть) невозможно без теории, т. е. без физической модели изучаемого класса явлений и процессов. Более того, планирование и организация новых экспериментов также невозможны без теории. Отметим, что в реальных физических системах или устройствах обычно встречаются не с одним каким-либо классом явлений или процессов, уже хорошо изученных современной наукой. Наоборот, обычно имеют дело с взаимодействием различных классов явлений и процессов, а также с новыми, еще недостаточно изученными явле- ниями наряду с неполнотой и неопределенностью исходной инфор- мации. В результате всего этого физические модели реальных объек- тов обычно оказываются весьма сложными и не вполне определен- ными, что сильно осложняет (или делает невозможным) их анализ. Все это и приводит к необходимости создания расчетных моделей. Расчетная модель также описывает процесс в физически содер- жательных терминах, но в отличие от физической модели в ней от- брошены параметры и не учитываются факторы, которые в заданных условиях и границах не оказывают заметного влияния на ход процес- са. При переходе от физической к расчетной модели сложные мате- матические зависимости или соотношения должны быть заменены по возможности более простыми приближенными или аппроксимиру- ющими соотношениями. В частности, во многих случаях переменные величины могут заменяться их средними постоянными значениями, нелинейные соотношения – линейными и т. д. По недостаточно изу- ченным связям или параметрам системы в расчетной модели могут вводиться аппроксимирующие гипотезы. При этом необходимо, что- бы все аппроксимации обеспечивали большую надежность системы в данных условиях ее работы. Следует заметить, что при всех упро- щениях и отбрасывании различных несущественных факторов или малых параметров необходимо соблюдать крайнюю осторожность, так как малые параметры могут оказывать громадное влияние на устойчивость или качественный характер процесса. Например, в тео- рии малых колебаний упругих систем часто говорят, что нелинейные члены ввиду их малости могут быть отброшены. Это утверждение в таком общем виде неверно. Под математической моделью процесса следует понимать уравне- ния и другие соотношения, приведенные в расчетной модели, алго- ритмы решения уравнений, составленные на их основе компьютерные 8 программы, схемы имитационного моделирования и т. д. При этом необходимо стремиться к эффективным математическим моделям. Это означает, что алгоритмы для решения уравнений должны быть по воз- можности простыми (но не в ущерб необходимой точности), носить универсальный характер, допускающий их удобное применение при различных граничных условиях, разнообразном характере внешних воздействий и т. д. Например, в механике и мехатронике такие алго- ритмы и схемы решения должны иметь адаптацию к различным усло- виям закрепления и нагружения объекта исследований. После того как введены математические объекты, называемые ма- тематическими моделями, возможно их изучение при помощи различ- ных формальных языков. Совершенные модели в современной техни- ке создаются путем последовательных приближений. Процесс усовер- шенствования моделей следует рассматривать как процесс их посте- пенного обогащения, доработки или проработки. Как правило, обычно начинают с простых моделей, которые могут сильно отличаться от действительности, а затем постепенно приходят к более совершенным моделям, точнее отражающим особенности изучаемого процесса. Методы решения определяющих уравнений математической модели В зависимости от физических и конструктивных особенностей объекта исследований, особенностей взаимодействия с окружающей средой, а также цели и требуемой точности исследований решение поставленных задач может быть получено различными методами. Аналитические методы. Получение решения построением в яв- ном виде функциональных аналитических зависимостей позволяет оценивать влияние на объект исследования различных факторов, од- нако это возможно лишь для узкого класса объектов исследований при простейших зависимостях в краевых условиях. Для обобщения результатов наиболее целесообразным является представление урав- нений в безразмерных относительных величинах. Переход к безраз- мерной форме обычно осуществляется методом теории размерности. На данный момент наиболее разработаны и описаны в литературе аналитические методы решения всех линейных задач [1]. Краткий обзор этих методов, выполненных на основе источников [2–4], пред- ставлен структурной схемой на рис. 1.2. 9 Рис. 1.2. Общая схема аналитических методов решения задач При помощи различных подстановок в аналитических методах используются приемы сведения исходной задачи к более простым, например: решение уравнения в частных производных – к решению обыкновенных дифференциальных уравнений, решение обыкновен- ных дифференциальных уравнений – к решению алгебраических и т. д. Использование каждого из рассматриваемых на схеме рис. 1.2 – Грина – для Ханкеля–Фурье – цилиндрических Шт ур ма –Л иу ви лл я 10 методов ограничивается видом области и краевых условий и исполь- зуемой системой координат. Так, например, метод разделения пере- менных (метод Фурье) можно использовать для классических задач при однородных граничных условиях; метод Дюамеля пригоден для неоднородных граничных условий, однако требует предварительного решения вспомогательной задачи; метод функций Грина снимает указанные ограничения, но построение этих функций требует опре- деленной изобретательности и в некоторых случаях трудновыполни- мо; метод конформных отображений используется для любой формы области, однако лишь для стационарных задач; в методах интеграль- ных преобразований требуется большая трудоемкость при обратном переходе от изображений к оригиналам. Заметим, что помимо точ- ных аналитических методов для предельных оценок находят приме- нение приближенные аналитические приемы некоторых из них, ос- нованные на специальных исследованиях условий улучшенной схо- димости рядов в точных аналитических решениях и последующем ограничении этих рядов. Все указанные выше аналитические методы в той или иной мере используются на практике для получения предельных оценок и пред- варительной проработки конструктивных решений, а также для обес- печения эксплуатационного контроля и управления. Отметим, что методы решения прямых нелинейных задач разра- ботаны слабо. Наиболее распространенным приемом является лине- аризация, т. е. сведение, где это возможно, нелинейной задачи к ли- нейной, например, с помощью интегральных преобразований или путем разложения нелинейных зависимостей в ряд Тейлора в окрест- ности среднеинтегрального значения с сохранением только линей- ной или какой-либо другой части разложения. Наибольшие трудно- сти возникают при решении нелинейных нестационарных задач в двух- и трехмерной постановке для тел с разрывными коэффици- ентами по пространству и переменными характеристиками во вре- мени. Подчеркнем, что аналитическая теория решения этих задач до сих пор не разработана. Численные методы. Внедрение компьютерных технологий спо- собствовало появлению значительного количества численных реше- ний, основанных на конечно-разностных, конечно-элементных, ва- риационных и других подходах. В настоящее время эти направле- ния активно развиваются, что стало возможным благодаря быстрому 11 совершенствованию вычислительной техники за последние десяти- летия. Это обусловлено потребностями многих областей фундамен- тальных и прикладных научных исследований, в том числе и меха- троники. Вообще говоря, интенсивное развитие компьютерных тех- нологий и разработка численных методов решения задач представ- ляют собой взаимно влияющие процессы, протекающие практически синхронно. В то же время одновременное существование большого количества методов говорит о том, что они не идеальны и работа в этом направлении должна быть продолжена. Схематический обзор наиболее развитых численных методов, на- ходящих применение в практике инженерных расчетов для деталь- ных исследований вариантов конструкций и режимов их эксплуата- ции, представлен на рис. 1.3. Кратко рассмотрим наиболее попу- лярные из них. Метод конечных разностей (МКР) основан на замене производ- ных в дифференциальных уравнениях и граничных условиях прибли- женными значениями, выраженными через разности значений функ- ций в конечном множестве отдельных дискретных точек (узлах) ис- следуемой области, например с помощью простейших соотношений линейной аппроксимации между соседними узлами и моментами времени. В общем случае для конечно-разностной аппроксимации используются и старшие члены ряда Тейлора при разложении иско- мых величин в узлах сеточной области или вводится необходимое число вещественных параметров (весов). В последнее время в отли- чие от традиционно применяемого фиксированного веса, дающего различные модификации расчетных схем (явная схема, неявная схе- ма, промежуточные схемы Кранка-Николса и др.), предложены схе- мы с «плавающим» весом, изменяющимся от точки к точке прост- ранственно-временной сетки. Специальные вопросы МКР, касающи- еся принципов построения эффективных сеток, условий сходимости и устойчивости решений, подробно рассматриваются в ряде моно- графий, в отдельных разделах многих книг и учебников, а также в большом количестве журнальных статей. Следует отметить только, что более простые и наглядные явные схемы требуют определенных соотношений между величинами шагов по времени и пространству для обеспечения сходимости и устойчивости решения, а более слож- ные неявные схемы являются абсолютно устойчивыми. 12 Рис. 1.3. Общая схема численных методов для решения краевых задач На На На Промежу- точные схемы Неявные схемы решения 13 Несмотря на относительную простоту, универсальность и возмож- ность получения априорных оценок точности, практическое использо- вание МКР для задач с разрывными коэффициентами, разрывными граничными условиями и сложными границами области затруднено. Методы аппроксимирующих функций позволяют искать прибли- женное решение краевой задачи в аналитической форме, например: 0 1 . n i i i x x A x     Аппроксимирующая система координатных функций ix строится из условия точного удовлетворения либо граничным условиям, либо дифференциальному уравнению, а входящие неизвестные парамет- ры iA определяются распространенными методами коллокаций (сов- падения), среднеквадратического приближения, методом Галеркина. Вариационные методы требуют замены дифференциальной крае- вой задачи некоторой вариационной, причем для построения функ- ционала разработаны специальные методы (энергетический, наимень- ших квадратов и др.). Это дает возможность, так же как и в методе аппроксимирующих функций, искать приближенные решения в ана- литической форме, применяя широко известные прямые методы, та- кие как метод Релея–Ритца, метод взвешенных невязок, метод Галер- кина. При применении двух последних методов в случае сложных граничных условий возникают значительные математические труд- ности, связанные с выбором системы координатных функций. Метод конечных элементов (МКЭ) объединяет конечно-разно- стный и вариационный методы в единый вариационно-разностный метод, основанный на аппроксимации непрерывной функции (напри- мер, перемещений, напряжений, скоростей, температур и т. д.) дис- кретной моделью (рис. 1.4), которая строится на множестве кусоч- но-непрерывных функций (конечных элементов). В качестве функ- ций элемента чаще всего применяется полином, порядок которого зависит от числа используемых в каждом узле элемента данных о непрерывной функции. Значения функций в узлах элементов на ранней стадии развития МКЭ определялись с помощью минимиза- ции функционала, связанного с рассматриваемым дифференциаль- ным уравнением. В дальнейшем более общие теоретические обос- 14 нования позволили исключить необходимость вариационной фор- мулировки физических задач и для вывода систем уравнений, опре- деляющих узловые значения, использовать методы взвешенных не- вязок, например, метод Галеркина. Одной из интересных неклас- сических модификаций вариационно-разностного метода является алгоритм, предложенный Л. А. Оганесяном. Рис. 1.4. Схема основной концепции МКЭ Наиболее важные преимущества МКЭ заключаются в возможно- сти решения задач с разрывными коэффициентами и разрывными граничными условиями, а также использования любых криволи- нейных элементов с переменными размерами для точного описания границ тела. Основным недостатком МКЭ является сложность по- лучения априорных оценок. История возникновения. 1950-е – возникновение МКЭ как чис- ленной процедуры решения задач строительной механики; 1960-е – связь МКЭ с процедурой минимизации позволила широко исполь- зовать его при решении задач в других областях техники; 1970-е – метод конечных элементов из численной процедуры решения задач 15 строительной механики превратился в общий метод численного ре- шения дифференциальных уравнений или систем дифференциаль- ных уравнений; 1980-е – разрабатываются графические пре-/пост- процессоры, решатели для нелинейных задач; 1990-е – инструмен- тальные средства МКЭ интегрируются в программное обеспечение систем автоматизированного проектирования. Основная идея метода конечных элементов состоит в том, что любую непрерывную величину можно аппроксимировать дискрет- ной моделью, состоящей из отдельных элементов (участков). На каж- дом из этих элементов исследуемая непрерывная величина аппрок- симируется кусочно-непрерывной функцией, которая строится на значениях исследуемой непрерывной величины в конечном числе точек рассматриваемого элемента. В общем случае непрерывная ве- личина заранее неизвестна, и нужно определить значения этой ве- личины в некоторых внутренних точках области. Дискретную мо- дель, однако, очень легко построить, если сначала предположить, что известны числовые значения этой величины в некоторых внут- ренних точках области (в дальнейшем эти точки мы назовем «узла- ми»). После этого можно перейти к общему случаю. Основные понятия. Чаще всего при построении дискретной мо- дели непрерывной величины поступают следующим образом. Об- ласть определения непрерывной величины разбивается на конечное число подобластей, называемых элементами. Эти элементы имеют общие узловые точки и в совокупности аппроксимируют форму об- ласти. В рассматриваемой области фиксируется конечное число то- чек. Эти точки называются узловыми точками или просто узлами. Значение непрерывной величины в каждой узловой точке перво- начально считается известным, однако необходимо помнить, что в действительности эти значения еще предстоит определить путем наложения на них дополнительных ограничений в зависимости oт физической сущности задачи. Используя значения исследуемой не- прерывной величины в узловых точках и ту или иную аппроксими- рующую функцию, определяют значение исследуемой величины внутри области. Аппроксимирующие функции чаще всего выбирают- ся в виде линейных, квадратических или кубических полиномов. Для каждого элемента можно подбирать свой полином, но полино- мы подбираются таким образом, чтобы сохранить непрерывность ве- 16 личины вдоль границ элемента. Этот полином, связанный с данным элементом, называют «функцией элемента». Как следует из основной концепции МКЭ, вся модель конструк- ции (или отдельной ее части) делится на множество конечных эле- ментов, соединенных между собой в вершинах (узлах), рис. 1.5, а, б. Силы также действуют в узлах. Конечный элемент не является «абсолютно жестким» телом. Имеются несколько наиболее употре- бительных типов конечных элементов (рис. 1.5, в): брус (А), стер- жень (В), тонкая пластина или оболочка (С), двумерное или трех- мерное тело (D). Естественно, что при построении модели могут быть использованы не один, а несколько типов элементов. Рис. 1.5. Представление конечных элементов Основные этапы практической реализации. Как было отмечено ранее, согласно МКЭ модель конструкции сложной формы подраз- деляется на более мелкие части (конечные элементы) сравнительно простой формы, в пределах которых ищется приближенное реше- ние. Результатом такого моделирования обычно является поле на- пряжений и смещений в целой конструкции. Таким образом, решение задачи с применением МКЭ состоит из следующих основных этапов:  создание геометрии модели, пригодной для МКЭ;  разбиение модели на сетку конечных элементов;  приложение к модели граничных условий;  численное решение системы уравнений (автоматически);  анализ результатов. б а в 17 Процесс конечно-элементного анализа включает в себя опреде- ленное количество последовательных шагов, базирующихся в ос- новном на двух подходах: математическом и физическом (рис. 1.6). Причем эти два подхода не являются взаимоисключающими, а до- полняют друг друга. Рис. 1.6. Схема конечно-элементного анализа Одним из способов настройки модели является дальнейшее уточ- нение модели. Основным источником ошибок является переход от фи- зической модели к математической: ошибки дискретизации и ошибки вычисления. Полная ошибка включает ошибки решения и погрешно- сти дискретизации. Ошибки (или погрешности), возникающие в мате- матической модели, устраняются на этапе верификации (проверка со- ответствия модели определяющим зависимостям и расчетным уравне- ниям, которые в нее заложены). При физическом подходе суммарная погрешность основана на результатах моделирования в сравнении с экспериментальными данными. Такая проверка соответствия модели процессам в реальном моделируемом объекте называется валидацией. Таким образом, для МКЭ кратко можно выделить следующие три ос- новных шага: идеализация, дискретизация и решение. Отметим, что каждый конечный элемент имеет собственную раз- мерность. Они могут описываться одной, двумя или тремя коорди- натами в зависимости от решаемой задачи (например, 1D, 2D или 3D). Причем время рассматривается как дополнительная координа- та. Каждый элемент описывается точками (узлы или узловые точ- ки). Геометрия КЭ определяется расположением узлов (рис. 1.7). 18 а б Рис. 1.7. Геометрия конечных элементов для двумерного (a) и трехмерного (б) пространства Общая схема алгоритма расчета по МКЭ задач механики и мехатроники 1. Дискретизация рассматриваемой области, т. е. замена непре- рывной среды совокупностью конечных элементов заданной формы, которые соединены в узлах непрерывным числом связей. Обычно при построении модели руководствуются предварительным опытом. В местах высоких концентраций величин конечно-элементную сет- ку сгущают. 2. Выбор вариационного принципа определения основных неиз- вестных функций, через которые устанавливаются остальные неиз- вестные. В задачах механики используются следующие основные принципы: принцип Лагранжа или принцип минимума потенциальной энер- гии (варьируются перемещения); принцип Кастельяно или принцип минимума дополнительной работы (варьируются напряжения); смешанный принцип Рейснера (варьируются напряжения и пе- ремещения); гибридный принцип Ху-Вашицу (варьируются напряжения, пе- ремещения, деформации). 3. Выбор аппроксимирующей функции. При этом предполагается, что перемещения внутри элемента могут быть выражены через пе- 19 ремещения его узлов. Связь описывается с помощью функции фор- мы. Такие функции должны удовлетворять определенным критери- ям, которые сводят результаты расчетов к точному решению с уве- личением числа конечных элементов (критерию полноты при стрем- лении функции к нулю; критерию совместности). 4. Реализация вариационного принципа. На этом этапе осуществ- ляется вычисление матриц жесткости элементов и построение гло- бальной матрицы жесткости для системы алгебраических уравнений и вектора узловых усилий. Глобальная матрица жесткости может быть получена несколькими методами, например, методом сложе- ния плоскостей или методом конгруэнтного преобразования. 5. Учет граничных условий. Полученная на предыдущем шаге матрица является вырожденной, так как часть уравнений системы окажутся взаимозависимыми. Корректировка этой матрицы приво- дит к невырожденной системе алгебраических уравнений. 6. Решение системы алгебраических уравнений. Особенность мат- рицы жесткости состоит в том, что она имеет ленточный вид. 7. Определение неизвестных характеристик (например, дефор- маций и перемещений). 20 Лабораторная работа № 2 ОПРЕДЕЛЕНИЕ ПЕРЕМЕЩЕНИЙ И ВЫЧИСЛЕНИЕ УГЛОВ РАЗВОРОТА ПРОИЗВОЛЬНЫХ ЭЛЕМЕНТОВ КОНСТРУКЦИЙ Цель работы: отработка численной реализации основных эта- пов МКЭ в пакете MathCAD. З а д а н и е 1. Методом конечных элементов выполнить расчет упругого из- гиба некоторого элемента конструкции с первоначальной прямой формой, один из концов которого жестко закреплен, а другой конец отогнут под действием некоторой силы (например: консольной бал- ки, стойки автомобиля, продольного сечения капота и т. д.). 2. Построить графические зависимости. 3. Проанализировать найденное численное решение, сравнив его с аналитическим. Варианты заданий приведены в таблице. Номера вариантов задания № варианта х1, мм х2, мм х3, мм I, мм 4 E, кг/мм2 u3, мм u3, мм 1 0 200 400 250 20 000 –20 25 2 0 300 600 100 25 000 –10 10 3 0 400 800 200 12 000 –15 15 4 0 250 500 150 10 000 –30 30 5 0 350 700 200 20 000 –10 10 6 0 150 300 150 15 000 –20 25 7 0 400 800 200 25 000 –30 30 8 0 50 100 100 30 000 –40 40 9 0 200 400 200 12 000 25 –25 10 0 300 600 150 10 000 10 –10 11 0 400 800 200 20 000 35 –35 12 0 250 500 150 15 000 30 –30 13 0 350 700 200 25 000 10 –10 21 Окончание таблицы № варианта х1, мм х2, мм х3, мм I, мм 4 E, кг/мм2 u3, мм u3, мм 14 0 150 300 150 30 000 20 –20 15 0 400 800 200 20 000 50 –50 16 0 50 100 100 15 000 40 –40 17 0 150 300 200 25 000 –10 10 18 0 400 800 150 30 000 –20 20 19 0 50 100 200 12 000 –50 50 20 0 200 400 300 10 000 –40 40 Основные понятия. Краткие теоретические сведения В лабораторной работе № 1 мы выяснили, что существует некий достаточно универсальный численный метод конечных элементов, который можно эффективно использовать для решения дифферен- циальных уравнений. В настоящей работе приведем пример его эф- фективного использования для расчетов классических задач меха- ники. С целью понять, как этот метод работает, знакомство с МКЭ проведем в достаточно простой форме, не углубляясь в математиче- ские аспекты, а лишь обозначая некоторые ключевые вопросы и наиболее важные проблемы его реализации [5–8]. В практике расчетов, связанных с обработкой эксперименталь- ных данных, вычислением ( )f x , разработкой вычислительных ме- тодов, часто встречаются следующие две основные ситуации: 1) как установить вид функции ( )y f x , если она неизвестна. При этом предполагается, что задана таблица ее значений  ( , ), 1,i ix y i m , которая получена либо из экспериментальных измерений, либо из сложных расчетов; 2) как упростить вычисление известной функции ( )f x или же ее характеристик ( ( ),max ( ))f x f x , если ( )f x слишком сложная. Ответы на эти вопросы даются теорией аппроксимации функций, основная задача которой состоит в нахождении функции ( )y x  , близкой (т. е. аппроксимирующей) в некотором нормированном про- странстве к исходной функции ( )y f x . Функцию ( )x выбирают 22 такой, чтобы она была удобной для последующих расчетов. Основ- ной подход к решению этой задачи заключается в том, что функцию ( )x выбирают зависящей от нескольких свободных параметров 1 2( , , ..., ),nс с с с т. е. 1( ) ( , , ..., ) ( , ),ny x x c c x c       значения которых подбираются из некоторого условия близости ( )f x и ( )x . Обоснование способов нахождения удачного вида функциональной зависимости ( , )x c  и подбора параметров c со- ставляет основную задачу теории аппроксимации функций. В рамках теории аппроксимации функции различают экстрапо- ляцию и интерполяцию. Под экстраполяцией понимается распространение результатов, полученных из наблюдений над одной частью некоторого явления, на другую его часть. Экстраполяция функции – продолжение функ- ции за пределы ее области определения, при котором продолженная функция (как правило, аналитическая) принадлежит заданному клас- су функций. Экстраполяция функций обычно происходит с помощью формул, в которых используется информация о поведении функции в некотором конечном наборе точек (в узлах экстраполяции), при- надлежащих ее области определения. Таким образом, экстраполяция – особый тип аппроксимации, при котором функция аппроксимируется вне заданного интервала, а не между заданными значениями. Задача интерполирования функций состоит в приближенной замене функции  f x более простой интерполирующей функцией  F x , значения которой в узлах интерполирования , 1, 2, ..., ,jx j n совпа- дают с соответствующими значениями  f x . На практике чаще всего интерполируют функции  f x , заданные таблично, в точках  1, 2, ...,jx j n , если необходимо узнать  f x при jx x . Обычно  F x отыскивают в виде обобщенного многочлена     1 , n n i i i P x c x    23 где  1, 2, ...,i i n  – линейно независимая система функций, а ic – действительные коэффициенты, определяемые из линейной алгеб- раической системы, например:     1 , 1, 2, ..., . n i k k k i c x f x y k n      Понятие интерполирования функций иногда употребляется в ка- честве противопоставления понятию экстраполирования, когда кон- структивно восстанавливаются (возможно, приближенно) значения функций в областях их определений. Отказавшись на время от ма- тематической терминологии, рассмотрим следующий пример: пусть есть некая абстрактная задача, например: найти некоторую зависи- мость y от x в интервале от a до b. Можно поступить двояко – ис- кать аналитический вид зависимости в виде функции y = y(x), т. е. в виде некоторой формулы, как, например, когда берется интеграл, или искать функцию в виде набора точек с некоторой нужной или заданной точностью (рис. 2.1). Рис. 2.1. Графическая иллюстрация аппроксимации произвольной функции Так, например, на рис. 2.1 показана некоторая искомая функция y = y(x), вместо которой найдена ломаная 1–2–3–4–5–6 без особого ущерба для точности. Понятно, что чем в большем числе точек ис- кать значение функции, тем точнее будет результат ее представле- ния ломаной. А сами методы и подходы, которые вместо аналитиче- ской зависимости находят искомую функцию (или много функций) в виде конечного числа точек (чисел), называются численными. 24 Для определенности на рис. 2.2 исследуемая балка уже разбита (или представлена) в виде двух (не будем усложнять) конечных эле- ментов, первый из которых имеет узлы 1 и 2, а второй – узлы 2 и 3. То есть узел 2 является общим для соседних конечных элементов. Рис. 2.2. Общая схема нагружения и разбиения произвольного объекта исследований Прежде чем приступить к решению конкретной задачи, остано- вимся на базовых свойствах конечного элемента. Для этого рас- смотрим некоторый абстрактный конечный элемент с произволь- ными узлами i и j (рис. 2.3). Далее конечным элементом будем называть такую часть элемента конструкции, чтобы перемещения любой его точки можно было бы выразить через перемещения его узлов без значительного ущерба для точности решения. Как легко заметить, при изгибе каждая точка конечного элемента с узлами ij имеет две степени свободы (точка может перемещаться по вертика- ли на расстояние u, а поперечное сечение в точке может развер- нуться на угол θ). Рис. 2.3. Схема двухузлового конечного элемента Так как конечный элемент представлен двумя узлами, то он име- ет четыре степени свободы. Значит, функцию перемещения в каж- 25 дой точке с координатой x внутри элемента можно представить по- линомом с четырьмя коэффициентами: 2 3 0 1 2 3 .u x x x        (2.1) Возникает вопрос: а почему полиномом, ведь истинная функция перемещения может быть сколь угодно сложной? Ответ: так проще. Кстати, полиномом первой степени (отрезком прямой) представлена функция y(x) внутри каждого конечного элемента, что было показа- но на рис. 2.1. А если нужна точность больше, то всегда можно уве- личить количество конечных элементов. Из теории упругости известно, что функция угла разворота сече- ния есть первая производная от функции перемещения, т. е. 2 1 2 3 d 2 3 . d u x x x         Положим, что начало системы координат расположено в узле i конечного элемента, а координатная ось x направлена вправо. Тогда выразим коэффициенты полинома 0 1 2 3, , ,    через перемеще- ния узлов i и j, т. е. , , , .i i j ju u  При x = 0, т. е. в узле i: 0 ,iu   1.i   Аналогично в узле j, т. е. при x = xj = l, где l = xj – xi – длина ко- нечного элемента, и с учетом полученных выше значений 0 и 1 2 3 2 3 2 3 2 3 , 2 . j i i j j j j i j j u u x x x x x               (2.2) Решая полученную систему (2.2), окончательно получаем значе- ния оставшихся коэффициентов полинома (2.1): 26     2 32 41 3 ;j i i j il u u l ll             23 41 2 .j i j i il l u u ll         Подставляя полученные выражения для 0 1 2 3, , ,    в формулу для u, получаем следующее равенство: 2 3 2 3 2 3 21 3 2 2i i x x x xu u x ll l l                     (2.3) 2 3 2 3 2 3 23 2 .j j x x x xu ll l l                   Формулой (2.3) выражается перемещение в любой точке конеч- ного элемента через перемещения в его узлах – основных неизвест- ных задачи. Подобным образом, разбивая конструкцию на конеч- ные элементы, вычисляем искомую функцию перемещений лишь в некоторых точках, а между ними представляем ее приближенно (аппроксимируя). Теперь, зная длину элемента l, перемещения в уз- лах конечного элемента , , ,i i j ju u  и координату какой-либо точ- ки x, по формуле (2.3) можно найти перемещение u этой точки. Громоздкую запись формулы (2.3) можно устранить, записав ее в матричном виде. Обозначим вектор узловых перемещений конеч- ного элемента  e в виде   i ie j j u u           . (2.4) 27 Введем вектор формы элемента  N в виде   2 3 2 3 2 3 2 2 3 2 3 2 3 2 1 3 2 2 . 3 2 x x l l x xx l lN x x l l x x l l                   Тогда выражение для перемещения u можно компактно записать в виде    ,T eu N  (2.5) где верхний индекс Т обозначает операцию транспонирования мат- рицы. Далее, считая деформации малыми, будем использовать их отно- сительную величину  , обозначая меру (степень) деформации об- щим выражением 0 0 0 ,L LL L L    (2.6) где L – конечная длина; L0 – исходная длина (до деформирования); L – отклонение вследствие деформации. Исследуемая конструкция имеет некоторую толщину (будем счи- тать ее постоянной по длине конечного элемента). Для наглядности на рис. 2.4 малый фрагмент элемента длиной dx до изгиба показан слева, а после изгиба – справа. 28 Рис. 2.4. Схема изгиба исследуемого элемента Как известно, при изгибе некоторое материальное волокно (ось) АВ не изменит своей длины, т. е. АВ = А1В1, а радиус его кривизны составит величину ρ. Тогда в соответствии с (2.6) относительная деформация произвольного волокна CD 1 1 .С D CD CD   (2.7) Из рис. 2.4 видно, что 1 1 d ,CD AB A B     а 1 1 ( )d .C D y    После подстановки в (2.7) получаем ( )d d . d y y         Как известно из математики, для малых величин углов разворота сечения справедливо выражение 2 22 222 3 d 1 dd . d d1 d u ux x u x           29 В таком случае получаем окончательное равенство для опреде- ления возникающих деформаций: 2 2 d . d uy x   (2.8) Формула (2.8) с учетом (2.5) примет вид        22d ,d T e ey N y Bx     где вектор  B   2 3 2 2 3 2 6 6 4 6 . 6 6 2 6 T x l l x l lB x l l x l l                   По теории упругости напряжения (силы, отнесенные к единице сечения) связаны с деформациями (по закону Гука) через модуль упругости (модуль Юнга) E. Тогда в нашем случае можно записать       .eE Ey B     (2.9) Из курса механики известно, что изгибающий момент связан с на- пряжениями в сечении, с учетом осевого момента инерции сечения, выражением 1 ,M I y   30 откуда получаем вектор узловых усилий в виде     .e eM EI B  (2.10) Таким образом, имеем: 1) вектор узловых перемещений  e конечного элемента (2.4); 2) вектор узловых сил  eM конечного элемента (2.10); 3) вектор относительных деформаций  (2.6); 4) вектор напряжений в сечении   (2.9). Заметим, что векторы   ,   и  eM выражаются через ком- поненты вектора  e , которые и являются неизвестными для дан- ной задачи. Общий принцип отыскания решения: если к некоторому узлу (или узлам) сетки конечных элементов (конечно-элементному ана- логу) приложить внешние силы или, что одно и то же, задать им некоторые перемещения, известные, например, из измерений дефор- мации рассматриваемой конструкции, то истинные перемещения остальных узлов будут такими, которые обеспечивают минимум полной энергии деформации. Минимизация полной энергии равносильна решению системы дифференциальных уравнений, описывающих деформацию балки. Осталось вычислить полную энергию. Произведение напряже- ний на соответствующие деформации дает удельную энергию де- формации (на единицу объема). Чтобы вычислить энергию дефор- мации по всему объему конечно-элементного аналога, надо проин- тегрировать (просуммировать) удельную энергию деформации по объему конечного элемента и просуммировать по всем элементам конечно-элементного аналога. Произведение узловых сил на соот- ветствующие узловые перемещения дают работу внешних сил. Пол- ная же энергия деформации есть разность энергии деформации всех конечных элементов конечно-элементного аналога с работой внеш- них сил. Чтобы найти ее минимум, надо записать выражение для полной энергии и продифференцировать его по каждому перемеще- 31 нию. Совокупность частных производных даст систему линейных уравнений в виде     ,K M  (2.11) где  K – матрица жесткости конечно-элементного аналога;   – вектор узловых перемещений конечно-элементного аналога;  M – вектор узловых сил конечно-элементного аналога. Формула (2.11) показывает, что энергия деформации конечно- элементного аналога равна работе внешних сил (работа внутренних сил рана работе внешних сил). В нашем случае выражение для матрицы жесткости конечного элемента eK   окончательно примет следующий вид: 3 2 3 2 2 2 3 2 3 2 2 2 12 6 12 6 6 4 6 2 . 12 6 12 6 6 2 6 4 e l l l l l ll lK EI l l l l l ll l                    Как видно, матрица жесткости конечного элемента eK   сим- метрична относительно главной диагонали. Используя обозначения узлов конечного элемента i и j, структуру этой матрицы можно представить в виде   , ii ije ji jj k k K EI k k                  32 где  , , ,ii ij ji jjk k k k           – подматрицы размером 2  2, а подста- новка вместо i и j номеров узлов конечного элемента показывает место каждой подматрицы в глобальной матрице жесткости  K конечно-элементного аналога. Аналогично структуру вектора  e перемещений и вектора уз- ловых усилий  eM конечного элемента можно представить в виде двух подвекторов, а подстановка вместо i и j номеров узлов конеч- ного элемента показывает место каждого подвектора в глобальном векторе узловых перемещений   и узловых сил  M для конеч- но-элементного аналога. П р и м е р Пусть конечный элемент № 1 (см. рис. 2.2) имеет узлы 1 и 2 с координатами 1 0x  и 2 100x  мм, а конечный элемент № 2 име- ет узлы 2 и 3 с координатами 2 100x  мм и 3 200x  мм, т. е. длина l = 100 мм одинакова для обоих конечных элементов. Для обоих конечных элементов примем модуль упругости E = 20000 кг/мм2, осевой момент инерции сечения I = 100 мм4. Сначала вычисляем компоненты матрицы жесткости первого и второго конечных элементов. Они будут одинаковы и равны между собой: 1 2 24 1200 24 1200 1200 80 000 1200 40 000 . 24 1200 24 1200 1200 40 000 1200 80 000 K K                  Далее распределяем подматрицы этих матриц по глобальной мат- рице жесткости в соответствии с номерами узлов, как показано на рис. 2.5. Суммируя получаем глобальную матрицу жесткости раз- мером 6  6 в виде 33   24 1200 24 1200 0 0 1200 80 000 1200 40 000 0 0 24 1200 48 0 24 1200 . 1200 40 000 0 160 000 1200 40 000 0 0 24 1200 24 1200 0 0 1200 40 000 1200 80 000 K                 Рис. 2.5. Схема распределения матриц жесткости 1-го и 2-го КЭ по глобальной матрице жесткости Глобальный вектор узловых усилий заполнен нулями (силы пока не заданы):   0 0 0 0 0 0 M              . 34 А неизвестный глобальный вектор узловых перемещений имеет вид   1 1 2 2 3 3 u u u               , тем самым получаем систему линейных уравнений вида     .K M  (2.12) Решение (2.12) пока тривиально, так как не заданы силы или пе- ремещения (или, как обычно говорят, не заданы граничные условия). Возвращаясь к рис. 2.2, видим, что в узле № 1 имеет место жесткое закрепление балки, т. е. перемещение и угол разворота сечения там равны нулю, или 1 0u  , 1 0  . Поэтому обнуляем первую и вторую строки и первый и второй столбцы матрицы жесткости ставим на главные диагонали по 1. Получаем следующий вид матрицы   :K   1 0 0 0 0 0 0 1 0 0 0 0 0 0 48 0 24 1200 . 0 0 0 160 000 1200 40 000 0 0 24 1200 24 1200 0 0 1200 40 000 1200 80 000 K               Легко заметить, что заданные граничные условия для такой сис- темы уравнений соблюдаются тождественно. Теперь из решения системы уравнений всегда будет 1 0u  и 1 0  . 35 Зададим перемещение узла № 3 вниз на 10 мм, т. е. 3 10u   . Для этого сформируем некоторый вектор  R в виде   0 0 0 0 10 0 R              и получим измененный вектор  M по алгоритму       M M K R  , т. е. вычтем из вектора  M произведение   K R , а результат по- местим снова в вектор  M . Теперь вектор  M выглядит следу- ющим образом:   0 0 240 12000 240 12000 M              . Далее обнулим пятую строку и пятый столбец матрицы  K , по- ставим в главной диагонали на пятой строке 1, в векторе  M на пятой строке поставим –10. Граничное условие внесено, а матрица и вектор  M окончательно выглядят следующим образом: 36   0 0 240 12000 10 12000 M              ;   1 0 0 0 0 0 0 1 0 0 0 0 0 0 48 0 0 1200 . 0 0 0 160 000 0 40 000 0 0 0 0 1 0 0 0 1200 40 000 0 80 000 K             Теперь можно решать систему уравнений. Решение нашей зада- чи имеет вид      1 0 0 3,125 0,056 10 0,075 K M                , или 1 0u  , 1 0  , 2 3,125u   мм, 2 0,056   рад, 3 10u   мм, 3 0,075   рад. Теперь, когда узловые перемещения найдены, можно подставить их в формулу для перемещения u и построить их график – форму деформированного объекта исследований (рис. 2.6). 37 Рис. 2.6. Изменение прогиба вдоль длины КЭ Сравнение полученного численного решения с точным аналити- ческим (упругая линия балки описывается известным уравнением третьей степени) дает полное совпадение результатов. Таким образом, в результате преобразований исходной системы уравнений фактически осталось только три неизвестных – переме- щение и угол разворота сечения в узле № 2 и угол разворота сече- ния в узле № 3. В процессе решения эти неизвестные найдены обес- печивающими минимум затрат энергии на деформацию. А любая иная их комбинация приведет к росту затраченной энергии. 38 Лабораторная работа № 3 МЕТОД ГАЛЕРКИНА Цель работы 1. Рассмотреть основные методы интегрирования линейных диф- ференциальных уравнений первого порядка. 2. Научиться с помощью метода Б. Г. Галеркина [9] решать обык- новенные дифференциальные уравнения численными методами; от- работать основные этапы указанного подхода. 3. Провести сравнение полученных приближенных численных решений с точными аналитическими результатами. З а д а н и е 1. Методом Галеркина решить задачу Коши для линейного обык- новенного дифференциального уравнения первого порядка ( ) ( )y p x y f x   на заданном отрезке  ,a b при условии, что ( ) .y a c 2. Найти приближенные функции численного решения исходной задачи (3) ( )y x и (4) ( )y x . 3. Получить точное аналитическое решение уравнения: найти функцию ( )y x . 4. Сравнить дискретные погрешности полученных решений (3) и (4) . 5. В пакете MathCAD построить графические зависимости ( )y x , (3) ( )y x и (4) ( )y x . 6. Рассмотреть сходимость вычислительного процесса. Варианты заданий приведены в табл. 3.1. 39 Таблица 3.1 Номера вариантов задания Номер варианта p(x) f(x) a b y(a) 1 –3 2xe 0 1 1 2 –2/x  2 /xe x x 1 2 1 3 –2/x x 1 2 1 4 –cos(x) sin(2 )x 0 1 1 5 –2/x 2 /xe x 1 2 1 6 –1/x ln 1x  1 2 1 7 tg(x) 2sin( )x 0 1 1 8 –1 cos( )x 0 1 0,5 9 (2x – 1)/x2 1 1 2 1 10 –2x 22 xxe 0 1 1 11 3/x x 1 2 1 12 –1/x /xe x 1 2 1 13 –1/x 3 / x 1 2 1 14 –x2 2x 0 1 2 15 1 x 0 1 1 16 1 xe 0 1 1 17 1,283 20,215( cos1,5 )x x 0,2 1,2 0,25 18 0,872 20,133( sin 2 )x x 0,2 1,2 0,25 19 2x 20,1x 0 1 0,8 20 2 x2 0 1 0,1 40 Основные понятия. Краткие теоретические сведения Пусть перед нами стоит задача решить уравнение ,Au f (3.1) где некоторый линейный оператор A определен на множестве D, плотном в некотором гильбертовом пространстве (линейном, век- торном в котором для любых двух элементов пространства опреде- лено скалярное произведение) H. Для его решения выбирается по- следовательность базисных функций { }j , обладающих следую- щими свойствами: 1) ,j D  1, ;j n 2) функции 1 2, , ..., n   линейно независимы при любом ;n 3) система базисных функций полна в H. Другими словами, данная полная система j является базисом, если каждый элемент пространства можно представить как линей- ную комбинацию элементов этой системы, и притом однозначно. Система гильбертова пространства является полной, если она по- рождает все пространство, т. е. если произвольный элемент прост- ранства может быть сколь угодно точно приближен по норме линей- ными комбинациями элементов этой системы. Теперь приближенное решение уравнения (3.1) строится в виде линейной комбинации базисных функций ( ) 1 ( ) ( ) nn j j j u x a x    (3.2) с постоянными коэффициентами ja , 1, .j n Коэффициенты ja определяются из условия, что выражение не- вязки (ошибка, погрешность в результате вычислений) ( )nR Au f  ортогонально (скалярное произведение двух элементов простран- ства равно нулю) в H ко всем функциям 1 2, , ..., .n   Таким обра- 41 зом получается следующая система линейных алгебраических урав- нений относительно :ja   1 , ( , ), n j k j kHj A a f      1, ,k n (3.3) где скалярное произведение определяется соотношением  , ( ) ( )d .j k j kH HA A х x x     Теперь после решения системы (3.3) приближенное решение за- дачи (3.1) находится по формуле (3.2). Этот метод впервые (в 1915 г.) был предложен русским ученым Борисом Григорьевичем Галерки- ным (1871–1945) и носит его имя. Кстати, родился этот ученый (наиболее известный своими результатами в области теории упру- гости) в Полоцке, а среднее образование получил в Минске. Указан- ный метод в литературе иногда называют методом Бубнова–Галер- кина, так как его применял в своих инженерных расчетах известный российский кораблестроитель, математик и механик И. Г. Бубнов. Позже (в 1942 г.) указанный метод был теоретически обоснован со- ветским математиком М.В. Келдышем. Одним из основных (и са- мых сложных) вопросов в методе Галеркина является вопрос о вы- боре базисных функций. В большом числе приложений метода Га- леркина в качестве базисных функций применяются полиномы низ- кого порядка. Ниже подробно разобран соответствующий пример. П р и м е р Методом Галеркина решить задачу Коши для дифференциально- го уравнения d 0 d y y x   (3.4) в области 0 1x  с граничным условием (0) 1y  . (3.5) 42 Р е ш е н и е Приближенное решение представим в виде суммы ( ) 1 ( ) 1 , nn j j j y x a x     (3.6) где первое слагаемое введено для удовлетворения граничному усло- вию (3.5). Базисные функции, таким образом, выберем из системы функций      2 21, , , ..., 0;1 ,nj x x x L   где  2 0;1L – пространство функций, суммируемых с квадратом на [0; 1]. Подставив функцию (3.6) в уравнение (3.4), получим невязку  1 1 1 . n j j j j R a jx x      Затем, приравнивая нулю скалярное произведение    1 11 1 1 1 10 0 , d d nk j j k k j j R x a jx x x x x x         , 1, ,k n получим систему линейных алгебраических уравнений относитель- но неизвестных коэффициентов ja ( 1, )j n в следующем матрич- ном виде: MA = B. (3.7) где элементы матриц М и В имеют следующий вид: 43  1 2 1, 0 1d , 1 j k j k k j jm jx x x j k j k          1 1 0 1d .kkb x x k   Легко заметить, что матрица М будет полнозаполненной. Далее, решая задачу при n = 3, получим следующий вектор-стол- бец коэффициентов :ja 1,0141 0,4225 0,2817 A        , подстановка которых в пробное решение (3.6) дает приближенное решение: (3) 2 3( ) 1 1,0141 0,4225 0,2817 .y x x x x    Нетрудно заметить, что задача (3.4), (3.5) имеет точное аналити- ческое решение ( ) .xy x e Сравнение полученного приближенного решения с точным в ше- сти точках 0,2ix i , 0,5,i  представлено в табл. 3.2. Таблица 3.2 Сравнение приближенного решения задачи (3.4), (3.5) с точным аналитическим x 0 0,2 0,4 0,6 0,8 1 (3) ( )y x 1 1,2220 1,4913 1,8214 2,2259 2,7183 (4) ( )y x 1 1,2214 1,4919 1,8221 2,2255 2,7183 ( )y x 1 1,2214 1,4918 1,8221 2,2251 2,7183 44 Используя данные табл. 3.2, вычислим дискретную погрешность:   1/26 2(3) (3) 1 ( ) ( ) 0,0013.i i i y x y x         Увеличив число слагаемых в (3.6) до n = 4 и разрешив матричное уравнение (3.7), получим, что  0,9990; 0,5095; 0,1399; 0,0699TA  , т. е. приближенное решение задачи (3.4), (3.5) теперь примет сле- дующий вид: (4) 2 3 4( ) 1 0,9990 0,5095 0,1399 0,0699 .y x x x x x     (3.8) Значения полученной функции (3.8) в отдельных точках отрезка [0; 1] также приведены в табл. 3.2. Вычисляя погрешность, получим   1/26 2(4) (4) 1 ( ) ( ) 0,00005.i i i y x y x         Таким образом, результаты проведенных вычислений показыва- ют, что с увеличением числа слагаемых n в приближенном пред- ставлении искомой функции погрешность быстро уменьшается. 45 Лабораторная работа № 4 РЕШЕНИЕ ЗАДАЧИ КОШИ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ Цель работы 1. Изучить основные этапы применения МКЭ к решению обык- новенных линейных дифференциальных уравнений (ОДУ) первого порядка. 2. Провести сравнение полученных приближенных (численных) решений с точными аналитическими результатами и данными, по- лученными по методу Б. Г. Галеркина. З а д а н и е 1. Методом конечных элементов решить задачу Коши для линей- ного ОДУ первого порядка ( ) ( )y p x y f x   на заданном отрезке  ,a b при условии ( ) .y a c 2. Найти приближенное решение ( ),y x разбивая отрезок  ,a b на 5 и 10 конечных элементов соответственно. Решив аналитически за- данное уравнение, найти функцию ( )y x . 3. Сравнить погрешности полученных решений (6) и (11) . 4. В пакете MathCAD построить графические зависимости ( )y x , ( )y x . 5. Проанализировать сходимость вычислительного процесса. Варианты заданий приведены в таблице. 46 Номера вариантов задания Номер варианта p(x) f(x) a b y(a) 1 –3 2xe 0 1 1 2 –2/x  2 /xe x x 1 2 1 3 –2/x x 1 2 1 4 –cos(x) sin(2 )x 0 1 1 5 –2/x 2 /xe x 1 2 1 6 –1/x ln 1x  1 2 1 7 tg(x) 2sin( )x 0 1 1 8 –1 cos( )x 0 1 0,5 9 (2x – 1)/x2 1 1 2 1 10 –2x 22 xxe 0 1 1 11 3/x x 1 2 1 12 –1/x /xe x 1 2 1 13 –1/x 3 / x 1 2 1 14 –x2 2x 0 1 2 15 1 x 0 1 1 16 1 xe 0 1 1 17 1,283 20,215( cos1,5 )x x 0,2 1,2 0,25 18 0,872 20,133( sin 2 )x x 0,2 1,2 0,25 19 2x 20,1x 0 1 0,8 20 2 x2 0 1 0,1 47 Краткие теоретические сведения Пусть требуется решить задачу Коши для линейного обыкновен- ного дифференциального уравнения первого порядка  , , 0F x y y  (4.1) на отрезке [a; b] с заданным граничным условием   .y a c (4.2) Применим к решению задачи метод Галеркина. В отличие от ла- бораторной работы № 1 выберем базисные функции в виде кусоч- но-линейных функций, систему которых обозначим  iN . Для этого разобьем отрезок [a; b] на участки (конечные элементы) точками 0 1 2, , , ... , ,nx a x x x b  которые называются узлами, и положим, что   1 0 11 00 1 , , 0, ; x x x x x x xN x x x       (4.3)   1 1 1 1 1 1 1 1 0, , , , , , 0, ; i i i i i i i i i i i i i x x x x x x x x x N x x x x x x x x x x                    1, 1,i n  (4.4)   11 1 1 0, , , ; n n n n n n n x x N x x x x x x x x          (4.5) 48 Функции ( )iN x , 1, ,i n линейно независимы и образуют полную систему функций в 2[ ; ].L a b Графики функций ( )iN x показаны на рис. 4.1. Рис. 4.1 (начало). Графическое представление функций ( )iN x 49 Рис. 4.1 (окончание) В таком случае приближенное решение задач (4.1), (4.2) пред- ставим в виде   0 ( ) ( ), n i i i y x y x N x     [ ; ],x a b где i – неизвестные коэффициенты. Найдем значения функции ( )y x в узлах конечных элементов ,jx x 0, .j n В силу того, что в точке jx x отлично от нуля только слагаемое, содержащее функцию ( )jN x , получим 0 ( ) ( ) ( ), n j i i j j j j i y x N x N x      0, .j n (4.6) 50 Кроме того, согласно (4.4) значение базисной функции ( )jN x в узле с номером j равно 1, т. е. ( ) 1.j jN x  Поэтому из (4.6) полу- чаем, что ( ) ,j j jy x y    т. е. неизвестные коэффициенты j совпадают со значениями функ- ции ( )y x в соответствующих узлах конечных элементов. Тогда функ- цию ( )y x можно записать в виде 0 ( ) ( ), n i i i y x y N x     [ ; ].x a b (4.7) Поскольку на каждом элементе 1[ ; ],k kx x 1, ,k n отличны от нуля только две базисные функции 1( )kN x и ( )kN x , то функцию ( )y x на отдельном элементе можно представить в виде 1 1( ) ( ) ( )k k k ky x y N x y N x     , 1[ ; ].k kx x x (4.8) Отсюда, в частности, следует, что 0 0 1 1( ) ( ) ( )y a y N a y N a    , но 0 ( ) 1,N a  1( ) 0N a  , значит, 0( )y a y c   . Таким образом, при- ближенное представление функции ( )y x в виде (4.7) обеспечивает точное выполнение граничного условия и введение дополнительно- го слагаемого не требуется. Для определения невязки подставим функцию (4.7) в уравнение (4.1), получим ( ) ( , , ),R x F x y y   здесь 0 d ( )d( ) , d d n i i i N xyy x y x x      [ ; ],x a b где согласно (4.3)–(4.5) производные базисных функций имеют сле- дующий вид: 51   0 10 1 0 1 1 , ,d d 0, ; x x xN x x x x x x       (4.9)   1 1 1 1 1 1 0, , 1 , , d 1d , , 0, ; i i i i ii i i i i i x x x x x x xN x x x x x x x x x                 1, 1,i n  (4.10)   1 1 1 0, , d 1 , .d n n n n n n x x N x x x xx x x         (4.11) Теперь, следуя методу Галеркина, приравняем к нулю скалярные произведения  ( ), ( ) ,kR x N x 0, ,k n тогда в силу линейности урав- нения (4.1) получим систему линейных алгебраических уравнений относительно ,iy 0, :i n  ( ), ( ) ( ) ( )d 0,bk k a R x N x R x N x x  0, .k n Учитывая, что 0y известно, первое уравнение системы можно исключить из рассмотрения, а в остальные подставить 0 .y c Уравнения системы можно упростить. Поскольку согласно (4.3)– (4.5) базисные функции ( ),kN x 1, 1,k n  отличны от нуля только на двух соседних элементах, содержащих узел kx , а функция ( )nN x – только на элементе 1[ ; ],n nx x интегрирование по отрезку [ ; ]a b мож- но заменить интегрированием по указанным элементам. Тогда сис- тема уравнений примет окончательный вид: 52 1 1 1 ( ) ( )d 0, 1, 1, ( ) ( )d 0. k k n n x k x x n x R x N x x k n R x N x x           (4.12) Следует отметить, что в силу (4.8), матрица этой системы линей- ных алгебраических уравнений имеет ленточную структуру. После решения системы (4.12) определяются значения аппрок- симирующей функции ( )y x в узлах ,ix x 1, .i n Затем, используя соотношение (4.7), можно найти приближенное значение искомой функции ( )y x в любой точке отрезка [ ; ]a b . Рассмотренный метод решения задачи (4.1), (4.2) называется ме- тодом конечных элементов. П р и м е р С помощью МКЭ решить задачу Коши для ОДУ первого порядка d 0 d y y x   в области 0 1x  с граничным условием (0) 1y  . Р е ш е н и е Разобьем отрезок [0; 1] точками 0 1 20, , , ... , 1nx x x x  на n ко- нечных элементов равной длины 1 / .x n  Приближенное решение задачи будем отыскивать в виде 0 ( ) ( ), n i i i y x y N x     где базисные функции ( )iN x задаются соотношениями (4.3)–(4.5). Согласно граничному условию 0 1.y  53 Невязка 0 d ( )d ( )( ) ( ) ( ) , d d n i i i i N xy xR x y x y N x x x           где производные d ( ) d iN x x определяются соотношениями (4.10). При- равняем нулю скалярные произведения:   1 00 d ( ) ( ), ( ) ( ) ( )d 0, d n i k i i k i N xR x N x y N x N x x x        1, .k n Заметим, что уравнение  0( ), ( ) 0R x N x  не рассматривается, так как известно, что 0 1.y  Учтем далее, что базисные функции ( ),kN x 1, 1,k n  отличны от нуля только на соседних элементах 1[ ; ]k kx x и 1[ ; ],k kx x  а функция ( )nN x – только на элементе 1[ ; ]n nx x , тогда эту систему линейных уравнений можно записать в виде 1 1 1 0 0 d ( ) ( ) ( )d 0, 1, ; d d ( ) ( ) ( )d 0. d k k n n x n i i i k ix x n i i i n ix N xy N x N x x k n x N xy N x N x x x                         Согласно (4.8) в суммах, которыми представляются функция ( )y x и ее производная на каждом элементе, сохраняются только два сла- гаемых. Этот факт, а также возможность поменять местами порядок интегрирования и суммирования в уравнениях позволяют записать систему в виде 1 1 1 1 1 1 d ( ) ( ) ( )d 0, 1, 1; d d ( ) ( ) ( )d 0 d k k n n xk i i i k i k x xn i i i n i n x N xy N x N x x k n x N xy N x N x x x                               54 или 1 1 1 1 1 1 1 1 1 1 1 1 1 d ( ) d ( ) ( ) ( )d ( ) ( )d d d d ( ) ( ) ( )d 0, 1, 1; d d ( ) d ( )( ) ( )d ( ) d d k k k k k k x x k k k k k k k k x x x k k k k x n n n n n n n N x N xy N x N x x y N x N x x x x N xy N x N x x k n x N x N xy N x N x x y N x x x                                                  1 1 ( )d 0. n n n n x x n x x N x x        Подставив выражения (4.3)–(4.5) и (4.9)–(4.11) в левые части урав- нений системы и осуществив интегрирование, окончательно получим 1 2 0 1 1 1 1 1 2 1 1 1 1 ; 3 2 6 2 6 1 2 1 0, 2, 1; 2 6 3 6 1 1 0. 2 6 2 3 k k k k k n n y y y x x y y y y y k n x x xy y                                                           Отметим, что матрица этой системы линейных уравнений имеет трехдиагональный вид, а значит, при численном решении ее можно использовать метод прогонки. Решая задачу при n = 5, т. е. разбивая отрезок от нуля до едини- цы [0;1] на пять конечных элементов и полагая, что 0,2x  , полу- чим следующее приближенное решение: 0 1 2( ) ( ) 1,2488 ( ) 1,4997 ( )y x N x N x N x    3 4 51,8557 ( ) 2,2441 ( ) 2,7620 ( ),N x N x N x   где числовые коэффициенты представляют собой значения функции ( )y x в узлах 0; 0,2; 0,4; 0,6; 0,8;1.ix  55 Вычисляя дискретную 2[0;1]L погрешность (см. лабораторную работу № 3), получим   1/26 2(6) (4) 1 ( ) ( ) 0,065.i i i y x y x         На рис. 4.2 сплошной линией показан график функции ,xy e пунктирной – график функции ( )y y x  , крестиками – вычисленные значения 0 ( ).y x Рис. 4.2. Сравнение результатов расчетов решения исходной задачи Повторив все вычисления, разбивая отрезок [0;1] на десять ко- нечных элементов точками 0;1 ,ix i  0;10,i  получим следующее приближенное решение: 0 1 2 3 4 5 6 7 8 9 10 ( ) ( ) 1,0987 ( ) 1,2205 ( ) 1,3428 ( ) 1,4899 ( ) 1,6409 ( ) 1,8190 ( ) 2,0049 ( ) 2,2210 ( ) 2,4495 ( ) 2,7120 ( ). y x N x N x N x N x N x N x N x N x N x N x N x               Погрешность вычисления составила   1/211 2(11) (4) 1 ( ) ( ) 0,020.i i i y x y x         56 Сравнение результатов двух вариантов расчетов показывает, что с ростом числа конечных элементов погрешность уменьшается. Од- нако следует подчеркнуть, что применение метода Галеркина с ба- зисными функциями в виде многочленов (см. лабораторную работу) дает более точные результаты. Лабораторная работа № 5 ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ Цель работы: научиться решать численно, методом конечных эле- ментов, краевые задачи для дифференциальных уравнений второго порядка. З а д а н и е 1. Методом конечных элементов решить краевую задачу для урав- нения Штурма–Лиувилля на заданном отрезке [ , ]a b при заданных граничных условиях. 2. Найти приближенное решение ( ),y x разбивая отрезок [ , ]a b на 10 и 20 конечных элементов. 3. Решить задачу с помощью компьютерного математического пакета MathCAD (Matlab, Maple, Mathematics) и сравнить получен- ные результаты, вычисляя дискретные (11) и (21) погрешности (см. лабораторную работу № 3). 4. Вычислить значения производной d d y x  в точках x a и x b . 5. Построить графические зависимости ( )y x , ( )y x . 6. Проанализировать сходимость вычислительного процесса. Варианты заданий приведены в таблице. 57 Номера вариантов задания № вар. p(x) q(x) f(x) a b y(a) y(b) 1 arctgx lnlnx lnx 3 7 0,49 –0,12 2 1 xe /3xe 2sin ( / 4)x 0 4 0,7 0 3 xe cos( / 3)x /2xe –2 1 –0,27 0,4 4 1+ x 1+ 2xe 3 1 x 0 1,5 2,6 –0,12 5 1 x 1 2 x x   x 1,5 3,5 0,29 –0,17 6 cos( /10)x 2 1 1x x  1 (1 )x 2,4 4 0,66 –0,36 7 21 sin x arctg( )x sin( )x 0,5 4 –1,5 0,64 8 1 / x 1 x 21 / x 1 7 1,2 –0,29 9 x 1 1 x 2 1 1 x 1,5 2,5 0 0,98 10 2 1 1x x  cos( /10)x cos( /10)x –2 0 2 0,1 11 2 1x x  12 x  x 2 3 –0,5 0,71 12 ln( 5)x  1 x 2 / (1 )x 0 2,5 0,78 –0,12 13 cos( / 2)x 2 – cos( )x sin( / 2)x 0 2 0,18 –0,37 14 21 / ( 1)xe  21 x 2x –1 3 0,49 –0,17 15 1 2 x x   cos( /10)x cos( )x π 4 0 0,74 16 2 /2xe 1 / x xe 1 2 –0,89 0,67 17 2 1 1xe  arctg x 21 x x  3 5 –0,61 0,49 58 Окончание таблицы № вар. p(x) q(x) f(x) a b y(a) y(b) 18 1 1 x sin xe x 2 4 0,95 0,88 19 2 1 1x  xe  arctg x –1 1 0,76 –0,20 20 xe 2 1 1 2cos x sin( )x 0 1 0,16 0,44 Основные понятия. Краткие теоретические сведения Пусть на отрезке [ , ]a b требуется найти решение уравнения Штур- ма–Лиувилля, представленного в виде  ( ) ( ) ( )p x y q x y f x   (5.1) и удовлетворяющего граничным условиям ( ) ,y a c ( ) .y b d (5.2) Предполагается, что ( )p x , ( )q x , ( )f x – непрерывные функции, известные изначально, c и d – постоянные. После дифференцирова- ния первого слагаемого в уравнении (5.1) его можно записать в виде 2 2 d d( ) ( ) ( ) ( ). dd y yp x p x q x y f x xx     Согласно методу конечных элементов разобьем отрезок [ , ]a b точками (узлами) ,i b ax a in   0; ,i n на заданное число n ко- нечных элементов и приближенное решение задачи (5.1), (5.2) бу- дем искать в виде 59 0 ( ) ( ) ( ), n i i i y x y x y N x     [ , ],x a b (5.3) где базисные функции ( )iN x задаются соотношениями (4.3)–(4.5), ,iy 0; ,i n – неизвестные значения аппроксимирующей функции ( )y x в узлах ix . Краевым условиям при x a и y b можно удовлетворить не- посредственно заданием соответствующих значений в крайних уз- лах. Однако удобнее составить уравнения метода Галеркина без пер- воначального задания граничных значений функции ( )y x , а затем учесть эти условия при решении окончательной системы уравне- ний. Таким образом, на данном этапе все значения 0 1, , ... , ny y y   бу- дем рассматривать как неизвестные. Отметим также, что поскольку функции ( ),iN x 1, 1,i n  отличны от нуля только на двух соседних элементах  1,i ix x и  1,i ix x  , то на каждом элементе функцию ( )y x можно записать в виде суммы только двух слагаемых, учитывающих вклад двух базисных функ- ций, отличных от нуля на этом элементе: 1 1( ) ( ) ( ),i i i iy x y N x y N x      1,i ix x x , 1, 1.i n  (5.4) Теперь определим невязку: 2 2 d d( ) ( ) ( ) ( ) ( ) dd y yR x p x p x q x y f x xx        и запишем ее скалярные произведения на базисные функции ( ),kN x 1, 1:i n    2 2d ( ), ( ) ( )d ( ) ( )dd d ( ) ( ) ( )d ( ) ( ) ( )d ( ) ( )d . d b b k k k a a b b b k k k a a a y xR N R x N x x p x N x x x y x p x N x x y x q x N x x f x N x x x                 (5.5) 60 Следует отметить, что в силу линейности функции ( )kN x произ- водные второго порядка от этих функций равны нулю. Проинтегри- ровав по частям в первом интеграле правой части выражения (5.5), можно понизить порядок производной, входящей в подынтеграль- ное выражение: 2 2 d ( ) d ( ) d ( )( ) ( )d ( ) ( ) ( ) ( ) d dd dd d( )d ( ) ( )d , 1, 1. d d d b k k k a b b k k a a y x y b y ap x N x x p b N b p a N a x xx Ny yp x x p x N x x k n x x x                (5.6) В соотношении (5.6) ( ) ( ) 0k kN a N b  при 1, 1k n  , 0 ( )N b  ( ) 0,nN a  0 ( ) ( ) 1,nN a N b  т. е. слагаемые, не содержащие ин- тегралы, сохраняются при 0k  и .k n Подставляя (5.6) в (5.5) и приравнивая скалярные произведения нулю, получим систему уравнений 0 0 0 dd ( ) d( ) ( )d ( )d ( )d ; d d d dd ( )d ( )d ( )d , 1, 1; d d dd ( ) d( ) ( )d ( )d ( )d . d d d b b b a a a b b b k k k a a a b b b n n n a a a Ny a yp a p x x yN q x x N f x x x x x Ny p x x yN q x x N f x x k n x x Ny b yp b p x x yN q x x N f x x x x x                             Учитывая, что функции ( ),kN x 1, 1,k n  отличны от нуля только на двух соседних элементах  1,k kx x и  1,k kx x  , а функции 0 ( )N x и ( )nN x – только на первом и последнем элементах соответственно, интегрирование по всему отрезку [ , ]a b можно заменить интегриро- ванием по указанным элементам. Тогда система примет вид 61 1 1 1 0 0 0 1 1 1 1 1 1 1 1 0 0 0 dd ( ) d( ) ( )d ( )d ( )d ; d d d dd ( )d ( )d ( )d , 1, 1; d d dd ( ) d( ) ( )d ( )d d d d k k k k k k n n n n n x x x x x x x x x k k k x x x x x n n n x x x Ny a yp a p x x yN q x x N f x x x x x Ny p x x yN q x x N f x x k n x x Ny b yp b p x x yN q x x N x x x                                    1 ( )d . nx f x x    Подставим в уравнения системы функцию ( )y x , учитывая при этом соотношение (5.4):   1 0 1 1 0 0 1 1 0 01 0 1 0 0 1 1 0 0 1 1 1 1 1 1 d ddd ( ) ( ) ( )d d d d d ( )d ( )d ; d d d ( )d d d d d d d ( )d d d d k k k k x x x x x x x k k k k k x x k k k k k x k k k N NNy a p a y y p x x x x x x y N y N N q x x N f x x N N Ny y p x x x x x N N Ny y p x x x x x y N y                                                  1 1 1 1 1 1 1 1 1 1 1 1 1 ( )d ( )d ( )d , 1, 1; d d dd ( ) ( ) ( )d d d d d ( )d ( )d . k k k k k k n n n n n n x k k x x x k k k k k k x x x n n n n n x x x n n n n n n x x N N q x x y N y N N q x x N f x x k n N N Ny b p b y y p x x x x x x y N y N N q x x N f x x                                              62 Представляя далее интегралы от сумм как суммы интегралов и вы- нося постоянные ,ky 0, ,k n за знаки интегралов, получим систему линейных алгебраических уравнений относительно :ky 1 1 0 0 1 1 1 0 0 0 1 1 2 20 0 0 01 1 1 0 0 1 1 1 dd ( ) ( ) ( )d ( )d d d dd ( )d ( )d ( )d ; d d d d ( )d ( )d d d d( d k k k k x x x x x x x x x x x x k k k k k x x k k Ny a p a y p x x N q x x x x NNy p x x N N q x x N f x x x x N Ny p x x N N q x x x x Ny x                                           1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 ( )d ( )d d d ( )d ( )d d d ( )d , 1, 1; d dd ( ) ( ) ( )d ( )d d d d k k k k k k k k k k n n n n x x k x x x x k k k k k x x x k x x x n n n n n x x p x x N q x x N Ny p x x N N q x x x x N f x x k n N Ny b p b y p x x N N q x x x x x                                            1 1 1 2 2d ( )d ( )d ( )d . d n n n n n n x x x n n n n x x x Ny p x x N q x x N f x x x                 Теперь из этой системы уравнений можно вычеркнуть первую и последнюю строки, получающиеся при применении метода Га- леркина к исходному уравнению с базисными функциями, отвеча- ющими узлам, где заданы значения решения, и использовать задан- ные граничные значения 0 0y y c  и n ny y d  в оставшихся уравнениях. Таким образом получим систему линейных алгебраи- ческих уравнений с трехдиагональной матрицей: (5.7) 63 2 2 2 1 1 1 2 2 1 1 1 1 0 0 2 21 2 1 1 1 2 01 2 1 1 0 1 1 1 1 d d d( )d ( )d ( )d d d d dd( )d ( )d ( )d ( )d ; d d d d ( )d ( )d d d k x x x x x x x x x x x x x x k k k k k x N N Ny p x x N q x x y p x x x x x NNN N q x x N f x x с p x x N N q x x x x N Ny p x x N N q x x x x                                        1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 1 2 2 1 d ( ( )d d d d ( )d ( )d ( )d d d ( )d , 2, 2; d d ( )d ( )d d d k k k k k k k k k k k k k x x x k k x x x x x k k k k k k x x x x k x n n n n n Ny p x x x N NN q x x y p x x N N q x x x x N f x x k n N Ny p x x N N q x x x                                                  1 1 2 2 2 2 2 1 1 2 21 1 1 1 1 1 d ( )d ( )d ( )d d d dd ( )d ( )d . d d n n n n n n n n n n n n n n x x x x x x x n n n n x x x x x n n n n x x x Ny p x x N q x x N f x x x N N p x x N N q x x x x                                              Численно решить эту систему уравнений можно методом про- гонки [10, 11]. После определения всех узловых значений аппроксимирующей функции ( )y x с помощью соотношения (5.3) можно найти прибли- женные значения искомой функции ( )y x в любой точке отрезка  ,a b . Заметим, что вычеркнутые из системы (5.7) уравнения можно использовать для определения значений производной d d y x  в точках x a и x b . 64 П р и м е р Реализация решения в пакете MathCAD 65 A0 0 16.525 X0 X2 xdN x 1( )( )2 p x( ) d 16.428 X0 X2 xN x 1( )( )2 q x( ) d 0.097 Aj 1 j 1 Xj 1 Xj 1 xdN x j( )( )2 p x( ) d Xj 1 Xj 1 xN x j( )( )2 q x( ) d Aj 1 j Xj Xj 1 xdN x j 1( ) dN x j( ) p x( ) d Xj Xj 1 xN x j 1( ) N x j( ) q x( ) d An 2 n 3 Xn 2 Xn 1 xdN x n 2( ) dN x n 1( ) p x( ) d Xn 2 Xn 1 xN x n 2( ) N x n 1( ) q x( ) d An 2 n 2 Xn 2 Xn xdN x n 1( )( )2 p x( ) d Xn 2 Xn xN x n 1( )( )2 q x( ) d B0 X0 X2 xf x( ) N x 1( ) d c X0 X1 xdN x 0( ) dN x 1( ) p x( ) d X0 X1 xN x 0( ) N x 1( ) q x( ) d    j 2 n 2 Bj 1 X j 1 X j 1 xf x( ) N x j( ) d Bn 2 Xn 2 Xn xf x( ) N x n 1( ) d d Xn 2 Xn xdN x n 1( ) dN x n( ) p x( ) d Xn 1 Xn xN x n 1( ) N x n( ) q x( ) d    A 16.525 8.342 0 0 0 0 0 0 0 8.342 17.125 8.635 0 0 0 0 0 0 0 8.635 17.705 8.919 0 0 0 0 0 0 0 8.919 18.265 9.193 0 0 0 0 0 0 0 9.193 18.809 9.46 0 0 0 0 0 0 0 9.46 19.338 9.72 0 0 0 0 0 0 0 9.72 19.852 9.973 0 0 0 0 0 0 0 9.973 20.354 10.22 0 0 0 0 0 0 0 10.22 20.843      B 2.592 0.276 0.29 0.303 0.316 0.329 0.341 0.352 1.415      66 Y lsolve A B( ) YT 0.365 0.412 0.432 0.427 0.396 0.339 0.255 0.144 2.615 10 3  yp x( ) c N x 0( ) 0 n 2 i Yi N x i 1( )    d N x n( ) Given 2x y x( )d d 2 p x( )( ) x p x( )d d x y x( )d d  q x( ) y x( ) f x( ) y a( ) c steps 105 y b( ) d y Odesolve x b steps( )  1 n 1 i y Xi  yp Xi  2    9.127 10 4 Сравнение результатов дано на рис. 5.3. 1.5 2 2.5 3 3.5 0.2 0 0.2 0.4 0.6 y x( ) yp x( ) x Рис. 5.3. Сравнение результатов 67 Лабораторная работа № 6 РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ В ЧАСТНЫХ ПРОИЗВОДНЫХ С ПОМОЩЬЮ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ НА ПРИМЕРЕ РЕШЕНИЯ КРАЕВОЙ ЗАДАЧИ ДЛЯ УРАВНЕНИЯ ЛАПЛАСА Цель работы: решить уравнение Лапласа 2 2 2 2 0 ( ) u uu f x x y       , (6.1) в области D с границей С при заданных граничных условиях 1( , ) ( , )u x y f x y , 1( , )x y C , (6.2) 2( , ) ( , ) u x y f x y n   , 2( , )x y C , (6.3) где 1 2C C C  , 1( , )f x y , 2 ( , )f x y – заданные функции; n – внешняя нормаль к части границы 2C . Решение реализовать численно методом конечных элементов. З а д а н и е 1. Методом конечных элементов решить смешанную краевую задачу для уравнения Лапласа в области ABCD, показанной ниже на рис. 6.1. Здесь AD и BC – дуги окружностей радиусов R и 1 соответ- ственно, φ – заданный угол. На рис. 6.1 показана также рекомендуемая нумерация конечных элементов. Граничные условия задаются в вариантах индивидуаль- ных заданий (приведены в таблице). 68 Рис. 6.1. Рекомендуемая нумерация КЭ при решении задачи Варианты индивидуальных граничных условий к заданию № варианта AB BC CD DA R φ 1 0nu  1nu   lnu R 1 /nu R 3 30 2 0nu  1nu   0nu  lnu R 3 30 3 0nu  0u  0nu  1 /nu R 3 30 4 lnu r 1nu   0nu  1 /nu R 3 30 5 lnu r 1nu   0nu  lnu R 3 30 6 0nu  0u  lnu r 1 /nu R 2,5 25 7 lnu r 0u  0nu  1 /nu R 2,5 25 8 0nu  1nu   lnu r lnu R 2,5 25 9 lnu r 1nu   lnu r 1 /nu R 2,5 25 10 0nu  0u  0nu  lnu R 2,5 25 11 0nu  1nu   lnu R 1 /nu R 2 10 69 Окончание таблицы № варианта AB BC CD DA R φ 12 0nu  1nu   0nu  lnu R 2 10 13 0nu  0u  0nu  1 /nu R 2 10 14 lnu r 1nu   0nu  1 /nu R 2 10 15 lnu r 1nu   0nu  lnu R 2 10 16 0nu  0u  lnu r 1 /nu R 1,5 15 17 lnu r 0u  0nu  1 /nu R 1,5 15 18 0nu  1nu   lnu r lnu R 1,5 15 19 lnu r 1nu   lnu r 1 /nu R 1,5 15 20 0nu  0u  0nu  lnu R 1,5 15 2. Вычислить погрешность 2 1 ( ( , ) ) m точ i i i i u x y u     , где m – число узлов, в которых найдено значение функции ( , );u x y точ ln( );u r r – радиус-вектор точки ( , ).x y 3. Найти значение функции 2 2x yV V V  в узле с координатами  0,(1 ) / 2R . Сравнить с точным значением точ 2 / (1 ).V R  Счи- тать, что x uV x   , y uV y   . Основные понятия. Краткие теоретические сведения Пусть требуется решить уравнение Лапласа (6.1) в области D с границей С при заданных граничных условиях в виде условий Ди- рихле (6.2) и Неймана (6.3). 70 Для решения задачи методом конечных элементов область D раз- бивают на отдельные участки – конечные элементы (КЭ). Чаще всего в качестве КЭ выбирают треугольные или четырехугольные области. В соответствии с этим и конечные элементы называются треугольны- ми или четырехугольными элементами. На каждом элементе неизвест- ную функцию аппроксимируют многочленом некоторой степени от переменных x и y. В зависимости от степени многочлена конечные элементы называются линейными, квадратичными и т. д. Рассмотрим применение МКЭ с линейными треугольными эле- ментами к решению задачи (6.1)–(6.3). Разобьем область D (рис. 6.2) на треугольные подобласти и про- нумеруем все получившиеся треугольные элементы: 1 2, , ..., MD D D и их вершины 1 2, , ..., , , , ..., .lP P P P P P   Эту нумерацию вершин бу- дем называть глобальной, а сами вершины треугольников, как и ра- нее, – узлами конечных элементов. D C Pα Pβ Pγ Dj Рис. 6.2. Сравнение результатов расчетов решения исходной задачи Рассмотрим один произвольный элемент jD с вершинами , ,P P P   и границей jC . Кроме глобальной нумерации введем локальную нумерацию узлов. Пусть узел P с глобальным номером  имеет локальный номер 1, узел P – номер 2 и узел P – номер 3. Локаль- ную нумерацию узлов следует проводить на каждом элементе про- тив часовой стрелки (рис. 6.3). 71 Рис. 6.3. Локальная нумерация узлов треугольного конечного элемента Искомую функцию на элементе jD аппроксимируем многочле- ном первой степени ( ) ( ) ( ) ( ) 1 2 3( , ) ( , ) j j j ju x y u x y x y       , ( , ) .jx y D (6.4) В дальнейшем верхние индексы (j) будем опускать. Обозначим узловые значения неизвестной функции 1 1 1( , )u x y u , 2 2 2( , )u x y u , 3 3 3( , )u x y u , тогда из выражения (6.4) можно получить 1 1 1 2 1 3 2 1 2 2 2 3 3 1 3 2 3 3 , , . u x y u x y u x y                   (6.5) Здесь и далее вместо знаков приближенного равенства исполь- зуются знаки точного равенства. Разрешив систему (6.5) относительно коэффициентов , 1,3i i  , получим 1 2 3 2 3 1 3 1 2 1 1 3 2 2 1 3 3 2 1 2 1 2 3 3 2 2 3 1 1 3 3 1 2 2 1 3 ( ) ( ) ( ) , 2 ( ) ( ) ( ) , 2 ( ) ( ) ( ) , 2 j j j u y y u y y u y y S u x x u x x u x x S u x y x y u x y x y u x y x y S                   (x1, y1) (x3, y3) (x2, y2) 72 где 1 1 2 2 3 3 1 1 1 2 1 j x y S x y x y  – площадь треугольника jD . Подставив найденные коэффициенты i в (6.5), получим       1 1 1 1 2 2 2 2 3 3 3 3 1 1( , ) 2 u x y u a x b y c u a x b y c u a x b y c S          , где использованы следующие обозначения: ,i k na y y  ,i k nb x x  ,i k n n kc x y x y  (6.6) а индексы , ,i k n изменяются по правилу круговой перестановки чисел 1, 2, 3. Обозначим  1( , ) 2i i i ij N x y a x b y c S    , 1, 2, 3.i  В таком случае функцию ( , )u x y можно записать в виде 3 1 ( , ) ( , ),i i i u x y u N x y    ( , ) .jx y D (6.7) Отметим, что функции ( , )iN x y обладают следующими свойст- вами: ( , ) 1,i i iN x y  1,3,i  ( , ) 0,i k kN x y  , 1,3, .i k i k  В системе координат , , ix y N функции ( , )iN x y представляют собой плоские треугольники с вершинами в точке ( , ,1)i ix y и в двух смежных узло- вых точках ( , )k kx y и ( , ),n nx y , , 1,3, .i k n i k n   На стороне тре- угольника между узлами ( , )k kx y и ( , )n nx y функции ( , )iN x y обра- щаются в нуль. Такие функции ( , )iN x y называются локальными ба- зисными функциями. График функции 2 ( , )N x y показан на рис. 6.4. 73 Рис. 6.4. Пример графика локальной базисной функции при k = 2 После подстановки функции (6.6) в уравнение (6.1) получим не- вязку 3 1 ( , ) ( , ),i i i R u x y u N x y      ( , ) .jx y D Для минимизации невязки на элементе jD применим метод Га- леркина: найдем скалярные произведения невязки на базисные функ- ции и приравняем их нулю: 3 1 0, j i i k i D u N N d      1,3.k  (6.8) Следует отметить, что поскольку ( , ),iN x y 1,3,i  – линейные функции, то их частные производные второго порядка обращаются в нуль. Для понижения порядка производных, входящих в (6.8), применим известную формулу Грина [5]: d d d , j j j i i k k i k D C D NN N N s N N n          , 1,3,i k  где iN n   – производная по направлению внешней нормали к гра- нице jC . Последний интеграл в правой части этого равенства можно вы- числить: (x1, y1) (x2, y2) (x3, y3) N2(x1, y1) 74 ( )d , 4 j ji k i k i k ki jD a a b bN N g S      где величины , , ,i k i ka a b b определены соотношениями (6.6); jS – площадь элемента jD . Тогда система уравнений (6.8) примет вид 3 3( ) 1 1 d , j j i i i kki i i C Nu g u N s n      1,3.k  (6.9) Правые части уравнений системы можно преобразовать: 3 3 ( ) 1 1 d d d d . j j j ji i i k i k k k i iC C C N N uu N s u s N s N s h n n n               Далее интеграл по границе jC элемента jD представим в виде суммы интегралов по всем трем сторонам треугольника: (2) (3) (1) (1) (2) (3) d d d d . j k k k k C u u u uN s N s N s N s n n n n             Один из этих интегралов, не содержащий узел с номером k , об- ратится в нуль, а для приближенного вычисления двух других мож- но вынести из-под знаков интеграла среднее значение нормальной производной nu . Например, при 2k  будем иметь (2) (2) (1,2)( ) 2 cp 2 (1) (1) (3) (3) (2,3)( ) 2 cp 2 (2) (2) (1) 2 (3) d d ; d d ; d 0, j n j n u N s u N s n u N s u N s n u N s n            75 где (1,2)( ) cp jnu и (2,3)( ) cp jnu – средние значения производной un   на сто- ронах элемента jD между указанными узлами. Не составляет труда вычислить интегралы ( )(2) (1,2) 2 (1) d , 2 jl N s  ( )(3) (2,3) 2 (2) d . 2 jl N s  Здесь ( )(1,2)jl и ( )(2,3)jl – длины сторон треугольника jD между уз- лами с локальными номерами 1, 2 и 2, 3 соответственно. Таким образом, получим, что  ( ) ( ) ( )(1,2)( ) (2,3)( ) cp cp2 (1,2) (2,3)12j j jj jn nh u l u l  . Теперь систему (6.9) можно записать в виде 3 ( ) ( ) 1 ,j ji ki k i u g h   1,3,k  или в матричном виде ( ) ,jGU H где  ( )jkiG g – матрица системы;  ( ) ( ) ( )( ) 1 2 3, , Tj j jjU u u u – матрица-столбец неизвестных;  ( ) ( ) ( )1 2 3, , Tj j jH h h h – матрица столбец правых частей. Следует отметить, что матрица G – симметричная. Полученная система уравнений в качестве неизвестных содер- жит только узловые значения искомой функции ( , )u x y на элемен- те jD . Такие системы линейных алгебраических уравнений запи- сываются для каждого элемента и называются локальными систе- мами уравнений. 76 Поскольку одни и те же узлы конечных элементов могут принад- лежать нескольким смежным элементам, то оказывается, что одни и те же неизвестные войдут в несколько локальных систем. Для того чтобы избежать многократного определения одних и тех же неиз- вестных, все локальные системы уравнений объединяют в одну – глобальную систему линейных алгебраических уравнений, например, ,FU Z где  , , 1, ,kiF f k i L  – глобальная матрица системы;  1 2, , ..., TjU u u u – матрица-столбец правых частей уравнений; L – количество узлов конечных элементов во всей области D. Для построения матрицы этой системы суммируются все эле- менты матриц локальных систем, соответствующие одним и тем же узловым значениям неизвестной функции. При этом можно исполь- зовать следующий прием, который проиллюстрируем на примере. Пусть область D имеет вид, изображенный на рис. 6.5. Разобьем ее на треугольные конечные элементы. Нумерация элементов (но- мера заключены в квадратики) и узлов (глобальные номера заклю- чены в кружочки, а локальные номера узлов показаны внутри каж- дого элемента) также показана на рис. 6.5. Рис. 6.5. Нумерация элементов и узлов при разбиении области на треугольные элементы 77 Из рис. 6.5 видно, что элементы матрицы глобальной системы уравнений  , , 1,9,kiF f k i  можно выразить через элементы ло- кальных матриц следующим образом: (1) 11 33 ;f g (1) 21 13 ,f g (1) (2) (3)22 11 11 33 ;f g g g   31 0,f  (3)32 13 ,f g (3) (4)33 11 33 ;f g g  (1) 41 23 ,f g (1) (2)42 12 13 ,f g g  43 0,f  (1) (2) (5)44 22 33 22 ;f g g g   51 0,f  (2) (3)52 12 23 ,f g g  (3) (4)53 12 23 ,f g g  (2) (5)54 23 23 ;f g g  (2) (3) (4) (5) (6) (7) 55 22 22 22 33 33 22 ;f g g g g g g      61 62 64 0,f f f   (4)63 13 ,f g (4) (7)65 12 23 ,f g g  (4) (7) (8)66 11 33 33 ;f g g g   71 72 73 76 0,f f f f    (5)74 12 ,f g (5) (6)75 13 23 ,f g g  (5) (6)77 11 22 ;f g g  81 82 83 84 0,f f f f    (6) (2)85 13 12 ,f g g  (7) (8)86 13 23 ,f g g  (6)87 12 ,f g (6) (7) (8) 88 11 11 22 ;f g g g   91 92 93 94 95 97 0,f f f f f f      (8)96 33 ,f g (8)98 12 ,f g (8)99 11 .f g Отметим, что индексы элемента глобальной матрицы совпадают с глобальными номерами узлов, взаимосвязь которых учитывает этот элемент. В силу симметричности матриц локальных систем уравнений сим- метричной будет и матрица глобальной системы линейных уравне- ний (поэтому выше записаны только элементы, стоящие на главной диагонали и ниже ее). Кроме того, глобальная матрица имеет разре- женный вид (т. е. имеет много нулевых элементов) и ленточную структуру, причем ширина ленты зависит от нумерации узлов. Заме- тим, что если максимальная разница в номерах узлов любого элемен- та равна m – 1, то матрица F будет иметь ширину ленты (2m – 1). В рассмотренном выше примере ширина ленты равна 7. При решении задачи в более сложных областях можно составить таблицы соответствия глобальных и локальных узлов для всех эле- ментов и с их помощью сформировать матрицу глобальной системы уравнений [4]. 78 Аналогично формируются и правые части глобальной системы уравнений: (1) 1 3 ;z h (1) (2) (3)2 1 1 3 ;z h h h   (3) (4)3 1 3 ;z h h  (1) (2) (5)4 2 3 2 ;z h h h   (2) (3) (4) (5) (6) (7) 5 2 2 2 3 3 2 ;z h h h h h h      (4) (7) (8)6 1 3 3 ;z h h h   (6.10) (5) (6) 7 1 2 ;z h h  (6) (7) (8)8 1 1 2 ;z h h h   (8)9 1 .z h При вычислении элементов  , 1,9,kz k  следует помнить, что при суммировании интегралов, взятых по одной и той же внутрен- ней стороне треугольника, но в противоположных направлениях, сумма их обращается в нуль. Тогда правая часть глобальной систе- мы будет содержать только интегралы по внешним границам обла- сти. Например, 5 0,z   (3) (1)(1) (3) (1) (3)(1,3)( ) (3,1)(3)2 cp cp1 3 (3,1) (3,1) (1) (3) 1d d . 2 j n n u uz N s N s u l u l n n        Следующий этап решения задачи состоит в учете граничных условий (6.2), (6.3). Условия Неймана (6.3) можно учесть непосред- ственно в правых частях (6.10) при формировании правых частей системы уравнений. Условия же Дирихле (6.2) задают значения са- мой неизвестной функции, поэтому можно подставить в систему эти заданные значения и исключить из рассмотрения соответству- ющие уравнения. После решения задачи из исключенных уравнений можно найти нормальную производную в узлах на соответствую- щих участках границы. Разрешив полученную после введения граничных условий систему уравнений, можно найти значения неизвестной функции во внутрен- них узлах и тех граничных узловых точках, в которых задано гра- ничное условие Неймана. Краевые задачи вида (6.1)–(6.3) часто возникают в различных обла- стях науки и техники. Например, в гидродинамике идеальной жидко- сти функция ( , )u x y имеет смысл потенциала скоростей, который свя- 79 зан со скоростью течения жидкости соотношением grad .V u В таком случае для определения компонент скорости перемещения точек эле- мента в этом случае можно воспользоваться соотношениями 3( ) 1 1 ; 2 j x i i ij uV a u x S     3( ) 1 1 , 2 j y i i ij uV b u y S     где коэффициенты , ,i ia b 1,3,i  определены соотношениями (6.6). Отметим, что значения компонент скорости оказываются посто- янными. Скорость перемещения узловой точки с глобальным номе- ром  определяется как среднее арифметическое значений скоро- сти элементов, имеющих этот общий узел: ( ) ( ) 1 1 ;jx x j V V      ( ) ( ) 1 1 j y y j V V      , где  – число элементов, имеющих общий узел .P Пример реализации численного решения задачи в MathCAD Рассматривается случай следующих граничных условий: : 0,nAB u  : 0,BC u  : ln ,CD u r : 1/ ,nDA u R 1,5,R  о15 . 80 81 x 6 3 0 y6 3 y2 2 x7 1 x1 1 y7 1 y1 1 x 7 2 0 y 7 2 y 2 2 x 7 3 0 y7 3 1 x8 1 x3 1 y8 1 y3 1 x8 2 x1 1 y8 2 y1 1 x 8 3 0 y8 3 1 x 0.324 0.324 0.259 0 0.388 0.324 0.324 0.259 0 0 0 0 0 0.388 0 0.324 0.388 0 0.324 0.259 0 0 0 0    y 1.207 1.207 0.966 1 1.449 1.207 1.207 0.966 1.5 1.25 1.25 1.25 1.5 1.449 1.25 1.207 1.449 1.5 1.207 0.966 1.25 1.25 1 1    i 1 8 ai 1 y i 2 y i 3 b i 1 x i 2 x i 3 c i 1 xi 2 yi 3 xi 3 yi 2 a i 2 y i 3 y i 1 bi 2 x i 3 x i 1 c i 2 xi 3 yi 1 xi 1 yi 3 c i 3 xi 1 yi 2 xi 2 yi 1 a i 3 y i 1 y i 2 bi 3 x i 1 x i 2 Si xi 1 xi 2 xi 3 yi 1 yi 2 yi 3 1 1 1   2  k 1 3 j 1 3 g i j k( ) ai k ai j bi k bi j 4 Si  f1 1 g 1 3 3( ) f2 1 g 1 1 3( ) f2 2 g 1 1 1( ) g 2 1 1( ) g 3 3 3( ) f3 1 0 f3 2 g 3 1 3( ) f3 3 g 3 1 1( ) g 4 3 3( ) f4 1 g 1 2 3( ) f4 2 g 1 1 2( ) g 2 1 3( ) f 4 3 0 82 f4 4 g 1 2 2( ) g 2 3 3( ) g 5 2 2( ) f 5 1 0 f5 2 g 2 1 2( ) g 3 2 3( ) f5 3 g 3 1 2( ) g 4 2 3( ) f5 4 g 2 2 3( ) g 5 2 3( ) f5 5 g 2 2 2( ) g 3 2 2( ) g 4 2 2( ) g 5 3 3( ) g 6 3 3( ) g 7 2 2( ) f 6 1 0 f 6 2 0 f6 3 g 4 1 3( ) f 6 4 0 f6 5 g 4 1 2( ) g 7 2 3( ) f6 6 g 4 1 1( ) g 7 3 3( ) g 8 3 3( ) f 7 1 0 f 7 2 0 f 7 3 0 f7 4 g 5 1 2( ) f7 5 g 5 1 3( ) g 6 2 3( ) f 7 6 0 f7 7 g 5 1 1( ) g 6 2 2( ) f 8 1 0 f 8 2 0 f 8 3 0 f 8 4 0 f8 5 g 6 1 3( ) g 7 1 2( ) f8 6 g 7 1 3( ) g 8 2 3( ) f8 7 g 6 1 2( ) f8 8 g 6 1 1( ) g 7 1 1( ) g 8 2 2( ) f 9 1 0 f 9 2 0 f 9 3 0 f 9 4 0 f 9 5 0 f9 6 g 8 1 3( ) f 9 7 0 f9 8 g 8 1 2( ) f9 9 g 8 1 1( ) i 1 9 j 1 9 f i j f j i f 0.98 0.724 0 0.256 0 0 0 0 0 0.724 2.089 0.592 0 0.773 0 0 0 0 0 0.592 1.141 0 0 0.549 0 0 0 0.256 0 0 1.96 1.448 0 0.256 0 0 0 0.773 0 1.448 4.179 1.185 0 0.773 0 0 0 0.549 0 1.185 2.282 0 0 0.549 0 0 0 0.256 0 0 0.98 0.724 0 0 0 0 0 0.773 0 0.724 2.089 0.592 0 0 0 0 0 0.549 0 0.592 1.141    l i j k( ) xi j xi k 2 yi j yi k 2 r i j( ) x i j 2 y i j 2 u7 ln r 1 3( )( ) u8 ln r 1 1( )( ) u9 ln r 3 1( )( ) u 3 0 u 6 0 83 uT 0 0 0 0 0 0 0.405 0.223 0( ) н е и з в е с т н ы е : u1 u2 u4 u5 un1 3 1 0 un3 3 1 0 un1 3 2 1 R  un5 2 1 1 z1 1 2 un13 2 l 1 3 2( ) un13 1 l 1 3 1( )  z2 1 2 un13 1 l 1 3 1( ) un33 1 l 3 3 1( )  z2 0 z4 1 2 un13 2 l 1 3 2( ) un52 1 l 5 2 1( )  z3 0 z5 0 z4 0.065 A f1 1 f2 1 f4 1 f5 1 f1 2 f2 2 f4 2 f5 2 f1 4 f2 4 f4 4 f5 4 f1 5 f2 5 f4 5 f5 5    B1 z1 f1 3 u3 6 9 j f1 j uj    B2 z2 f2 3 u3 6 9 j f2 j uj    B3 z4 f4 3 u3 6 9 j f4 j uj    B4 z5 f5 3 u3 6 9 j f5 j uj    q lsolve A B( ) u1 q1 u2 q2 u4 q3 u5 q4 uT 0.27 0.136 0 0.14 0.115 0 0.405 0.223 0( ) 84 ut2 ln R 1( ) 2   ut1 ln R( ) ut3 ln 1( ) ut4 ln R( ) ut5 ln R 1( ) 2   ut6 ln 1( ) ut7 ln R( ) ut8 ln R 1( ) 2   ut9 ln 1( ) utT 0.405 0.223 0 0.405 0.223 0 0.405 0.223 0( )  1 9 i uti ui 2    0.329 t1 1 u2 t 1 2 u 4 t1 3 u1 t2 1 u2 t2 2 u 5 t 2 3 u 4 t3 1 u 3 t 3 2 u 5 t3 3 u2 t4 1 u6 t4 2 u 5 t 4 3 u 3 t5 1 u7 t 5 2 u 4 t5 3 u5 t6 1 u8 t6 2 u 7 t 6 3 u 5 t7 1 u 8 t 7 2 u 5 t7 3 u6 t8 1 u9 t8 2 u 8 t 8 3 u 6 k 1 8 vxk 1 3 j ak j tk j   2 Sk  vyk 1 3 j bk j tk j   2 Sk  VX 2 7 k vxk  6  VY 2 7 k vyk  6  V VX2 VY2 V 0.434 V 2 1 R 0.366 i 1 8 Ni 1 2 Si 1 3 j ai j xi j bi j yi j ci j    85 Литература и методическое обеспечение 1. Полянин, А. Д. Справочник по линейным уравнениям матема- тической физики / А. Д. Полянин. – М. : ФИЗМАТЛИТ, 2001. – 576 с. 2. Методы конечных элементов : методические указания к лабо- раторным работам / сост. Н. А. Димитриева; Чувашский ун-т. – Че- боксары, 2004. – 32 с. 3. Тихонов, А. Н. Математическая физика / А. Н. Тихонов, А. А. Самарский. – М. : МГУ, 1999. – 799 с. 4. Зенкевич, О. Конечные элементы и аппроксимация / О. Зенке- вич, К. Морган. – М. : Мир, 1986. – 318 с. 5. Галлагер, Р. МКЭ. Основы / Р. Галлагер. – М. : Мир, 1984. – 428 с. 6. Зенкевич, О. Метод конечных элементов в технике / О. Зенке- вич. – М. : Мир, 1975. – 473 с. 7. Оден, Дж. Конечные элементы в нелинейной механике сплош- ных сред /Дж. Оден. – М. : Мир, 1976. – 464 с. 8. Чигарев, А. В. ANSYS для инженеров : справ. пособие / А. В. Чи- гарев, А. С. Кравчук, А. Ф. Смалюк. – М. : Машиностроение ; Ма- шиностроение-1, 2004. – 511 с. 9. Флетчер, К. Численные методы на основе метода Галеркина / К. Флетчер. – М. : Мир, 1988. – 352 с. 10. Годунов, С. К. Решение систем линейных уравнений / С. К. Го- дунов. – Новосибирск : Наука; Сиб. отд-ние, 1980. – 177 с. 11. Годунов, С. К. Разностные схемы: Введ. в теорию / С. К. Го- дунов, В. С. Рябенький. – М. : Наука, 1977. – 439 с. Программные пакеты и системы компьютерного моделирования Mathcad, Mathematica, Maple, MATLAB. 86 ПРИЛОЖЕНИЯ ПРИЛОЖЕНИЕ А Рекомендации по оформлению отчетов о лабораторных и практических работах 1. Отчет должен содержать следующие структурные части: ти- тульный лист (приложение Б); оглавление; цель работы; основную часть, разбитую на теоретическую и практическую главы; описание использованных методов, компьютерных программ, оборудования, материалов и т. д., а также сущность и основные результаты лабора- торного исследования; краткие выводы; список использованных ли- тературных и других источников; приложения (при необходимости). 2. Титульный лист отчета оформляется по форме согласно при- ложению Б. 3. В раздел «Приложения» включается вспомогательный матери- ал. В этот раздел включаются: промежуточные математические до- казательства, формулы и расчеты, оценки погрешности измерений; исходные тексты компьютерных программ и краткое их описание, таблицы и т. д. 4. Отчет печатается с использованием компьютера и принтера на одной стороне листа белой бумаги формата А4 (210  297 мм). При необходимости допускается представлять таблицы и иллюстрации на листах формата А3 (297  420 мм). 5. Набор текста осуществляется с использованием текстового ре- дактора Word. При этом рекомендуется использовать шрифты типа Times New Roman размером 13–14 пунктов, межстрочный интервал должен составлять 12–15 пунктов. Устанавливаются следующие раз- меры полей: верхнего и нижнего 20 мм, левого 30 мм, правого 10 мм. Шрифт печати должен быть прямым, светлого начертания, четким, черного цвета, одинаковым по всему объему текста отчета. Разре- шается использовать компьютерные возможности акцентирования внимания на определениях, терминах, теоремах, важных особенно- стях, применяя разное начертание шрифта. 6. Нумерация страниц дается арабскими цифрами. Первой стра- ницей отчета является титульный лист, который включают в общую нумерацию страниц. На титульном листе номер страницы не ставят, на последующих листах номер проставляют в центре нижней части листа без точки в конце. 87 7. Иллюстрации (скриншоты, фотографии, рисунки, чертежи, схе- мы, диаграммы, графики) и таблицы служат для наглядного пред- ставления характеристик объектов исследования, полученных теоре- тических и (или) данных численных экспериментов (расчетов). 8. Иллюстрации должны быть выполнены с помощью компью- терных технологий либо чернилами на белой непрозрачной бумаге. Качество иллюстраций должно обеспечивать возможность их чет- кого копирования. Допускается в качестве иллюстраций использо- вать распечатки с компьютерных программ, а также иллюстрации в цветном исполнении. 9. Иллюстрации, как правило, имеют наименование и поясни- тельные данные (подрисуночный текст), располагаемые по центру страницы. Пояснительные данные помещают под иллюстрацией, а со следующей строки – слово «Рис.», номер и наименование ил- люстрации, отделяя «точкой» номер от наименования. Точку в кон- це нумерации и наименований иллюстраций не ставят. 10. Цифровой материал отчета оформляют в виде таблиц. При оформлении таблиц необходимо руководствоваться следующими правилами: в таблице допускается применять шрифт на 1–2 пункта меньший, чем в тексте отчета; не следует включать в таблицу графу «Номер по порядку». 11. Формулы и уравнения в отчете (если их более одной) нуме- руют в пределах текста отчета. При оформлении формул и уравне- ний необходимо соблюдать следующие правила: формулы и урав- нения следует выделять из текста в отдельную строку. 12. Ссылки на источники в тексте отчета осуществляются путем приведения номера. Содержание сведений об источниках должно соответствовать ГОСТ 7.11–2004. 13. Приложения оформляют в конце отчета либо в виде отдель- ной части (книги), располагая их в порядке появления ссылок в тек- сте отчета. Каждое приложение следует начинать с нового листа с указанием в правом верхнем углу слова «ПРИЛОЖЕНИЕ», напеча- танного прописными буквами. Приложение должно иметь содержа- тельный заголовок, который размещается с новой строки по центру листа с прописной буквы. Приложения обозначают заглавными бук- вами русского алфавита, начиная с А, например: ПРИЛОЖЕНИЕ А, ПРИЛОЖЕНИЕ Б, ПРИЛОЖЕНИЕ В. 88 ПРИЛОЖЕНИЕ Б Пример оформления титульного листа отчета о лабораторных работах Министерство образования Республики Беларусь Белорусский национальный технический университет Машиностроительный факультет Кафедра «Теоретическая механика» ОТЧЕТ о лабораторной работе № 6 «Решение дифференциального уравнения в частных производных с помощью МКЭ. Численная реализация расчета краевой задачи для уравнения Лапласа» Выполнил: студент группы 103911 Иванов Иван Иванович дата, подпись Принял: преподаватель, ФИО дата, подпись Отметка о защите работы ________ дата, подпись Минск, 20__ 89 СОДЕРЖАНИЕ Предисловие ............................................................................................ 3 Лабораторная работа № 1 Построение физических, расчетных и математических моделей физических явлений и технических процессов. Методы решения определяющих уравнений математической модели. Метод конечных элементов. Основные этапы практической реализации ............................................................................................... 4 Лабораторная работа № 2 Определение перемещений и вычисление углов разворота произвольных элементов конструкций .............................................. 20 Лабораторная работа № 3 Метод Галеркина .................................................................................. 38 Лабораторная работа № 4 Решение задачи Коши методом конечных элементов ...................... 45 Лабораторная работа № 5 Численное решение краевой задачи методом конечных элементов ............................................................................................... 56 Лабораторная работа № 6 Решение дифференциального уравнения в частных производных с помощью МКЭ на примере решения краевой задачи для уравнения Лапласа ......................................................................... 67 Литература и методическое обеспечение ....................................... 85 ПРИЛОЖЕНИЯ .................................................................................. 86 Приложение А. Рекомендации по оформлению отчетов о лабораторных и практических работах ........................................... 86 Приложение Б. Пример оформления титульного листа отчета о лабораторных работах ....................................................................... 88 90 Учебное издание ШИРВЕЛЬ Павел Иванович ОСНОВЫ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ В МЕХАТРОНИКЕ Учебно-методическое пособие для студентов технических специальностей высших учебных заведений В 2 частях Часть 1 Редактор Т. Н. Микулик Компьютерная верстка Н. А. Школьниковой Подписано в печать 12.06.2015. Формат 6084 1/16. Бумага офсетная. Ризография. Усл. печ. л. 5,23. Уч.-изд. л. 4,09. Тираж 50. Заказ 440. Издатель и полиграфическое исполнение: Белорусский национальный технический университет. Свидетельство о государственной регистрации издателя, изготовителя, распространителя печатных изданий № 1/173 от 12.02.2014. Пр. Независимости, 65. 220013, г. Минск.