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

Вид материала

Содержание

Подобный материал:

1 2 3 4 5 6 7 8 9 10 ... 27

Отсюда следует, что в момент времени t=(j+1) температура ФМ составит:

Т_n_,_j₊₁ = Т_n_,_j{	S	-	A	+1} - Т_n-1,j	S	(8.8)
	hС_ФМ		С_ФМ		hС_ФМ

Заменим в уравнении (8.2) вторые производные конечными разностями (8.6):

	Т_i+1,j–2Т_i,j+Т_i-1,j	= C	Т_i,j+1– Т_i,j
	h²		

Пусть  =(/h²C). Можно показать, что условие устойчивости решения, то есть возможность получения точных решений за определенное количество шагов, примет вид:

	h²C	=	1	 
	2		2

Далее примем =0.2/. Тогда значение рабочей температуры в момент времени t=(j+1) в i –м сечении стержня примет вид:

T_i,j+1 = 0.2T_i+1,j + 0,6 T_i,j+ 0.2T_i-1,j

(8.9)

По формуле (8.8) и (8.9) можно рассчитать температуру ФМ в сечениях вывода последовательно в моменты времени t_j = 0, t_j₊₁ = , t_j₊₂ = 2, … . В общем случае, в произвольный k-й момент времени t_k= (t_k_-1 + ). Программа на языке Паскаль вычисления температурного поля в стержне методом конечных разностей (МКР), действующая по формуле (8.9), может иметь следующий вид:

CONST n=10; Delta=0.1; L=20; Type mass=array[1..n]of LongInt;

Lambda=400; A=0.1; Cfm=6; C=4000000; S=0.000001; h=0.001;

VAR Md,Mt:mass; i,t:integer; b:boolean; k1,k2:Real;

BEGIN

k1:=((S*lambda*Delta)/(h*Cfm))-(delta*A/Cfm)+1; k2:=S*lambda*Delta/(h*Cfm);

For i:=1 to n Do Mt[i]:=293; Mt[1]:=600; For i:=1 to n Do MD[i]:= Mt[i];

Repeat For i:=2 to n-1 Do

Md[i]:=round(0.2*(Mt[i+1]+Mt[i-1])+0.6*Mt[i]);

md[n]:=round(k1*mt[n] - k2*mt[n-1]); For i:=2 to n Do Mt[i]:=Md[i];

Until false; END.

По формуле (8.9) аналогичным образом может быть составлена программа, реализующая вычисление значений температуры в дискретные моменты времени в последнем n-м сечении стержня, к которому подсоединена фиктивная масса.

Ошибки ограничения уменьшаются при h  0 и   0, то есть решение разностных уравнений (8.8) и (8.9) асимптотически приближается или сходится к решению ДУ (8.1) и (8.2). Следовательно, с целью повышения точности решений следует уменьшать шаг дискретизации по длине и по времени, что, в свою очередь, ведет к возрастанию количества основных циклов в алгоритме вычислений (пропорционально m x n), то есть – к увеличению времени счета.

Другим важным свойством МКР является устойчивость решений. Устойчивость метода означает, что любые ошибки (округления или ограничения) не возрастают в ходе вычисления решения. Доказано, что указанный процесс сходится и устойчив, если (<0.5) или (<1/[2]). Тем самым накладываются ограничения на выбор шага  по времени при выбранном шаге h разбиения вывода сечениями.

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

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

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

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

Рассмотрим систему из 4-х линейных алгебраических уравнений вида:

а⁺¹₁₁х₁ + а⁺¹₁₂ х₂ + а⁺¹₁₃ х₃+ а⁺¹₁₄ х₄ = b⁺¹₁

а⁺¹₂₁ х₁+ а⁺¹₂₂ х₂ + а⁺¹₂₃ х₃+ а⁺¹₂₄ х₄ = b⁺¹₂

а⁺¹₃₁ х₁ + а⁺¹₃₂ х₂ + а⁺¹₃₃ х₃+ а⁺¹₃₄ х₄ = b⁺¹₃

а⁺¹₄₁ х₁+ а⁺¹₄₂ х₂ + а⁺¹₄₃ х₃+ а⁺¹₄₄ х₄ = b⁺¹₄

Выразив из первого уравнения переменную x₁, имеем:

x₁=	[	b⁺¹₁		-a⁺¹₁₂		-a⁺¹₁₃		-a⁺¹₁₄	]		[1 x₂ x₃ x₄ ] ^т
x₁=	[	a⁺¹₁₁		a⁺¹₁₁		a⁺¹₁₁		a⁺¹₁₁	]		[1 x₂ x₃ x₄ ] ^т

Подставляя полученное выражение для х₁ во 2-е, 3-е и 4-е уравнения и приводя подобные члены, приходим к системе:

а⁺¹₁₁х₁ + а⁺¹₁₂ х₂ + а⁺¹₁₃х₃+ а⁺¹₁₄ х₄ = b⁺¹₁

0 + (а₂₂–а₂₁[а₁₂/а₁₁])х₂+(а₂₃–а₂₁[а₁₃/а₁₁])х₃+(а₂₄–а₂₁[а₁₄/а₁₁])х₄ = b₂–b₁(а₂₁/а₁₁)

0 + (а₃₂–а₃₁[а₁₂/а₁₁])х₂+(а₃₃–а₃₁[а₁₃/а₁₁])х₃+(а₃₄–а₃₁[а₁₄/а₁₁])х₄ = b₃–b₁(а₃₁/а₁₁)

0 + (а₄₂–а₄₁[а₁₂/а₁₁])х₂+(а₄₃–а₃₁[а₁₄/а₁₁])х₃+(а₄₄–а₄₁[а₁₄/а₁₁])х₄ = b₄–b₁(а₄₁/а₁₁)

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

a_pq⁺⁽^k⁺¹⁾= a_pq⁺^k- a_pk⁺^k( a_kq⁺^k/ a_kk⁺^k) ; b_p⁺⁽^k⁺¹⁾= b_p⁺^k- b_k⁺^k( a_pk⁺^k/ a_kk⁺^k) (14.1)

Тогда последнюю систему можно переписать в виде:

а⁺¹₁₁х₁ + а⁺¹₁₂ х₂ + а⁺¹₁₃х₃+ а⁺¹₁₄ х₄ = b⁺¹₁

0 + а⁺²₂₂ х₂ + а⁺²₂₃х₃+ а⁺²₂₄ х₄ = b⁺²₂

0 + а⁺²₃₂ х₂ + а⁺²₃₃х₃+ а⁺²₃₄ х₄ = b⁺²₃

0 + а⁺²₄₂ х₂ + а⁺²₄₃х₃+ а⁺²₄₄ х₄ = b⁺²₄

Выразив из второго уравнения переменную x₂, имеем:

X₂=	[	b⁺²₂		-a⁺²₂₃		-a⁺²₂₄	]		[1 x₃ x₄ ] ^т
		a⁺²₂₂		a⁺²₂₂		a⁺²₂₂

Подставляя полученное выражение для х₂ в 3-е и 4-е уравнения и приводя подобные члены, приходим к системе:

а⁺¹₁₁х₁ + а⁺¹₁₂ х₂ + а⁺¹₁₃х₃+ а⁺¹₁₄ х₄ = b⁺¹₁

а⁺²₂₂ х₂ + а⁺²₂₃х₃+ а⁺²₂₄ х₄ = b⁺²₂

(а₃₃⁺²–а₃₂⁺² [а₂₃⁺²/а₂₂⁺²])х₃+ (а₃₄⁺²–а₃₂⁺² [а₂₄⁺²/а₂₂⁺²])х₄ = b⁺²₃–b₂⁺²(а₃₂⁺²/а₂₂⁺²)

(а₄₃⁺²–а₄₂⁺² [а₂₃⁺²/а₂₂⁺²])х₃+ (а₄₄⁺²–а₄₂⁺² [а₂₄⁺²/а₂₂⁺²])х₄ = b⁺²₄–b₂⁺²(а₄₂⁺²/а₂₂⁺²)

Коэффициент при неизвестной х3 в третьем уравнении логично было бы обозначить как а₃₃⁺³. Попробуем получить его формально, используя первую формулу в выражении (14.1). С этой целью обозначим p=3 (№ строки) , q=3 (№ столбца), k=2 (номер текущей итерации) и подставим эти индексы в (14.1):

a₃₃⁺⁽²⁺¹⁾= a₃₃⁺²- a₃₂⁺²( a₂₃⁺²/ a₂₂⁺²) ;

Получили очевидное совпадение результатов. Вычислим аналогично остальные коэффициенты при неизвестных в третьем и четвертом уравнениях:

a₃₄⁺⁽²⁺¹⁾= a₃₄⁺²- a₃₂⁺²( a₂₄⁺²/ a₂₂⁺²) =a₃₄⁺³

a₄₃⁺⁽²⁺¹⁾= a₄₃⁺²- a₄₂⁺²( a₂₃⁺²/ a₂₂⁺²) =a₄₃⁺³

a₄₄⁺⁽²⁺¹⁾= a₄₃⁺²- a₄₂⁺²( a₂₃⁺²/ a₂₂⁺²) =a₄₄⁺³

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

где: b₃⁺³= b⁺²₃–b₂⁺²(а₃₂⁺²/а₂₂⁺²) и b₄⁺³= b₄⁺²–b₂⁺²(а₄₂⁺²/а₂₂⁺²).

После третьей итерации система примет вид:

а₁₁⁺¹х₁+ а₁₂⁺¹х₂+ а₁₃⁺¹х₃+ а₁₄⁺¹х₄= b₁⁺¹

0+ а₂₂⁺²х₂+ а₂₃⁺²х₃+ а₂₄⁺²х₄= b₂⁺²

0+ 0+ а₃₃⁺³х₃+ а₃₄⁺³х₄= b₃⁺³

0+ 0+ 0+ а₄₄⁺⁴х₄= b₄⁺⁴

где: a₄₄⁺⁴= a₄₄⁺³–a₄₃⁺³(а₃₄⁺³/а₃₃⁺³) и b₄⁺⁴= b₄⁺³–b₃⁺³(а₄₃⁺³/а₃₃⁺³).

Решение полученной системы выполняем методом обратной прогонки. Из четвертого уравнения вычисляем неизвестную Х4:

х₄= b₄⁺⁴ /а₄₄⁺⁴

Из третьего уравнения вычисляем неизвестную Х3:

х₃= [b₃⁺³-(а₃₄⁺³х₄)]/а₃₃⁺³

Из второго уравнения вычисляем неизвестную Х2:

х₂+ = [b₂⁺² – (а₂₄⁺²х₄ + а₂₃⁺²х3)]/а₂₂⁺²

Наконец, из первого уравнения вычисляем неизвестную Х1:

х₁ = [b₁⁺¹ – (а₁₄⁺¹х₄+а₁₃⁺¹х₃+а₁₂⁺¹х₂)]/а₁₁⁺¹

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

Однако, непосредственно перед решением система должна быть преобразована, поскольку, как правило, некоторые компоненты неизвестного вектора Ф (в программе – это вектор Xr) узловых значений известны. Так, в большинстве задач теории поля некоторые граничные значения искомой величины заданы; во всех задачах теории упругости должны быть фиксированы некоторые перемещения с тем, чтобы исключить перемещение среды жесткого тела.

То есть, элементы матриц [K] и {F} системы: [K]{Ф}={F} необходимо преобразовать, чтобы ее решение после преобразование давало правильный результат. При этом желательно не изменять программу решения самой системы, поскольку это повлечет за собой трудности при программировании.

Пусть в исходной системе задано фиксированное значение р-й переменной (Х_р=Q). Преобразование системы проводим по шагам:

коэффициенты р-й строки, кроме диагонального коэффициента, равного а_pp⁺¹, приравниваем нулю;
свободный член в р-й строке заменяем произведением: (а_pp⁺¹Q);
уравнения, содержащие переменную Х_р, преобразуем, вычитая из обеих частей каждого из них произведение (а_qp⁺¹Q), где q – номер строки (qp).

Проиллюстрируем это на примере системы уравнений:

46,6T₁ – 21,7T₂ + 0 + 0 = 1000

- 21,7T₁ + 93,2T₂ – 21,7T₃+ 0 = 2000

0 - 21,7T₂ + 93,2T₃ – 21,7T₄ = 2000

0 + 0 – 21,7T₃ + 56,6T₄ = 1400

(14.2)

Blog

Содержание

S

Тi+1,j–2Тi,j+Тi-1,j

Тi,j+1 – Тi,j

h2C

1

Т_i+1,j–2Т_i,j+Т_i-1,j

Т_i,j+1– Т_i,j

h²C