МОДЕЛИРОВАНИЕ КОНСТРУКЦИЙ В ANSYS

Линейный плоский треугольный элемент

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

Этот элемент имеет ряд отличительных особенностей:

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

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

Точность решения, обеспечиваемая данным элементом, не может быть повышена путем добавления внутренних степеней свободы.

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

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

ЗАДАЧИ ТРЕХМЕРНОГО НАПРЯЖЕННОГО И ДЕФОРМИРОВАННОГО СОСТОЯНИЯ

ТРЕХМЕРНОЕ НАПРЯЖЕННОЕ И ДЕФОРМИРОВАННОЕ СОСТОЯНИЕ НА ОСНОВЕ МКЭ - Часть 1

1. Напряжённое состояние в точке

Вспомним, как нами было введено понятие напряжения. Рассмотрим тело, находящееся под действием системы уравновешенных сил (рис.48).

 Рис.1

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

Рис. 2

 Выделим малую площадку , содержащую точку . Внешнюю нормаль этой площадки обозначим .

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

 Величина  1 зависит от размеров площадки. Чтобы избавиться от влияния размеров площадки deltaA , перейдём к пределу и будем стягивать площадку к точке B

  

                                                                                  (1)

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

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

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

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

Докажем это важнейшее положение.

Нас интересует напряжённое состояние в точке . Выделим в окрестности этой точки малый прямоугольный параллелепипед (рис.50). Размеры параллелепипеда настолько малы, что напряжённое состояние в пределах параллелепипеда можно считать однородным, и что грани параллелепипеда являются площадками, проходящими через точку B и имеющими внешними нормалями оси X,Y,Z

 Рис. 3

Полное напряжение можно разложить на три составляющие, направленные по координатным осям. Всего будем иметь 9 компонент напряжённого состояния: три нормальных и шесть касательных напряжений. Нормальные напряжения обозначим sizma и припишем индекс, указывающий внешнюю нормаль.

Например: sizma— нормальное напряжение по площадке с внешней нормалью X.

Касательные напряжения обозначаются Tau с двумя индексами. Первый индекс указывает площадку, второй – направление напряжения.

Например Tau— касательное напряжение по площадке с внешней нормалью X, параллельное оси Y.

Нормальное напряжение считается положительным, если оно направлено по внешней нормали, т.е. является растягивающим.

Касательные напряжения считаются положительными, если при положительной внешней нормали они направлены в сторону положительных координатных осей.

Совершенно очевидно, что по противоположным граням параллелепипеда действуют равные по величине и противоположные по направлению напряжения.

Заметим также, что хотя на рис.50 компоненты напряжённого состояния показаны в виде векторов, они являются величинами скалярными.

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

Применение метода Аппроксимации в оптимизации конструкций

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

В Maple встроенная процедура pade(f, x = a, [n, m]) позволяет найти аппроксимацию Паде типа (n, m) в точке x = a функции f(x). Если при обращении к процедуре pade запросить доступную для пользователя информацию, указав, например, значение параметра infolevel[pade]:= 1, то будет получено сообщение об используемом для вычисления аппроксимации алгоритме. В случае, когда коэффициенты ряда Тейлора функции f(x) не содержат параметров и чисел с плавающей точкой, используется ≪быстрый≫ алгоритм Кабая и Чоя[2]. В остальных случаях нахождение аппроксимации проводится по алгоритму Геддеса [3],использующему символьные вычисления.

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

Пример 1

Найдем аппроксимацию Паде π2,3(f1) типа (2, 3) функции f1(x) = x+1,0001 (x+1,999)(x−2,001) в точке x = 0. В этом случае ядро матрицы (1) неодномерно, однако, как видно из примера, найденные знаменатель и числитель не содержат ≪лишних≫ нулей.

restart;

with(numapprox);

padef1:= pade

x+1,0001

(x+1,999)(x−2,001), x = 0, [2, 3];

−0,2500250626−0,2500000625x1,000000000+0,0005000002006x−0,2500000625x2

factor(denom(padef1), complex); solve(numer(padef1) = 0);

0,2500000625(x + 1,999000000)(x − 2,001000000)−1,000100000.

 Пример 2

Найдем аппроксимацию Паде π4,5(f2) типа (4, 5) функции f2(x) = (x−3,001)(x+1,9999) (x2+1)(x+4,0001) вточке x = 0. В этом случае ядро матрицы (1) неодномерно, и найденные знаменатель и числитель аппроксимации Паде имеют ≪лишние≫ нули (они выделены жирным шрифтом).

padef2:= pade

(x−3,001)(x+1,9999)

(x2+1)(x+4,0001) , x = 0, [4, 5];

−1,500387465−0,3753104091x−0,4121913908x

2−0,08614086896x3+0,1068576988x41,000000000+0,3333333332x+1,448275862x2+0,4401910336x3+0,4482758628x4+0,1068577004x5

factor(denom(padef2), complex); factor(numer(padef2), complex);

0,1068577004(x + 4,000100004)(x + 0, 09748654036 + 1, 526433054I) (x + 0, 09748654036 − 1, 526433054I)(x + 1,000000001I)(x − 1,000000001I) 0,1068576988(x + 1,999900006)(x + 0, 09748653975 + 1, 526433060I) (x + 0, 09748653975 − 1, 526433060I)(x − 3,001000018).

Отметим, что в обоих приведенных примерах мы находимся в так называемом сингулярном блоке таблицы Паде для аппроксимируемой рациональной функции и, согласно теории аппроксимаций Паде, π2,3(f1), π4,5(f2) должны были оказаться равными рациональным дробям f1(x) и f2(x).

Найденная с помощью Maple аппроксимация π2,3(f1) оказалась равной аппроксимируемой функции, и в этом случае не появилось лишних нулей числителя и знаменателя. Однако аппроксимация π4,5(f2) не совпала с f2(x), у числителя и знаменателя аппроксимации появились близкие нули – дуплеты Фруассара. Дело в том, что из неодномерного ядра матрицы Tn+1 (1) в первом случае (случайно) был выбран ≪правильный≫ знаменатель. Во втором же примере взятый знаменатель не совпал со знаменателем рациональной дроби f2(x).

Похожая картина наблюдается и при попытке найти π2,3(f1), π4,5(f2) с помощью системы Mathematica. Как показывают приведенные ниже примеры, при этом опять появляются дуплеты Фруассара.

Пример 3

Вновь, как в примере 1, найдем аппроксимацию Паде π2,3(f1) типа (2, 3) функции f1(x) = x+1,0001(x+1,999)(x−2,001) в точке x = 0. Появляющиеся при этом дуплеты Фруассара выделены жирным шрифтом.

In[1]:= PadeApproximant[

x+1,0001

(x+1,999)(x−2,001), {x, 0, {2, 3}}]

Out[1]= −0,250025−0,222305x+0,0276927x 2

1,000000000000000−0,110271x−0,250055x2+0,0276927x3

In[2]:= Roots[Denominator[Out[1]] == 0, x]

Out[2]= x == −1, 999||x == 2, 001||x == 9, 02765

In[3]:= Roots[Numerator[Out[1]] == 0, x]

Out[3]= x == −1, 0001||x == 9, 02765

Пример 4

Как в примере 2, найдем аппроксимацию Паде π4,5(f2) типа (4, 5) функции f2(x) =(x−3,001)(x+1,9999)(x2+1)(x+4,0001) в точке x = 0.

In[4]:= PadeApproximant[

(x−3,001)(x+1,9999)((x2+1)(x+4,0001) , {x, 0, {4, 5}}]

Out[4]= −1,50039−25,4213x−6,51318x

2+3,7662x

3+0,42731x4

1,00000000000000+17,0263x+6,90326x2+17,4536x3+5,90326x4+0,42731x5

In[5]:= Roots[Denominator[Out[4]] == 0, x]

Out[5]:= x == −9, 75487||x == −4, 0001||x == −0, 0599743||x == −1, 0I||x == +1, 0I

In[6]:= Roots[Numerator[Out[4]] == 0, x]

Out[6]= x == −9, 75487||x == −1, 9999||x == −0, 0599743||x == 3, 001

Итак, числитель и знаменатель найденной по алгоритму, реализованному в Maple и Mathematica аппроксимации Паде, могут оказаться не взаимно простыми. 

Применение Метода градиентного спуска на С++, С#, Fotran, Delphi, Mathlab

Применение Метода градиентного спуска на С++, С#, Fotran, Delphi, Mathlab

Рассмотрим функцию f, считая для определенности, что она зависит от трех переменных x,y,z. Вычислим ее частные производные дf/дх, дf/ду, дf/дz и образуем с их помощью вектор, который называют градиентом функции:

grad f(x, у, z) = дf (х, у,z) /дх*i+дf( x, у, z)/ду*j+дf(x, y,z)/дг*k.

Здесь i, j, k - единичные векторы, параллельные координатным осям. Частные производные характеризуют изменение функции f по каждой независимой переменной в отдельности. Образованный с их помощью вектор градиента дает общее представление о поведении функции в окрестности точки (х, у,z). Направление этого вектора является направлением наиболее быстрого возрастания функции в данной точке. Противоположное ему направление, которое часто называют антиградиентным, представляет собой направление наиболее быстрого убывания функции. Модуль градиента 
                          ______________________                                  __ 
|grad (х, у,z)| =Ц (дf/дх (х, у,z))² +(дf/ду( x, у, z))²+(дf/дг(x, y,z))².

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

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

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

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

дf ~  f(x₁, ...,x_i+Dx_i, ..., x_n) - f(x₁, ..., x_i, ..., x_n)

 Отметим, что при таких расчетах Dx_i ,нельзя брать слишком малым, а значения функции нужно вычислять с достаточно высокой степенью точности, иначе при вычислении разности

 Df(x₁, ...,x_i+Dx_i, ..., x_n) - f(x₁, ..., x_i, ..., x_n)

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

Сравнение рис. 1 и 3 показывает, насколько более эффективным является метод градиентного спуска.

Метод покоординатного спуска С++,С#,Fotran,Delphi

Метод покоординатного спуска С++,С#,Fotran,Delphi

Пусть нужно найти наименьшее значение целевой функции u=f(M)=f(x1, x2, . . . ,xn). Здесь через М обозначена точка n-мерного пространства с координатами x1, x2, . . . ,xn: M=(x1, x2, . . . ,xn).Выберем какую-нибудь начальную точку М0=(x10, x20, . . . ,xn0) и рассмотрим функцию f при фиксированных значениях всех переменных, кроме первой: f(x1, x20,x30, . . . ,xn0 ). Тогда она превратится в функцию одной переменной  x1 . Изменяя эту переменную, будем двигаться от начальной точки x1=x10  в сторону убывания функции, пока не дойдем до ее минимума при x1=x11, после которого она начинает возрастать. Точку с координатами ( x11, x20,x30, . . . ,xn0) обозначим через М1, при этом f(M0) >= f(M1).

Фиксируем теперь переменные: x1=x11, x3= x30, . . . ,xn=xn0 и рассмотрим функцию f как функцию одной переменной  x2: f(x11, x22, x30 . . . ,xn0). Изменяя  x2 , будем опять двигаться от начального значения x2=x20   в сторону убывания функции, пока не дойдем до минимума при x2=x21 .Точку с координатами {x11, x21, x30 . . . xn0} обозначим через М2, при этом   f(M1) >=f(M2). 

Проведем такую же минимизацию целевой функции по переменным   x3, x4,  . . . ,xn. Дойдя до переменной xn, снова вернемся к x1 и продолжим процесс. Эта процедура вполне оправдывает название метода. С ее помощью мы построим последовательность точек М0,М1,М2, . . . , которой соответствует монотонная последовательность значений функцииf(M0) >= f (M1)>= f(M2) >=  ...    Обрывая ее на некотором шаге k можно приближенно принять  значение функции f(Mk) за ее наименьшее значение в рассматриваемой области.

Рис. 2. Поиск наименьшего значе- ния функции методом покоординат- ного спуска.

Отметим, что данный метод сводит задачу поиска наименьшего значения функции нескольких переменных к многократному решению одномерных задач оптимизации. Если целевая функция   f(x1, x2, ... ,xn)  задана явной формулой и является дифференцируемой, то мы можем вычисллить ее частные производные и использовать их для

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

На рис. 2 изображены линии уровня некоторой функции двух переменных             u= f (х, у). Вдоль этих линий функция сохраняет постоянные значения, равные 1, 3, 5, 7, 9. Показана траектория поиска ее наименьшего значения, которое достигается в точке О, с помощью метода покоординатного спуска. При этом нужно ясно понимать, что рисунок служит только для иллюстрации метода. Когда мы приступаем к решению реальной задачи оптимизации, такого рисунка, содержащего в себе готовый ответ, у нас, конечно, нет.

 Пусть требуется решить задачу (2):

f(x) -->min,  х ОRn.                                      (2)

В двумерном пространтсве R2. Решение задачи (2) методом покоординатного спуска, иначе называемого методом Гаусса - Зейделя, производят по следующей общей схеме. 
Выбирают произвольно начальную точку х(0)   из   области определения функции f(х). Приближения х(k) определяются соотношениями (3): 
x(k+1)=x(k)+t(k)S(k)      (k=0,1,2, ...), 
где вектор направления спуска s(k)- это единичный вектор, совпадающий с каким-либо координатным направлением (например, если  S(k) параллелен х1, то S(k)= {1,0,0,...,0}, если он параллелен  x2, то S(k)={0, 1, 0, . . . ,0} и т.д.) ; величина t(k)   является решением задачи одномерной минимизации: 
f(x(k)+ts(k)) -->  min,   t ОR1,    (k=0,1,2, ...), 
и может определяться, в частности, методом сканирования. 
Детальная реализация общей схемы в двумерном случае R2 дает траекторий приближения к точке х* методом покоординатного спуска, состоящую из звеньев ломаной, соединяющих точки х(k), х~(k),x(k+1)   (k=0, 1,  2,)       (рис.2).    При k=0, исходя  из начальной точки х(0)=(x1(0),x2(0)), находят точку х~(0)= (x1~(0),x2(0)), минимума функции одной переменной f(x1,x2(0)); при этом f(x~(0))<=f(x(0)).Затем находят точку минимума x(1)  функции f (x1~(0),x2) по второй  координате. Далее делают следующий шаг вычислений при k=1. Полагают, что исходной точкой расчета является  х(1).   Фиксируя вторую  координату   точки  х(1),   находят   точку   минимума х~(1)= (x1~(1),x2(1)), функции  f(x1,x2(1)) одной переменной x(1); при  этом f(x~(1))<=f(x(1))<=f(x(0)). Точку х(2)   получают, минимизируя целевую функцию f(x1~(1),x2), вновь по коорданате х2, фиксируя координату x1~(1) ,точки x(1) , и т.д. 
Условием прекращения вычислительной процедуры при достижении заданной точности e может служить неравенство 
||x(k+1) - x(k) ||<e           (4) 
  
Блок-схема поиска минимума функции двух переменных методом покоординатного спуска.