Двумерное уравнение теплопроводности

Форум пользователей пакета Mathematica

Модератор: Admin

daniel33
Сообщения: 1
Зарегистрирован: Пн май 15, 2017 10:59 pm

Двумерное уравнение теплопроводности

Сообщение daniel33 » Пн май 15, 2017 11:10 pm

Добрый день! Мне нужно написать солвер для следующего уравнения теплопроводности:
$$\partial_t u = (u^2\partial_x u)_x + (u^2\partial_y u)_y$$
С такими начальными условиями, что функция u везде равна нулю, кроме диска в середине того квадрата, в котором я решаю уравнение. Внутри этого диска значение температуры, например, 5. Или 10. Я использую явный метод:
$$\frac{u_{i,j}^{k+1} - u_{i,j}^k}{\tau} = \frac{F_{i +1/2,j}^k - F_{i-1/2,j}^k} {h} + \frac{F_{i,j+1/2}^k - F_{i,j-1}^k}{h}$$, где
$$F_{i+1/2,j}^k = (u_{i+1/2,j}^k)^2(\frac{u_{i+1}^k - u_{i-1}^k}{h})$$
и $$u_{i+1/2}^k = \frac{u_{i+1}^k + u_{i}^k}{2}$$

На самом деле, я не очень хорошо знаю, как задать диск и аппроксимировать гран условия на нем, чтобы все было хорошо.
Я попытался задать большую температуру в какой-то одной точке в середине диска, но моя схема почем-то не видит этой самой точки.
Код такой:

Код: Выделить всё

nn = 10;
mm = 10000;
h = 1./nn;
tau = 1./mm;
v0[x_, y_] := 1.;
T[x_, y_, t_] :=
  Piecewise[ {{5., And[x ==  (nn/2), y  ==  (nn/2)]}, {0,
     And[(x !=  nn/2), (y !=  nn/2)]}}, 0];



Table[{x[i] = i*h, y[j] = j*h}, {i, 0, nn}, {j, 0, nn}];
Table[t[i] = i*t, {i, 0, mm}];
Table[u[i, j, 0] = v0[x[i], y[j]], {i, 0, nn}, {j, 0, nn}];
Table[u[0, j, k] = T[0, y[j], t[k]], {j, 0, nn - 1}, {k, 0, mm}];
Table[u[nn, j, k] = T[nn, y[j], t[k]], {j, 0, nn - 1}, {k, 0, mm}];
Table[u[i, 0, k] = T[x[i], 0, t[k]], {i, 0, nn - 1}, {k, 0, mm}];
Table[u[i, nn, k] = T[x[i], nn, t[k]], {i, 0, nn - 1}, {k, 0, mm}];

Do[ u[i, j, k + 1] =
   u[i, j, k] +
    tau/h*(      (((u[i + 1, j, k])^2 + (u[i, j, k])^2)/
          2)     *    ((u[i + 1, j, k] - u[i, j, k])/
          h) - (((u[i, j, k])^2 + (u[i - 1, j, k])^2)/
          2)  *((u[i, j, k] - u[i - 1, j, k])/
          h) +   (   ((u[i, j + 1, k])^2 + (u[i, j, k])^2)/
          2  )   *((u[i, j + 1, k] - u[i, j, k])/
          h) -   (((u[i, j, k])^2 + (u[i, j - 1, k]))^2/
          2)*((u[i, j, k] - u[i, j - 1, k])/h)), {k, 0, mm}, {i, 1,
   nn - 1}, {j, 1, nn - 1}];



g[k_] := Table[u[i, j, k], {i, 0, nn}, {j, 0, nn}];
p = Table[
  ListPlot3D[g[k], PlotRange -> {{0, nn}, {0, nn}, {0, 6}}], {k, 0,
   200, 50}]


Ну я смотрю на начало процесса, потому что там уже все понятно.
Пытался просто после гран условий написать строчку u[nn/2,nn/2,0] = 5;
При этом мне выдает ошибку оверфлоу. Хотя, если поставить, например, не 5, а 2, то считать будет.

Не подскажите, как сделать так, чтобы это все работало. А еще лучше, как это сделать с диском
Заранее благодарен!

Вернуться в «Mathematica»

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 1 гость