Размещено на http://www.allbest.ru/

СОДЕРЖАНИЕ

Введение

1. Численные методы моделирования

1.1 Понятие метода конечных элементов

1.2 Описание основных типов конечных элементов и построение сетки

1.3 Общий алгоритм метода конечных элементов

1.4 Обзор программных средств конечно-элементного анализа

2. Алгоритмический анализ задачи

2.1 Постановка задачи

2.2 Описание математической модели

2.3 Графическая схема алгоритма

3. Программная реализация поставленной задачи

3.1 Реализация модели конструкции в пакете ANSYS

3.2 Реализация на языке программирования C#

3.3 Исследование полученных результатов

3.4 Реализация интерфейса пользователя

Заключение

Список использованных источников

Приложение А Моделируемая конструкция

Приложение Б Исходный код программы на C#

Приложение В Реализация задачи в MathCad

Приложение Г Результаты моделирования в Ansys

Приложение Д Результаты разработанного приложения



ВВЕДЕНИЕ

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

Наибольшее распространение при анализе напряжений и деформаций в инженерных конструкциях получил метод конечных элементов (МКЭ), а также конечных разностей (МКР). Распространенной системой, которая использует метод конечных элементов, является ANSYS.

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

Целью курсовой работы является моделирование напряжённо-деформированного состояния конструкции в пакете ANSYS и написание программного приложения на языке высокого уровня, выполняющего аналогичный расчёт конструкции методом конечных элементов. В настоящее время существует достаточно много программных продуктов для решения отдельных классов задач, основанных на методе конечных элементов. Особенностью ANSYS является чрезвычайно широкий спектр задач, которые он в состоянии решать. Сюда входят задачи: расчетов на прочность (как линейные, так и нелинейные), теплообмена, гидродинамики, смешанные и даже акустики. Актуальность темы заключается в том, что расчёты подобного рода распространены в нынешнем веке и являются вполне востребованными при реализации сложных конструкций.



1. ЧИСЛЕННЫЕ МЕТОДЫ МОДЕЛИРОВАНИЯ

1.1 Понятие метода конечных элементов

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

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

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

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

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

Для этого требуется решить систему линейных алгебраических уравнений вида

KЫ= (1)

Здесь - вектор внешних сил, а есть узловой вектор начальных сил, который имеет место, например, при учете начальных температурных напряжений. Матрица K называется матрицей жесткости.

Порядок этой системы равен глобальному числу степеней свободы системы. По вычисленным таким образом перемещениям определяются напряжения и деформации. Физический смысл векторов U и F определяются областью применения МКЭ в соответствии с таблицей 1.1.



Таблица 1.1 - Физический смысл векторов Ы и F? в различных конечно-элементных приложениях

Область применения	Ы	F?
Механика деформируемого твердого тела	Перемещение	Сила
Теплоперенос	Теплопроводность	Тепловой поток
Гидромеханика	Скорость	Поток
Электростатика	Электрический потенциал	Плотность заряда
Магнитостатика	Магнитный потенциал	Интенсивность магнитного поля

Основные этапы решения задач с применением МКЭ могут быть представлены в виде схемы (рисунок 1.1).

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

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

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

Численное решение системы уравнений равновесия выполняется, как правило, автоматически с использованием ЭВМ.

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

Рисунок 1.1 - Основные этапы решения задачи с применением МКЭ

1.2 Описание основных типов конечных элементов и построение сетки

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

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

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

2) Узловые точки. Узловые точки, или узлы предназначены для описания геометрии элемента и задания физических степеней свободы (числа неизвестных функций). Как правило, узлы находятся в угловых или крайних точках элемента. Иногда вводят дополнительные узлы, расположенные внутри элемента. Число узлов связано с порядком аппроксимации, который обеспечивает данный конечный элемент. Элементы, имеющие только угловые узлы, называются линейными и обеспечивают линейную интерполяцию геометрии и функций. Элементы, имеющие дополнительные узлы на своих границах между угловыми точками, могут обеспечивать квадратичную или даже кубическую интерполяцию (рисунок 1.2). При наличии современных автоматических генераторов конечно-элементных сеток часто бывает проще и удобнее разбить конструкцию на большое число линейных элементов простой формы, чем использовать элементы высокого порядка, требующие для построения сетки значительной работы вручную. В то же время квадратичная и кубическая интерполяции обеспечивают более высокую точность расчета.

3) Геометрия элемента. Геометрия элемента определяется расположением узловых точек. Большинство элементов, используемых в расчетах, имеют простую геометрическую форму. Например, в одномерном случае элементы обычно представляют собой прямолинейные отрезки или сегменты кривых линий; в двумерном случае элементы имеют трех- или четырехстороннюю форму (см. рисунок 1.2); в трехмерных задачах наиболее распространены такие геометрические фигуры, как тетраэдры, гексаэдры (см. рисунок 1.2) и призмы.

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

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

6) Свойства сечения. К свойствам сечения относятся площади и моменты инерции одномерных конечных элементов (балки, стержни), а также толщины двумерных конечных элементов (пластины, оболочки).



Рисунок 1.2 - Конечные элементы для одно-, двух- и трехмерных задач с разным порядком аппроксимации.

моделирование пакет интерфейс сетка

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

- линейные элементы требуют более мелкой сетки, чем элементы более высокого порядка (с промежуточными узлами);

- упорядоченная сетка предпочтительнее произвольной;

- прямоугольная сетка с четырьмя узлами (рисунок 1.4, а) более выгодна, чем сетка с треугольными элементами (рисунок 1.4, б);

- сетка треугольных элементов с промежуточными узлами (рисунок 1.4, в) имеет точность, близкую к сетке прямоугольных элементов с четырьмя узлами (см. рисунок 1.4, а);

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

- аппроксимация смещений кубическим полиномом (элемент третьего порядка) не требует мелкой сетки.

Рисунок 1.3 - Примеры свободного и регулярного разбиения

Рисунок 1.4 - Различные типы конечно-элементных сеток

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

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

Рисунок 1.5 - Допустимые размеры элементов

Форма конечных элементов также влияет на точность вычислений. Следует избегать слишком вытянутых элементов, так как элементы с примерно одинаковыми размерами сторон дают меньшую ошибку (рисунок 1.6, а). Запрещается использовать четырехугольные элементы с углами, превышающими 180° (рисунок 1.6, б).

Рисунок 1.6 - Примеры форм элементов

1.3 Общий алгоритм метода конечных элементов

Последовательность процедур алгоритма МКЭ может быть представлена в следующем виде:

1. Дискретизация рассматриваемой области.

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

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

2. Выбор вариационного принципа.

Выбор вариационного принципа определяет основные неизвестные функции, через которые впоследствии устанавливаются остальные неизвестные. В задачах механики деформируемого твердого тела используются следующие вариационные принципы: принцип Лагранжа, в соответствии с которым варьируются перемещения; принцип Кастильяно (варьируются напряжения), принцип Рейсснера (варьируются перемещения и напряжения), принцип Ху-Вашицы (варьируются перемещения, напряжения и деформации). В практических расчетах чаще всего используется принцип Лагранжа.

3. Выбор аппроксимирующих функций.

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

- критерию полноты: при стремлении размеров элемента к нулю выбранные функции формы должны обеспечить любые простые значения;

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

При выполнении этих критериев с увеличением числа конечных элементов, моделирующих конструкцию, результаты расчета монотонно сходятся к точному решению.

4. Реализация вариационного принципа.

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

- методом непосредственного сложения жесткостей;

- методом конгруэнтного преобразования;

- при помощи конечно-разностных операторов.

5. Учет граничных условий.

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

6. Решение системы алгебраических уравнений.

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

7. Определение деформаций и напряжений.

После определения узловых перемещений в соответствии с известными соотношениями теории упругости могут быть определены деформации и напряжения. [3, стр. 25]



1.4 Обзор программных средств конечно-элементного анализа

Современные комплексы программ, в которых используется МКЭ, позволяют получать приближенные численные решения при расчете конструкций на статические и динамические нагрузки для широкого класса материалов с различными механическими характеристиками и поведением. Расчет конструкций на статические нагрузки может производиться с учетом физической и геометрической нелинейности, температурных полей, взаимодействия с другими средами. Производится расчет критических нагрузок, при которых конструкция или ее элементы теряют устойчивость, поведения конструкции после потери устойчивости. МКЭ позволяет также определить нагрузки, при которых происходит разрушение конструкции. Основными динамическими задачами являются: расчет собственных колебаний конструкции; динамический отклик на нагрузку, зависящую от времени; распространение волн. [2, стр. 56] Наибольшее распространение из программных средств конечно-элементного анализа получили: ABAQUS, ADINA, ASKA/DYNAN, ANSYS, MARC, MSC/NASTRAN, EUFEMI, COSMOS, HERCULE, MODULEF, SAP-7, LS-DYNA.

ANSYS

ANSYS - это универсальный конечноэлементный программный пакет (разработчиком которого является компания ANSYS Inc.), который позволяет решать в единой пользовательской среде (и, что очень важно, на одной и той же конечноэлементной модели) широкий круг задач в областях: прочности; тепла; гидрогазодинамики; электромагнетизма; междисциплинарного связанного анализа, объединяющего все четыре типа; оптимизации конструкции на основе всех вышеприведенных типов анализа.

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

Также существуют встраиваемые в CAD-системы расчетные модули ANSYS, позволяющие реализовать преимущества технологии "сквозного проектирования".

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

MSC.Nastran/MSC.Patran

Профессиональная CAE система применяемая при решении задач механики деформируемого твердого тела, теплообмена, гидро- и газодинамики, электро- и магнитостатики, а также других областей физики. MSC.Nastran обеспечивает полный набор расчетов, включая расчет напряженно - деформированного состояния, собственных частот и форм колебаний, анализ устойчивости, решение задач теплопередачи, исследование установившихся и неустановившихся процессов, акустических явлений, нелинейных статических процессов, нелинейных динамических переходных процессов, расчет критических частот и вибраций роторных машин, анализ частотных характеристик при воздействии случайных нагрузок, спектральный анализ и исследование аэроупругости. Предусмотрена возможность моделирования практически всех типов материалов, включая композитные и гиперупругие. Расширенные функции включают технологию суперэлементов (подконструкций), модальный синтез и макроязык DMAP для создания пользовательских приложений. Основу MSC.Nastran составляют отработанная технология элементов и надежные численные методы. Программа позволяет одновременно применять в одной и той же модели h- и p-элементы для достижения точности расчета при минимальных компьютерных ресурсах. Элементы супер высокого порядка аппроксимации - p-элементы - хорошо отражают криволинейную геометрию конструкции и обеспечивают высокую точность при детальном расчете напряжений. Эти элементы автоматически адаптируются к желаемому уровню точности. Численные методы разреженных матриц, используемые при любом типе расчетов, резко повышают скорость вычислений и минимизируют объем требуемой дисковой памяти, что повышает эффективность обработки данных.

LS-DYNA

Многоцелевой конечно-элементный расчетный комплекс для проведения динамических расчетов явным методом. Применяется для расчета высоконелинейных динамических переходных процессов. В отличие от широко известных программ ABAQUS, ANSYS, MSC.Nastran, LS-DYNA изначально создавалась для решения динамических нелинейных задач. Далее была создана фирма Livermore Software Technology Corporation (LSTC), которая стала разрабатывать коммерческие версии программы, ориентированные на автомобильную промышленность. В настоящее время в LS-DYNA реализованы эффективные методы решения перечисленных задач, в том числе явный и неявный методы интегрирования, многокомпонентная гидродинамика (Mulimaterial Eulerian Hydrodynamics), вычислительная гидродинамика несжимаемых потоков, бессеточный метод сглаженных частиц (SPH - Smoothed Particle Hydrodynamics), бессеточный метод, основанный на методе Галеркина (EFG - Element Free Galerkin method). В LS-DYNA встроены процедуры автоматического перестроения и сглаживания конечно-элементной сетки при вырождении элементов - произвольные лагранжево-эйлеровы сетки (Arbitrary Lagrangian-Euleran), высокоэффективные алгоритмы решения контактных задач, широкий набор моделей материалов, возможности пользовательского программирования, а также процедуры лагранжево-эйлерового связывания и расчета многокомпонентных течений сжимаемых сред на подвижных эйлеровых сетках.

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

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



2. АЛГОРИТМИЧЕСКИЙ АНАЛИЗ ЗАДАЧИ

2.1 Постановка задачи

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

Разработанное Windows-приложение, должно выполнять следующие действия:

- разбивать плоскую конструкцию на конечные элементы;

- графически изображать конструкцию до и после приложения нагрузки.

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

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



Рисунок 2 - Вид моделируемой конструкции.

2.2 Описание математической модели

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

Распределённая нагрузка действует в точках по осям OX и OY P=20000 Н/м. Материал, из которого сделана конструкция - сталь. Модуль упругости для данного материала равен 2•10¹¹ Па, а коэффициент Пуассона равен 0,27.

Данные с размерами детали представлены в таблице 2.1. Вид детали с размерами находятся в приложении А.

Таблица 2.1 - Геометрические характеристики детали

Ширина a, мм	Длина b, мм	Диаметр, D1, мм	Радиус, R1, мм
8	16	1	0.5

Перемещения каждого узла имеют два компонента, формула (2.1):

(2.1)



шесть компонентов перемещений узлов элемента образуют вектор перемещений {д}:

(2.2)

Перемещение любой точки внутри конечного элемента определяется соотношениями (2.3) и (2.4):

(2.3)

(2.4)

При объединении (2.3) и (2.4) в одно уравнение получается следующее соотношение:

(2.5)

Деформации и перемещения связаны между собой следующим образом:

(2.6)

При подстановке (2.5) в (2.6) получается соотношение (2.7):

(2.7)

Соотношение (2.7) можно представить в виде:

(2.8)

где [В] называется градиентная матрица вида (2.9):

(2.9)

Функции формы линейно зависят от координат x, y, и следовательно, градиентная матрица не зависит от координат точки внутри конечного элемента, и деформации и напряжения внутри конечного элемента в этом случае постоянны.

При плоском деформированном состоянии в изотропном материале матрица упругих постоянных [D] определяется по формуле (2.10):

(2.10)

где Е - модуль упругости, - коэффициент Пуассона.

Матрица жесткости конечного элемента имеет вид:



(2.11)

где h^e - толщина, А^e - площадь элемента.

Уравнение равновесия i -ого узла имеет вид:

(2.12)

Для учета условий закрепления существует следующий метод. Пусть имеется некоторая система N уравнений (2.13):

(2.13)

В случае, когда одна из опор неподвижна, т.е. U_i=0, используют следующую процедуру. Пусть U₂=0, тогда:

(2.14)

то есть соответствующие строка и столбец задаются нулевыми, а диагональный элемент - единичным. Соответственно, приравнивается нулю и F₂ [1].

Для решения полученной системы выбираем метод Гаусса.

Прямой ход: путём элементарных преобразований над строками систему приводят к ступенчатой или треугольной форме, либо устанавливают, что система несовместна. Выбирается разрешающая строка k-ая , где k = 0…n - 1, и для каждой следующей строки выполняется преобразование элементов

(2.15)

для i = k+1, k+2 … n-1; j = k+1,k+2 … n.

2. Обратный ход: осуществляется определение значений неизвестных. Из последнего уравнения преобразованной системы вычисляется значение переменной х_n, после этого из предпоследнего уравнения становится возможным определение переменной x_n-1 и так далее [3].

2.3 Графическая схема алгоритма

Графическая схема алгоритма приведена на рисунке 4. В блоке 1 происходит ввод исходных данных. Далее в блоке 2 идет построение конечно-элементной сетки. На основании введённых данных в блоке 3 и 4 строятся локальные и глобальная матрицы жесткости. Для получения системы алгебраических уравнений в блоке 4 формируется вектор свободных членов. В блоке 5 полученная система решается методом Гаусса. В блоке 6 происходит вывод результатов на экран и в файл на жёстком диске.



Рисунок 6 - Графическая схема алгоритма



3. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ ПОСТАВЛЕННОЙ ЗАДАЧИ

3.1 Реализация модели конструкции в пакете ANSYS

На первом шаге моделирования построим исследуемую деталь:

Задаем ключевые точки:

Main Menu > Preprocessor > Modeling > Keypoints > In Active CS...

1 (0;0), 2 (8;0), 3 (16.5), 4 (16;8), 5 (8;8), 6 (8;0)

Соединяем полученные точки линиями, смотреть рисунок 3.1:

Preprocessor > Modeling > Create > Lines > Lines > Straight Line...

Рисунок 3.1 - Соединение линиями

Создаем плоскость:

Preprocessor > Modeling > Create > Areas > Arbitrary > By Lines

Вырезание отверстия:

Preprocessor > Modeling > Create > Areas > Circle > Solid Circle

Выбираем тип конечного элемента:

Preprocessor > Element Type > Add/Edit/Delete... > Add...

Задаем характеристики материала, смотреть рисунок 3.2:

Preprocessor > Material props > Materials Modes… > Structural > Linear > Elastic > Isotropic.

Рисунок 3.2 - Ввод модуля Юнга и коэффициента Пуассона

Preprocessor > Material props > Materials Modes… > Structural > Density.

Рисунок 3.3 - Ввод плотности материала

Построение конечно-элементной сетки с треугольными элементами смотреть рисунок 3.4:

Preprocessor > Meshing > MeshTool



Рисунок 3.4 - Конечно-элементная сетка

Закрепляем и нагружаем деталь в соответствии с заданием смотреть рисунок 3.5:

Sulution > Define Loads > Apply > Structural > Displacement > On Lines

Sulution > Define Loads > Apply > Structural > Force/Moment > On Nodes



Рисунок 3.5 - Нагрузка и закрепление детали

Рассчитываем деталь:

Sulution > Solve > Current LS

Вывод результатов (см.приложение Г):

General Postproc > Plot Results > Deformed Shape

General Postproc > Plot Results > Contour Plot > Nodal Solu > DOF Solution…

General Postproc > List Results > Nodal Solution > DOF Solution…

3.2 Реализация задачи на языке программирования C#

Для реализации математической модели был выбран язык программирования C#. Листинг программы приведен в приложении Б.

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



Таблица 3.1. Структура класса Node

Имя переменной	Тип	Комментарий
_x	double	Координата узла по оси OX
_y	double	Координата узла по оси OY
_Fx	double	Горизонтальная составляющая силы в узле
_Fy	double	Вертикальная составляющая силы в узле
_Px	double	Признак перемещения по оси OX
_Py	double	Признак перемещения по оси OY
_u	double	Смещение по оси OX
_v	double	Смещение по оси OY

Таблица 3.2. Структура класса Element

Имя переменной	Тип	Комментарий
_density	double	плотность материала элемента
_elasticity	double	модуль Юнга
_factor_Puason	double	коэффициент Пуассона

Для решения СЛАУ был реализован метод Гаусса.

В программе реализованы следующие варианты вывода результатов расчета (приложение Д):

Вывод результата по каждому узлу, рисунок Д.2. При наведении курсора мыши на узел на всплывающей панели отображаются перемещение узла по осям OX, OY, а также силы по оси OX и OY.

Графический вывод деформации детали рисунок Д.3.

Графический вывод перемещений в узлах по оси OX, рисунок Д.4.

Графический вывод перемещений в узлах по оси OY, рисунок Д.5.

Графический вывод суммы перемещений в узлах по оси OX и OY, рисунок Д.6.

Для проверки корректности работы программы был разработан метод предоставления исходных данных для пакета MathCad. Реализация в MathCad приведена в приложении В.



3.3 Исследование полученных результатов

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

Для проверки результатов, поставленная задача была реализована в пакете ANSYS и MathCad. Для сравнения были выбраны основные точки рассчитываемой конструкции, рисунок 13.

№ узла	ANSYS	Разработанное приложение	MathCad
	UX,m	UY,m	UX,m	UY,m	UX,m	UY,m
1	0	0	0	0	0	0
2	0.14283E-06	-0.11956E-06	0.10485E-06	-0.09956E-06	0.10485E-06	-0.09956E-06
3	0.16194E-05	0.15692E-05	0.12295E-05	0.13691E-05	0.12295E-05	0.13691E-05
4	0.77953E-06	0.10403E-04	0.57957E-06	0.08453E-04	0.57957E-06	0.08453E-04
5	-0.20460E-05	0.10833E-04	-0.20434E-05	0.08833E-04	-0.20434E-05	0.08833E-04
6	0	0	0	0	0	0

Анализ результатов показывает, что рассчитанные величины практически одинаковы, однако перемещения наблюдаются в некоторых точках, однако их величина незначительна и не влияет на прочность конструкции. Процент сходимости результатов разработанной программы и ANSYS составил порядка 80%. А результаты c Mathcad сошлись на 95%.

3.4 Реализация интерфейса пользователя

Программа моделирования напряженно-деформированного состояния детали формируется на языке высоко уровня С#. Интерфейс программы очень простой и интуитивно понятен.

Вначале вводится сила, которая будет воздействовать на деталь.

Рисунок 3.6 - Отображение формы

При нажатии на «Сетка» наносится сетка конечных элементов и выводится на экран еще 2 кнопки: «Решить» и «MathCad». При нажатии на «MathCad» мы получаем решение в системе MathCad, а при нажатии на «Решить» производится расчет в С#. Графическая оболочка представлена на рисунке 3.7:



Рисунок 3.7 - Нанесение сетки

Появляется меню, где отображаются основные действия работы над деталью. Результат выполнения представлен на рисунке 3.8:

Рисунок 3.8 - Операции с деталью



ЗАКЛЮЧЕНИЕ

В ходе выполнения курсового проекта были определены размеры деталей соединяемых с помощью заклёпочных соединений методом конечных элементов в пакете моделирования ANSYS с помощью разработанного программного обеспечения C#.

На начальном этапе конструкция была смоделирована в пакете ANSYS. Далее был произведен расчет с помощью системы MathCad. Результаты расчётов ANSYS и MathCad практически полностью совпадают.

Кроме этого было разработано приложение реализующее расчёт модели с применением метода конечных элементов. Полученные результаты практически полностью совпадают с результатами ANSYS и MathCad. Допустимые погрешности находятся в пределах нормы.

Приложение позволяет довольно быстро рассчитывать размеры деталей соединяемых с помощью заклёпочных соединений. Данный продукт вполне может использоваться при расчёте более сложных задач.



СПИСОК ИСПОЛЬЗОВАНЫХ ИСТОЧНИКОВ

Каплун, А. ANSYS в руках инженера: Практическое руководство / А. Каплун - М.: Едиториал УРСС, 2003. - 272 с.

Басов, К. ANSYS в примерах и задачах / К. Басов - М.: КомпьютерПресс, 2002. - 224 с.

Шалумов, А. Введение в ANSYS / А. Шалумов - Москва, 2002. 437 с.

Зенкевич, О. Конечные элементы и аппроксимация / О. Зенкевич - М.: Мир, 1986. 318 с.

Баженов, В. Численные методы в механике / В. Баженов - М.: КомпьютерПресс, 2004. 356с.

Молчанов, И. Основы метода конечных элементов / И. Молчанов - Киев: Наукова Думка, 1989. 272 с.



Приложение Б

Исходный код программы на С#

using System;

using System.Collections;

using System.ComponentModel;

using System.Data;

using System.Drawing;

using System.Drawing.Drawing2D;

using System.Linq;

using System.Text;

using System.Windows.Forms;

using System.IO;

namespace WindowsFormsApplication2

{

public partial class Form1 : Form

{

int select_node = -1; // номер узла над которым висит фокус мыши для отображения свойств (-1 если не один не выбран)

Color Color_background = Color.Black;

// для буферизации

private BufferedGraphicsContext context;

private BufferedGraphics grafx;

double max_size = 600; //максимальныая высота либо ширина детали

StreamReader read_files;

ArrayList Nodes = new ArrayList(0);

ArrayList Elements = new ArrayList(0);

double scale = 18000; // переменная для хранения масштаба (сколько пикселей на одну единицу размера детали)

int dx = 30; // смещение детали по Ox

int dy = -150; // смещение детали по Oy

ArrayList kol_given_displacement = new ArrayList(0);

bool plot_result = false;

Color[] color_plot; // массив содержащий цвет каждого элемента в соответствии со смещением

public Form1()

{

InitializeComponent();

}

private void button1_Click(object sender, EventArgs e)

{

// установка начальных значений (если мы разбиваем еще раз)

Nodes = new ArrayList(0);

Elements = new ArrayList(0);

plot_result = false;

// чтение файла с элементами

string str;

read_files = new StreamReader("ELIST.lis");

while ((str = read_files.ReadLine()) != null)

{

int l = 0;

int a = 0;

int b = 0;

int c = 0;

str = str.Trim();

string[] s = String.Join(",", str.Split('.')).Split();

for (int i = 0; i < s.Length; i++)

{

try

{

if (s[i] != "")

{

if (l == 8)

{

a = Convert.ToInt32(s[i]);

l++;

}

if (l == 7)

{

b = Convert.ToInt32(s[i]);

l++;

}

if (l == 6)

{

c = Convert.ToInt32(s[i]);

l++;

}

if (l < 6)

{

Convert.ToInt32(s[i]);

l++;

}

}

}

catch { }

}

if (l == 9)

{

Elements.Add(new Element(a, b, c, 7800, 2E+11, 0.27));

}

}

read_files.Close();

// Чтение файла с узлами

int number_node=1;

read_files = new StreamReader("NLIST.lis");

while ((str = read_files.ReadLine()) != null)

{

int l = 0;

double x=0;

double y=0;

str=str.Trim();

string[] s = String.Join(",", str.Split('.')).Split();

for (int i = 0; i < s.Length; i++)

{

try

{

if (s[i] != "")

{

if (l == 2)

{

y = Convert.ToDouble(s[i]);

l++;

}

if (l == 1)

{

x = Convert.ToDouble(s[i]);

l++;

}

if (l == 0)

{

Convert.ToInt32(s[i]);

l++;

}

}

}

catch { }

}

if (l == 3)

{

// проверяем на участие конкретного узла в формировании конечно элементной сетки

bool t = false; // переменная следящая за наличием узла в элементах

for (int i = 0; i < Elements.Count; i++)

{

Element temp = (Element)Elements[i];

if (temp.A == number_node)

{

temp.A = Nodes.Count + 1;

t = true;

}

if (temp.B == number_node)

{

temp.B = Nodes.Count + 1;

t = true;

}

if (temp.C == number_node)

{

temp.C = Nodes.Count + 1;

t = true;

}

}

if(t)

Nodes.Add(new Node(x, y, 0, 0, 0, 0, 0, 0));

number_node++;

}

}

read_files.Close();

first_value();

button3.Visible = true;

button2.Visible = true;

grafx.Graphics.Clear(Color_background);

DrawToBuffer(grafx.Graphics);

this.Refresh();

}

public void first_value()

{

int j = 0;

for (int i = 0; i < Nodes.Count; i++)

{

Node temp=(Node) Nodes[i];

if (Math.Pow(temp.x-2,2)+ Math.Pow(temp.y-2,2)<=0.26)

{

temp.Px = 1;

temp.Py = 1;

}

if (Math.Pow(temp.x - 6, 2) + Math.Pow(temp.y - 2, 2) <= 0.26)

{

temp.Px = 1;

temp.Py = 1;

}

if (Math.Pow(temp.x - 4, 2) + Math.Pow(temp.y - 6, 2) <= 0.26)

{

temp.Px = 1;

temp.Py = 1;

}

if ((temp.x >= 13 && temp.y == 8))

{

temp.Fx = -Convert.ToDouble(textBox1.Text);

j++;

}

if ((temp.x == 16 && temp.y > 5))

{

temp.Fx = -Convert.ToDouble(textBox1.Text);

j++;

}

Nodes[i] = temp;

}

}

private void panel1_Paint(object sender, PaintEventArgs e)

{

}

// расчет площади одного элемента

public double area_element(Element element)

{

Node A = (Node)Nodes[element.A-1];

Node B = (Node)Nodes[element.B-1];

Node C = (Node)Nodes[element.C-1];

return Math.Abs((((A.x - B.x) * (A.y - C.y)) - (A.x - C.x) * (A.y - B.y)) / 2);

}

public void DrawToBuffer(Graphics g)// создание рисунка в буфер

{

Brush plot_result_deformet = Brushes.Blue; // цвет фона на выводе деформации детали

Brush background = Brushes.YellowGreen; // цвет исходной детали (после загрузки формы, до разбиения)

Brush backgroubd_work = Brushes.Black; // цвет заливки рабочей области

Pen element_lines = new Pen(Brushes.White, 1);// цвет и размер конечно элементной сетки

Brush Node_point = Brushes.White;// цвет не нагруженных узлов

Brush Node_point_Force = Brushes.Red; // цвет нагруженных узлов

Brush Node_point_fixing = Brushes.Blue; // цвет жестко закрепленных узлов

Pen system_axes = new Pen(Brushes.Gray, 1);

Brush property = Brushes.LightYellow;

if (Elements.Count > 0 && Nodes.Count > 0)

{

background = Brushes.Green;

}

if (!plot_result)

{

GraphicsPath contour = new GraphicsPath();

contour.AddLines(new Point[]

{

/*1*/ new Point(Convert.ToInt32(scale * 0)+dx, panel1.Height-Convert.ToInt32(scale * 0)+dy),

/*2*/ new Point(Convert.ToInt32(scale * 8)+dx, panel1.Height-Convert.ToInt32(scale * 0)+dy),

/*3*/ new Point(Convert.ToInt32(scale * 16)+dx, panel1.Height-Convert.ToInt32(scale * 5)+dy),

/*32*/ new Point(Convert.ToInt32(scale * 16)+dx, panel1.Height-Convert.ToInt32(scale * 8)+ dy),

/*4*/ new Point(Convert.ToInt32(scale * 0)+dx, panel1.Height-Convert.ToInt32(scale * 8)+ dy)

});

GraphicsPath contour2 = new GraphicsPath();

Region r = new Region(contour);

g.FillRegion(background, r);

g.FillEllipse(backgroubd_work, new Rectangle(Convert.ToInt32(scale * 1.6) + dx, panel1.Height - Convert.ToInt32(scale * 2.4) + dy, Convert.ToInt32(scale * 0.8), Convert.ToInt32(scale * 0.8)));

g.FillEllipse(backgroubd_work, new Rectangle(Convert.ToInt32(scale * 5.6) + dx, panel1.Height - Convert.ToInt32(scale * 2.4) + dy, Convert.ToInt32(scale * 0.8), Convert.ToInt32(scale * 0.8)));

g.FillEllipse(backgroubd_work, new Rectangle(Convert.ToInt32(scale * 3.6) + dx, panel1.Height - Convert.ToInt32(scale * 6.4) + dy, Convert.ToInt32(scale * 0.8), Convert.ToInt32(scale * 0.8)));

}

int x = panel1.Width;

int y = panel1.Height;

g.DrawLine(system_axes, new Point(x - 60, y - 20), new Point(x - 60, y - 50));

g.DrawLine(system_axes, new Point(x - 60, y - 50), new Point(x - 55, y - 40));

g.DrawLine(system_axes, new Point(x - 60, y - 50), new Point(x - 65, y - 40));

g.DrawString("OY", Font, Brushes.Red, new PointF(x - 70, y - 65));

g.DrawLine(system_axes, new Point(x - 30, y - 20), new Point(x - 60, y - 20));

g.DrawLine(system_axes, new Point(x - 30, y - 20), new Point(x - 40, y - 15));

g.DrawLine(system_axes, new Point(x - 30, y - 20), new Point(x - 40, y - 25));

g.DrawString("OX", Font, Brushes.Red, new PointF(x - 25, y - 27));

if (Elements.Count > 0 && Nodes.Count > 0)

{

for (int i = 0; i < Elements.Count; i++)

{

Element temp = (Element)Elements[i];

Node A_ = (Node)Nodes[temp.A - 1];

Node B_ = (Node)Nodes[temp.B - 1];

Node C_ = (Node)Nodes[temp.C - 1];

Point A = new Point(Convert.ToInt32(scale * (A_.x)) + dx, panel1.Height - Convert.ToInt32(scale * (A_.y)) + dy);

Point B = new Point(Convert.ToInt32(scale * (B_.x)) + dx, panel1.Height - Convert.ToInt32(scale * (B_.y)) + dy);

Point C = new Point(Convert.ToInt32(scale * (C_.x)) + dx, panel1.Height - Convert.ToInt32(scale * (C_.y)) + dy);

if (!plot_result || (plot_result && radioButton1.Checked))

g.DrawLines(element_lines, new Point[] { A, B, C, A });

if (!plot_result)

{

if (A_.Fx != 0 || A_.Fy != 0)

g.FillEllipse(Node_point_Force, A.X - 3, A.Y - 3, 6, 6);

else

if (A_.Px == 1 || A_.Py == 1)

g.FillEllipse(Node_point_fixing, A.X - 3, A.Y - 3, 6, 6);

else

g.FillEllipse(Node_point, A.X - 3, A.Y - 3, 6, 6);

if (C_.Fx != 0 || C_.Fy != 0)

g.FillEllipse(Node_point_Force, C.X - 3, C.Y - 3, 6, 6);

else

if (C_.Px == 1 || C_.Py == 1)

g.FillEllipse(Node_point_fixing, C.X - 3, C.Y - 3, 6, 6);

else

g.FillEllipse(Node_point, C.X - 3, C.Y - 3, 6, 6);

if (B_.Fx != 0 || B_.Fy != 0)

g.FillEllipse(Node_point_Force, B.X - 3, B.Y - 3, 6, 6);

else

if (B_.Px == 1 || B_.Py == 1)

g.FillEllipse(Node_point_fixing, B.X - 3, B.Y - 3, 6, 6);

else

g.FillEllipse(Node_point, B.X - 3, B.Y - 3, 6, 6);

}

}

if (plot_result)

{

double max_x=-1;

double max_dx=-1;

for(int i=0; i<Nodes.Count;i++)

{

Node temp=(Node) Nodes[i];

if(Math.Abs(temp.u)>max_dx)

{

max_dx=Math.Abs(temp.u);

max_x=temp.x;

}

if(Math.Abs(temp.v)>max_dx)

{

max_dx=Math.Abs(temp.v);

max_x=temp.y;

}

}

int scale_result =Convert.ToInt32(max_x/(max_dx*10));//

for (int i = 0; i < Elements.Count; i++)

{

Element temp = (Element)Elements[i];

Node A_ = (Node)Nodes[temp.A - 1];

Node B_ = (Node)Nodes[temp.B - 1];

Node C_ = (Node)Nodes[temp.C - 1];

Point _A = new Point(Convert.ToInt32(scale * (A_.x + A_.u * scale_result)) + dx, panel1.Height - Convert.ToInt32(scale * (A_.y + A_.v * scale_result)) + dy);

Point _B = new Point(Convert.ToInt32(scale * (B_.x + B_.u * scale_result)) + dx, panel1.Height - Convert.ToInt32(scale * (B_.y + B_.v * scale_result)) + dy);

Point _C = new Point(Convert.ToInt32(scale * (C_.x + C_.u * scale_result)) + dx, panel1.Height - Convert.ToInt32(scale * (C_.y + C_.v * scale_result)) + dy);

if (radioButton1.Checked)

{

g.FillPolygon(plot_result_deformet, new Point[] { _A, _B, _C, _A });

g.DrawLines(element_lines, new Point[] { _A, _B, _C, _A });

}

if (radioButton2.Checked || radioButton3.Checked|| radioButton4.Checked)

{

GraphicsPath elem = new GraphicsPath();

elem.AddLines(new Point[] { _A, _B, _C, _A });

Region r = new Region(elem);

g.FillRegion(new SolidBrush(color_plot[i]), r);

}

}

}

}

}

private void Form1_Load(object sender, EventArgs e)

{

this.SetStyle(ControlStyles.AllPaintingInWmPaint | ControlStyles.UserPaint, true);

this.SetStyle(ControlStyles.Opaque, false);

context = BufferedGraphicsManager.Current;

grafx = context.Allocate(this.CreateGraphics(), new Rectangle(0, 0,panel1.Width,panel1.Height));

grafx.Graphics.SmoothingMode = SmoothingMode.HighQuality;

grafx.Graphics.PixelOffsetMode = PixelOffsetMode.HighQuality;

grafx.Graphics.Clear(Color_background);

scale /= max_size;

DrawToBuffer(grafx.Graphics);

this.Refresh();

}

private void panel1_MouseMove(object sender, MouseEventArgs e)

{

}

private void Form1_Paint(object sender, PaintEventArgs e)

{

grafx.Render(e.Graphics);

}

private void Form1_MouseMove(object sender, MouseEventArgs e)

{

int x = e.X;

int y = e.Y;

select_node = -1;

for (int i = 0; i < Nodes.Count; i++)

{

Node temp = (Node)Nodes[i];

if (x > (Convert.ToInt32(scale * temp.x) + dx - 5) && x < (Convert.ToInt32(scale * temp.x) + dx + 5) && y < panel1.Height - Convert.ToInt32(scale * temp.y) + dy + 5 && y > panel1.Height - Convert.ToInt32(scale * temp.y) + dy - 5&&!plot_result)

{

select_node = i;

if (panel1.Height - Convert.ToInt32(scale * temp.y) + dy + 120 < panel1.Height)

panel2.Location = new Point(Convert.ToInt32(scale * temp.x) + dx, panel1.Height - Convert.ToInt32(scale * temp.y) + dy);

else

panel2.Location = new Point(Convert.ToInt32(scale * temp.x) + dx, panel1.Height - Convert.ToInt32(scale * temp.y) + dy-120);

label1.Text = (i+1).ToString();

label6.Text = temp.Fx.ToString("E");

label7.Text = temp.Fy.ToString("E");

label8.Text = temp.u.ToString("E");

label9.Text = temp.v.ToString("E");

panel2.Visible = true;

i = Nodes.Count;

}

}

if (select_node == -1)

panel2.Visible = false;

}

private void label3_Click(object sender, EventArgs e)

{

}

private void button2_Click(object sender, EventArgs e)

{

kol_given_displacement = new ArrayList(0);

// формирование массива массовых сил для каждого из узлов

double[] FM = new double[Nodes.Count];

for (int i = 0; i < Elements.Count; i++)

{

Element temp = (Element)Elements[i];

double area= area_element(temp);

FM[temp.A - 1] -= (area * temp.density * 9.81) / 3;

FM[temp.B - 1] -= (area * temp.density * 9.81) / 3;

FM[temp.C - 1] -= (area * temp.density * 9.81) / 3;

}

//*********************************************************

// Узловые силы, перемещения и их признаки

//размещаются в векторах последовательными

//парами значений: на четных позициях вертикальные

//компоненты, на нечетных-- горизонтальные:

double[] F_c = new double[Nodes.Count*2];

double[] Delta_c = new double[Nodes.Count*2];

double[] S = new double[Nodes.Count*2];

for (int i = 0; i < Nodes.Count; i++)

{

Node temp=(Node) Nodes[i];

F_c[2*i]=temp.Fx;

F_c[2 * i + 1] = (temp.Fy);

Delta_c[2 * i] = temp.u;

Delta_c[2 * i+1] = temp.v;

if (temp.Px == 1)

kol_given_displacement.Add(2*i);

if (temp.Py == 1)

kol_given_displacement.Add(2 * i+1);

S[2 * i] = temp.Px;

S[2 * i + 1] = temp.Py;

}

//*********************************************************

// матрица жесткости системы

double[,] M0 = new double[Nodes.Count*2, Nodes.Count*2];

double[,] MJS = M0;

for (int n = 0; n < Elements.Count; n++)

{

Element temp=(Element) Elements[n];

int[] V = new int[6];

V[0] = 2 * (temp.A-1);

V[1] = 2 * (temp.A-1)+1;

V[2] = 2 * (temp.B - 1);

V[3] = 2 * (temp.B - 1) + 1;

V[4] = 2 * (temp.C - 1);

V[5] = 2 * (temp.C - 1) + 1;

double[,] D = new double[3, 3];

double[,] B = new double[3, 6];

double[,] B_t = new double[6, 3];

double area = area_element(temp);

matrix_D(temp.elasticity, temp.factor_Puason, D);

matrix_B(temp, B);

for (int i = 0; i < 3; i++)

for (int j = 0; j < 6; j++)

B_t[j, i] = B[i, j]*area;

double[,] D_B_t = new double[6, 3];

for (int i = 0; i < 6; i++)

for (int j = 0; j < 3; j++)

{

D_B_t[i, j] = D[0, j] * B_t[i, 0] + D[1, j] * B_t[i, 1] + D[2, j] * B_t[i, 2];

}

double[,] K_s = new double[6, 6];

for (int i = 0; i < 6; i++)

for (int j = 0; j < 6; j++)

{

K_s[i, j] = D_B_t[i, 0] * B[0, j] + D_B_t[i, 1] * B[1, j] + D_B_t[i, 2] * B[2, j];

}

for (int i = 0; i < 6; i++)

for (int j = 0; j < 6; j++)

{

MJS[V[i], V[j]] = MJS[V[i], V[j]] + K_s[i, j];

}

}

//************************************************************

// Формируем вектор сил

double[] F1 = new double[Nodes.Count*2];

for (int i = 0; i < Nodes.Count * 2; i++)

{

double sum = 0;

for (int j = 0; j < Nodes.Count * 2; j++)

sum += MJS[i, j] * Delta_c[j] * S[j];

F1[i] = F_c[i] - sum;

}

//************************************************************

//удаление строк и столбцов матрицы жесткости системы

//для которой задано перемещение и удаление элементов из

// вектора сил

double[,] K = new double[Nodes.Count * 2 - kol_given_displacement.Count, Nodes.Count * 2 - kol_given_displacement.Count];

double[] F = new double[Nodes.Count * 2 - kol_given_displacement.Count];

int ii=0;

int jj=0;

for (int i = 0; i < Nodes.Count * 2; i++)

{

bool t=true;

for (int k = 0; k < kol_given_displacement.Count; k++)

if (((int) kol_given_displacement[k]) == i)

{

t = false;

k = kol_given_displacement.Count;

}

if (t)

{

F[ii] = F1[i];

jj = 0;

for (int j = 0; j < Nodes.Count * 2; j++)

{

t = true;

for (int k = 0; k < kol_given_displacement.Count; k++)

if (((int) kol_given_displacement[k]) == j)

{

t = false;

k = kol_given_displacement.Count;

}

if (t)

{

K[ii, jj] = MJS[i, j];

jj++;

}

}

ii++;

}

}

//************************************************************

//решение СЛАУ методом гаусса

int razmer=Nodes.Count*2-kol_given_displacement.Count;

for (int i = 0; i < razmer-1; i++)

for(int j=i+1; j<razmer;j++)

{

double f = K[j, i] / K[i, i];

for (int k = i + 1; k < razmer; k++)

{

K[j, k] = K[j, k] - f * K[i, k];

}

f = K[j, i] / K[i, i];

F[j] = F[j] - f * F[i];

K[j, i] = 0;

}

double[] delta1 = new double[razmer];

delta1[razmer - 1] = F[razmer - 1] / K[razmer - 1, razmer - 1];

for (int i = razmer - 2; i >= 0; i--)

{

double summa = F[i];

for (int j = i + 1; j < razmer; j++)

summa = summa - delta1[j] * K[i, j];

delta1[i] = summa / K[i, i];

}

double[] delta = new double[Nodes.Count * 2];

int index=0;

for (int i = 0; i < Nodes.Count * 2; i++)