Помогите найти ошибку.

[brddrdtt]

С целью повышения квалификации, а так же по причине жуткого интереса, я решил разобраться в методе конечных элементов и самому написать программу. И написал. Но, где-то засела методическая ошибка и программа не работает должным образом. При этом она показывает качественно верный результат! Но не количественно. А именно: при увеличении количества элементов решение не только не сходится, но и изменяется согласно квадратному закону. Сама же картина деформирования отличается от верной только масштабом.
http://imageshack.us/photo/my-images/819/lj7d.jpg/
По оси х - размер элемента, по оси у - перемещение (е-6 ) (в данном случае правой нижней точки при закреплении левой и верхней стороны прямоугольника). Линия тренда говорит, что результат квадратично зависит от размера элемента. Abaqus показывает, что истинное перемещение для параметров по - умолчанию должно быть 12,05 е-6.

Программа рассчитывает перемещения узлов пластины в её плоскости при действии распределённых сил приложенных к этой пластине, действующих так же в её плоскости.
Структура программы: 1. при запуске предлагается ввести ширину, длину, толщину, модуль упругости, коэф. Пуассона и силы. А так же примерный размер треугольного элемента. 2. при нажатии elements size программа найдёт наиболее подходящий размер элемента. 3. create mesh - разбивка на сетку. 4. plot mesh - вывод сетки на экран, номеров узлов и номеров элементов. 5. element stiffness matrix - не участвующая в расчёте процедура - выводит на экран матрицу жесткости второго элемента. 6. global stiffness matrix - рассчитывает глобальную матрицу жёсткости. 7. up - left - right - down - запрещает соответствующим сторонам пластины перемещения. 8. loads vector формирует вектор нагрузок. 9. solution equations - решает систему уравнений методом Гаусса. 10. visualisation results - рисует деформированную пластину.

Я сначала грешил на процедуру определения матрицы жёсткости элемента, но сверка с другими источниками литературы (Р. Галлагер "Метод конечных элементов основы" стр. 134) показала верность процедуры. (Правда матрица жёсткости немного отличается от книги - цифры те же, но расположены в другом порядке в матрице - я это связываю с тем, что у Галлагера в другом порядке учтены перемещения узлов: U1U2U3V1V2V3 - у меня U1V1U2V2U3V3).

программа http://files.mail.ru/D0E5518FE4224FBFBB3C266A6E1DC903

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

[Theoristos]

Вы жесткость элемента масштабируете обратно пропорционально размерам?

[brddrdtt]

Матрица жёсткости элемента представляется в следующем виде.

$$[k] = tA[B]^T[D][B]$$

t -толщина пластины, А площадь элемента, [B] матрица градиентов, [D] матрица упругих констант.[k] - матрица жёсткости элемента.

$$[D] = \frac{E}{1-\mu^2}\left( \begin{array}{ccc} 1 & \mu & 0\\ \mu & 1 & 0\\ 0 & 0 & \frac{1-\mu}2 \end{array} \right) $$

E модуль упругости первого рода, $\mu$ - коэффициент Пуассона.

$$[B] = \frac1{2A}\left( \begin{array}{cccccc} y_j_k & 0 & y_k_i & 0 & y_i_j & 0\\ 0 & x_k_j & 0 & x_{ik} & 0 & x_{ji}\\ x_{kj} & y_{jk} & x_{ik} & y_{ki} & x_{ji} & y_{ij}\end{array} \right)$$

где $y_{ik} = y_i - y_k$ и т.д.
где $y_i $ координата У i-го узла элемента, и т.д.

Согласно Балабуху "Строительная механика ракет" стр. 98.
Я сверил с (Р. Галлагер "Метод конечных элементов основы" стр. 134) в маткаде вроде всё сошлось. Но у Галлагера мжэ немного по-другому сформирована.

[Theoristos]

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

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

[brddrdtt]

По - видимому, программа требует net.framework 4.5.

Что такое нормализация входных сил\ конечного результата?
Я силы представляю в виде глобального вектора столбца сил, который получается суммированием векторов столбцов сил элементов. Кстати, отсюда выходит, что строка глобального столбца, соответствующая некоторому узлу равна силе, приложенной к пластине, умноженной, на количество элементов, содержащих в себе этот узел. После этого суммирования строки глобального вектора сил, соответствующие узлам в которых запретили перемещения признаются неизвестными.
Получается система уравнений, содержащая n уравнений и n неизвестных. В каждом узле неизвестно либо перемещение, либо сила. Система решается методом Гаусса. Полученные перемещения без каких - либо "нормализаций" считаются истинными.
Это верно?

На выходных займусь всеми перепроверками в маткаде.

[Theoristos]

Я бы попробовал посмотреть вид матриц и решения при самом крупном разбиении.

[brddrdtt]

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

[brddrdtt]

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