Основы компьютерного моделирования физических и технических систем. Математическая модель системы

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ

УЧРЕЖДЕНИЕ ОБРАЗОВАНИЯ

ГОМЕЛЬСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ ИМЕНИ П. О. СУХОГО

Факультет автоматизированных и информационных систем

Кафедра «Информационные технологии»

дисциплина «Компьютерные системы конечноэлементных расчётов»

ОТЧЕТ ПО ЛАБОРАТОРНОЙ РАБОТЕ № 1

«Основы компьютерного моделирования физических и технических систем»

Выполнила: студентка гр. ИТ-31

Принял: преподаватель

Гомель 2012

Лабораторная работа №1

Основы компьютерного моделирования физических и технических систем

Задание.

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

- обеспечить ввод исходных данных с помощью GUI;

- отобразить в виде графиков результаты решения;

- все результаты решения сохранять в тестовые файлы

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

Необходимо:

- провести верификацию, разработанного программного обеспечения на примере решения задачи для однородного стержня;

- сравнить численные методы по скорости нахождения решения при заданной точности;

- сравнить численные методы по значениям найденных решений;

- найти максимальное отклонение найденных решений;

- построить усреднённое решение на основе найденных;

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

Текст программы МКР и МКЭ:

public partial class Form1 : Form

{

public float r = 0.1f, lam_med = 0.400f, lam_olovo = 0.067f, q;

public float c_med = 0.385f, c_olovo = 0.218f, ro_med = 8.920f, ro_olovo = 7.300f;

public float Tw = 100, h = 8.7f, Sa, L1, olovo_med, dt = 0.12f;

public int tr, xr, length, length_olovo, length_med;       

public Graphics g;

public float[] XX,X, B, f, T, PV;

public float[,] A,AA, C, K, P, BB;

public Form1()

{

InitializeComponent();

g = panel1.CreateGraphics();

g.Clear(Color.White);

listBox1.SelectedIndex = 0;

Sa = r * r * (float)Math.PI;            

}

public void Initial()

{

for (int i = 0; i < (xr + 1); i++)

{

A[i, i] = 1.0f;

B[i] = float.Parse(comboBox1.SelectedItem.ToString());

AA[0,i] = B[i];

}           

}

public void GLeft()

{

for (int i = 0; i < (tr + 1) * (xr + 1); i += xr + 1)

{               

A[i, i] = 1.0f;

B[i] = float.Parse(listBox1.SelectedItem.ToString());

}           

}

public void GRight()

{

for (int i = xr; i < (tr + 1) * (xr + 1); i += xr + 1)

{

A[i, i] = 1.0f;

B[i] = float.Parse(domainUpDown1.Text.ToString());

}           

}

private void button1_Click(object sender, EventArgs e)

{ 

if (radioButton1.Checked)

{

Initial();

GLeft();

GRight();

Centr();

StreamWriter ff = new StreamWriter("A.txt");

StreamWriter f2 = new StreamWriter("B.txt");

for (int j = 0; j < A.GetLength(0); j++)

{

for (int u = 0; u < A.GetLength(1); u++)

{

ff.Write(A[j,u]+"  ");

}

f2.WriteLine(B[j]);

ff.WriteLine();

}

f2.Close();

ff.Close();

ff = new StreamWriter("mkr.txt");

X= BMG(A,B);

for (int j = 0; j < X.Length; j++)

{

ff.WriteLine(X[j]);

if (j>0 && (j+1)%(xr+1) == 0)

ff.WriteLine();

}

ff.Close();

DrawMethod(X, 1);                   

}

if (radioButton2.Checked)

{

T = new float[BB.GetLength(0)];

Init();                   

Matrix();

AA[0,0] = Int32.Parse(listBox1.SelectedItem.ToString());

AA[0,xr] = Int32.Parse(domainUpDown1.Text);

StreamWriter fm = new StreamWriter("mke.txt");

for (int j = 0; j < xr + 1; j++)

{

fm.WriteLine(AA[0, j]);

XX[j] = AA[0, j];

}

fm.WriteLine();

for (int i = 1; i < tr + 1; i++)

{

AA[i, 0] = Int32.Parse(listBox1.SelectedItem.ToString());

AA[i, xr] = Int32.Parse(domainUpDown1.Text);

T = BMG(BB, VectorPTF(i));

T[0] = Int32.Parse(listBox1.SelectedItem.ToString());

T[xr] = Int32.Parse(domainUpDown1.Text);

for (int j = 0; j < xr + 1; j++)

{

XX[i * (xr + 1) + j] = T[j];

AA[i, j] = T[j];

fm.WriteLine(T[j]);

}

fm.WriteLine();

}

fm.Close();

DrawMethod(XX, 2);

}               

MessageBox.Show("Вычисления завершены");

}

public void Centr()

{

int num = 0;

float  lam1,lam2,lam;

for (int i = 1; i < tr + 1; i++)

{

for (int j = 1; j < xr; j++)

{

num = j + i * (xr + 1);

if (j*1.0f < olovo_med - 1)

{

lam1 = lam = lam2 = dt * lam_olovo / (c_olovo * ro_olovo * L1 * L1);

}

else

if (j*1.0f > olovo_med + 1)

{

lam1 = lam = lam2 = dt * lam_med / (c_med * ro_med * L1 * L1);

}

else

if (j*1.0f == olovo_med)

{

lam1 = dt * lam_olovo / (c_olovo * ro_olovo * L1 * L1);

lam2 = dt * lam_med / (c_med * ro_med * L1 * L1);

lam = (lam_med / (c_med * ro_med) + lam_olovo / (c_olovo * ro_olovo)) / (L1 * L1 * 2);

}

else

if (j*1.0f <= olovo_med)

{

lam1=lam = dt * lam_olovo / (c_olovo * ro_olovo * L1 * L1);

lam2 = dt * lam_med / (c_med * ro_med * L1 * L1);

}

else

{

lam2=lam = dt * lam_med / (c_med * ro_med * L1 * L1);

lam1 = dt * lam_olovo / (c_olovo * ro_olovo * L1 * L1);

}

A[num, num] = 1 - lam;

if (i == 1)

{

if (j == 1)

{

A[num, num + 1] = -lam2 / 2;

B[num] = lam2 / 2 * B[num - xr] - (lam - 1) * B[num - xr - 1] + lam1 / 2 *

B[num - xr - 2] + lam1 / 2 * B[num - 1];

}

else

{

A[num, num - 1] = -lam1 / 2;

if (j == xr - 1)

{

B[num] = lam1 / 2 * B[num - xr - 2] (lam - 1) * B[num - xr - 1] + lam2 / 2 * B[num - xr] +

lam2 / 2 * B[num + 1];

}

else

{

A[num, num + 1] = -lam2 / 2;

B[num] = lam2 / 2 * B[num - xr] - (lam - 1) * B[num - xr - 1] +

lam1 / 2 * B[num - xr - 2];

}

}

}

else

{

A[num, num - xr - 1] = lam - 1;

if (j == 1)

{

A[num, num + 1] = -lam2 / 2;

A[num, num - xr] = -lam2 / 2;

B[num] = lam1 / 2 * B[num - 1] + lam1 / 2 * B[num - xr - 2];

}

else

{

A[num, num - 1] = -lam1 / 2;

A[num, num - xr - 2] = -lam1 / 2;

if (j == xr - 1)

{

B[num] = lam2 / 2 * B[num + 1] + lam2 / 2 * B[num - xr];

}

else

{

A[num, num + 1] = -lam2 / 2;

A[num, num - xr] = -lam2 / 2;

B[num] = 0;

}

}

}

}

}

}

public void DrawMethod(float[] Y, int u)

{

float min = Y[0];

float max = Y[0];           

for (int i = 1; i < Y.Length; i++)

if (Y[i] < min)

min = Y[i];

else

if (Y[i] > max)

max = Y[i];

int sizex = panel1.Width;

int sizey = panel1.Height;

float k = sizex / (xr+1);

Pen pen;

if (u == 1)

pen = Pens.Red;

else

pen = Pens.Blue;                      

for (int j = 0; j < Y.Length - 1; j++)

{

g.DrawLine(pen, 2*j*sizex/(max-min), Y[j] * sizey / (max - min), 2*(j+1)*sizex/(max-min), Y[j + 1] * sizey / (max - min));

}

}

public float[] VectorPTF(int k)

{

float sum;

for (int i = 0; i < xr+1; i++)

{

sum = 0;

for (int j = 0; j < xr + 1; j++)

{

sum += P[i, j] * AA[k - 1, j];

}

PV[i] = sum - 2*f[i];

}

return PV;

}

public void MatrixBP()

{

BB[0, 0] = K[0, 0] + C[0, 0];

P[0, 0] = C[0, 0] - K[0, 0];

BB[0, 1] = K[0, 1] + C[0, 1];

P[0, 1] = C[0, 1] - K[0, 1];

for (int i = 1; i < xr; i++)

{

BB[i, i - 1] = K[i, i - 1] + C[i, i - 1];

P[i, i - 1] = C[i, i - 1] - K[i, i - 1];

BB[i, i] = K[i, i] + C[i, i];

P[i, i] = C[i, i] - K[i, i];

BB[i, i + 1] = K[i, i + 1] + C[i, i + 1];

P[i, i + 1] = C[i, i + 1] - K[i, i + 1];

}

BB[xr, xr] = K[xr, xr] + C[xr, xr];

P[xr, xr] = C[xr, xr] - K[xr, xr];

BB[xr, xr - 1] = K[xr, xr - 1] + C[xr, xr - 1];

P[xr, xr - 1] = C[xr, xr - 1] - K[xr, xr - 1];

}

public void VectorF()

{

float kf = -q*Sa/2;

float kf2 = h * Sa/2;

f[0] = kf;

f[f.Length - 1] = kf2;

}

public void MatrixK()

{

float koef;

koef = Sa * lam_olovo / L1;

K[0, 0] = koef;

K[0, 1] = -koef;

K[1, 0] = -koef;

K[1, 1] += koef;

for (int i = 1; i < K.GetLength(0) - 2; i++)

{

if (L1 * i == length_olovo)

koef = Sa * (lam_olovo + lam_med) / (L1 * 2);

else if (L1 * i > length_olovo)

koef = Sa * lam_med / L1;

K[i, i] += koef;

K[i, i + 1] = -koef;

K[i + 1, i] = -koef;

K[i + 1, i + 1] += koef;

}

K[xr, xr] = koef;

K[xr, xr - 1] = -koef;

K[xr - 1, xr] = -koef;

K[xr - 1, xr - 1] += koef;

}

public void MatrixC()

{

float kf = c_olovo * ro_olovo * L1 * Sa / (3*dt);

C[0, 0] = 2*kf;

C[0, 1] = kf;

C[1, 0] = kf;

C[1, 1] += 2 * kf;

for (int i = 1; i < C.GetLength(0) - 2; i++)

{

if (L1 * i == length_olovo)

kf = (c_olovo*ro_olovo+c_med*ro_med)*L1*Sa/(6*dt);

else if (L1 * i > length_olovo)

kf = c_med * ro_med * L1 * Sa / (3 * dt);

C[i, i] += 2 * kf;

C[i, i + 1] = kf;

C[i + 1, i] = kf;

C[i + 1, i + 1] += 2 * kf;

}

C[xr, xr] = 2*kf;

C[xr, xr - 1] = kf;

C[xr - 1, xr] = kf;

C[xr - 1, xr - 1] += 2 * kf;

}

public void Init()

{

for (int i = 0; i < xr + 1; i++)

{

for (int j = 0; j < xr + 1; j++)

{

if (i != 0)

A[i, j] = 0;

else

if(i<tr+1)

A[i, j] = (float)Int32.Parse(comboBox1.SelectedItem.ToString());

C[i, j] = 0;

K[i, j] = 0;

P[i, j] = 0;

BB[i, j] = 0;

}

f[i] = 0;

PV[i] = 0;

}

}

public float[] Vihs(float[,]m1,float[] m2)

{

float[] resultat;

PramHod(m1.GetLength(0),ref m1,ref m2);

resultat = ObratHod(ref m1, ref m2);

return resultat;

}

public void PramHod(int n,ref float[,]m1,ref float[] m2)

{

for (int kk = 0; kk < n; kk++)

{

float max = Math.Abs(m1[kk, kk]);

int nmax = kk;

for (int i = kk + 1; i < n; i++)

if (Math.Abs(m1[i, kk]) > max)

{

max = Math.Abs(m1[i, kk]);

nmax = i;

}

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

{

float c = m1[kk, j];

m1[kk, j] = m1[nmax, j];

m1[nmax, j] = c;

}

float v = m2[kk];

m2[kk] = m2[nmax];

m2[nmax] = v;

for (int i = kk + 1; i < n; i++)

{

float m = m1[i, kk] / m1[kk, kk];

for (int j = kk; j < n; j++)

m1[i, j] = m1[i, j] - m * m1[kk, j];

m2[i] =m2[i] - m * m2[kk];

}

}

for (int kk = 0; kk < n; kk++)           

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

m1[kk, i] = 0;           

}

public float[] ObratHod(ref float[,] a, ref float[] b)

{

float[] xx = new float[A.GetLength(0)];

int i, j; float sum, perem;

for (i = 0; i < a.GetLength(0); i++)

xx[i] = 0;

for (i = 0; i < a.GetLength(0); i++)

for (j = 0; j < i; j++)

a[i, j] = 0;

i = 0; j = 0;           

int n = a.GetLength(0);

int m = a.GetLength(1) + 1;

float[,] aa = new float[n, m];

for (i = 0; i < n; i++)

for (j = 0; j < n; j++)

aa[i, j] = a[i, j];

for (i = 0; i < n; i++)

{

aa[i, m - 1] = b[i];

xx[i] = 0;

}

int t = n;

for (i = n - 1; i >= 0; i--)

{

sum = 0;

for (int p = 0; p < m; )

{

if (aa[i, p] == 0)

p++;

else

{

j = m - 2;

while (j != p)

{

sum += aa[i, j] * xx[j];

j--;

}

perem = (aa[i, m - 1] - sum) / aa[i, p];

xx[t - 1] = perem;

t--;

break;

}

}

}

return xx;

}

private void button5_Click(object sender, EventArgs e)

{

Pogresh(X, XX);

}

private void button6_Click(object sender, EventArgs e)

{

tr = Int32.Parse(numericUpDown4.Value.ToString());

xr = Int32.Parse(numericUpDown3.Value.ToString());

length_med =  Int32.Parse(numericUpDown2.Value.ToString());

length_olovo =  Int32.Parse(numericUpDown1.Value.ToString());

length = 1;

L1 = length * 1.0f / xr;

q = 100;

olovo_med = length_olovo;

XX = new float[(xr + 1) * (tr + 1)];

A = new float[(xr + 1) * (tr + 1), (xr + 1) * (tr + 1)];

B = new float[(xr + 1) * (tr + 1)];

C = new float[(xr + 1), (xr + 1)];

AA = new float[(tr + 1), (xr + 1)];

P = new float[(xr + 1), (xr + 1)];

BB = new float[(xr + 1), (xr + 1)];

K = new float[(xr + 1), (xr + 1)];

PV = new float[(xr + 1)];

f = new float[(xr + 1)];

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

{

for (int j = 0; j < B.Length; j++)

A[i, j] = 0;

B[i] = 0;

}

button6.Enabled = false;

}     

}

Результат выполнения:

Выводы: во время выполнения лабораторной работы была построена