Промерзание грунта

[RAEman]

svv в сообщении #1020107 писал(а):

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

А точно! Руководитель что-то говорил насчет метода сглаживания $\delta-\text{функции}$. Не знаю это она или не она

[RAEman]

svv
Вот листинг программы якобы нормально работает

Используется синтаксис C

#include <math.h>
#include <stdio.h>
#include <stdlib.h> 
#define Pi 3.14159265358979
#define q_w 3.34*100000
#define ro_w 1000
#define n 12
#define N 100
#define lambda2 0.58
#define lambda1 0.5
#define delta 0.01
#define h 0.1
#define c_s 700
#define ro_s 4000
#define c_w 4190
#define ro_w 1000
#define c_ice 2000
#define ro_ice 900
#define lambda_s 900
#define lambda_w 900
#define lambda_ice 900
#define m 0.4
#define Tfaza 0
#define Tnach -3
double c_ro(double S_w, double M)
{
        double temp;
        temp = (1-M)*c_s*ro_s+c_w*ro_w*S_w*M+c_ice*ro_ice*(1-S_w)*M;
        return temp;
}
double _lambda(double S_w , double M)
{
        double temp;
        temp = (1 - M)*lambda_s + lambda_w*S_w*M + lambda_ice*(1 - S_w)*M;
        return temp;
}
double Tgran(double a,double t)
{
        return sin(Pi*a*t);
}
int main()
{
        FILE *f = fopen("1.txt", "w");
        int k,i,j;
        double tau = 1, t = 0, maxy = 1, alpha[N+1], beta[N+1], lambda, cro, y[N+1], yj[N+1],a , b, c, d, e;
        for (j = 0; j <= N; j++)
        { 
                y[j] = Tnach;
                yj[j] = y[j]; 
        }
        for (j = 1; j < N; j++)
        {
                tau = tau*1.1;
                if (tau>1e6) tau = 1e6;
                t = t + tau;
                k = 0;
                while ((maxy > 1e-4) && (k < 20))
                {
                        k = k + 1;
                        alpha[0] = 0; beta[0] = Tgran(1,t);
                        for (i = 1; i < n; i++)
                        {
                                if (y[i] < Tfaza) lambda = _lambda(0.0, m);
                                else lambda = _lambda(1.0, m);
                                if (y[i] < (Tfaza - delta)) cro = c_ro(0.0, m);
                                else if (y[i] > (Tfaza + delta)) cro = c_ro(1.0, m);
                                else cro = c_ro(0.0, m) + m*q_w*ro_w / (2 * delta);
                                a = lambda / (h * h);
                                b = a;
                                c = cro / tau;
                                d = c*yj[i];
                                e = a + b + c;
                                alpha[i] = b / (e - a*alpha[i - 1]);
                                beta[i] = (d + a*beta[i - 1]) / (e - a*alpha[i - 1]);
                        }
//                      y[n] = Tnach; //1 rod
                        y[n] = beta[n - 1] / (1 - alpha[n - 1]); //2 rod
                        maxy = y[n];
                        for (i = n - 1; i >= 0; i--)
                        {
                                a = y[i];
                                y[i] = y[i + 1] * alpha[i] + beta[i];
                                if (maxy < fabs(y[i] - a)) maxy = fabs(y[i] - a);
                        }
//                      printf("%d::",k);
                        for (int i = 1; i <= n; i++) fprintf(f,"%.5f |",y[i]);
                        fprintf(f,"\n");
                        yj[j] = y[j];
                }
        }
}
 

Однако ж, изменение N никак не влияет на изменение количества результатов. Можете подсказать в чем дело?

[RAEman]

svv
Оказывается вычисление yj[j] = y[j]; находилось внутри while{...}, Хотя по хорошему должно находится в конце цикла по дискретному времени for (j = 1; j < N; j++){...}.(И я опять задаюсь вопросом в чем измеряется это дискретное время. Видимо я сам задаю это время и мне надо его вычислить с помощью t которое растет до $10^6$).

[svv]

RAEman в сообщении #1020304 писал(а):

Оказывается вычисление yj[j] = y[j]; находилось внутри while{...}, Хотя по хорошему должно находится в конце цикла по дискретному времени for (j = 1; j < N; j++){...}.

Точно! Мне сразу не понравилось то место, где находился этот оператор.

RAEman в сообщении #1020304 писал(а):

И я опять задаюсь вопросом в чем измеряется это дискретное время. Видимо я сам задаю это время и мне надо его вычислить с помощью t которое растет до $10^6$.

Было бы понятно, если бы именно время t ограничивалось пределом 1e6. Но, посмотрите, там же ограничивается tau:
if (tau>1e6) tau = 1e6;
А tau — это не само время, а шаг по времени, прибавка. Это подтверждается операторами
t = t + tau; // получили следующий момент времени прибавкой шага
beta[0] = Tgran(1,t); // и использовали время для вычисления граничного значения
А для шага вообще непонятно, 1) зачем его увеличивать в геометрической прогрессии и 2) зачем его доводить до огромного значения 1e6. (Кстати, в программе все величины вычисляются в системе СИ, так что 1e6 — это миллион секунд, то есть одиннадцать с половиной суток.) Шаг по времени во многих методах постоянный, в других методах — адаптивный, он уменьшается там, где функция совершает частые колебания, и увеличивается на более плавных и спокойных участках. Но чтобы каждый раз тупо, независимо ни от чего, умножать шаг tau на 1.1 — этого я не понимаю.

Можно было бы подумать, что на самом деле ограничиваться должно время, а в программе ошибка. Но дело в том, что вот так делать
if (t>1e6) t = 1e6;
тоже нет никакого смысла. Если время достигло предельного значения, надо просто выходить из цикла, а не удерживать его искусственно на этом предельном значении.

[RAEman]

svv
Извините пожалуйста, я наверно Вам уже надоел с просьбами, но не могли бы вы мне разъяснить еще одну матмодель

$ \begin{cases} \frac{\partial (Sc)}{\partial t} = D \frac{\partial}{\partial x}(S \frac{\partial c}{\partial x})( \text{ну или можно упростить и приравнять к нулю} \frac{\partial (Sc)}{\partial t} = 0)\\ c\rho \frac{\partial T}{\partial t} = \frac{\partial}{\partial x}(\lambda \frac{\partial T}{\partial x}) - m \rho_w q_w \frac{\partial S}{\partial t}\\ T = -\alpha \cdot c\\ \end{cases}\\ T^{\ast}=-c_H \cdot \alpha\\ \alpha = 66.7(NaCl)\\ c - \text{концентрация соли} $
Я даже не знаю как тут сетку построить изза $ m \rho_w q_w \frac{\partial S}{\partial t}$

[RAEman]

Я так понимаю, что лучше будет если $\frac{\partial(Sc)}{\partial t}$ приравнять к нулю.
Также мне кажется, что сетка останется та же самая, что и в прошлой модели, но изменятся кое какие условия ибо появилась же концентрация соли и прочее
Вобщем ничего не понятно

[svv]

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

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

[RAEman]

svv
Спасибо, я уже разобрался: на самом деле это практически таже самая модель только с добавлением соляного раствора. Подставляя уравнения друг в друга мы получим что

$\frac{\partial S}{\partial t} = -\frac{S}{T} \frac{\partial S}{\partial t} $
где
$\\S = \frac{S_H T_H \alpha}{T}\\T_H=-c_H \alpha$

Cетка практически остается той же самой, за исключением того, что к $c\rho$ добавляется $m \rho_w q_w S_{\text{нач}} c_H \alpha$.
Т.о. у нас сетка получается такой

$y_i(\frac{c\rho + m \rho_w q_w S_{\text{нач}} c_H \alpha}{\tau}+\frac{2\lambda}{h^2}) = (\frac{c\rho + m \rho_w q_w S_{\text{нач}} c_H \alpha}{\tau} \tilde{y_i} + \lambda \frac{y_{i+1}+y_{i-1}}{h^2})$

Также у нас меняется немного среда где мы рассматриваем "промерзание":
В прошлой матмодели мы рассматривали лед -> вода-лед -> вода, а сейчас только вода-лед -> вода. Значит меняются эти условия

$\\ \begin{cases} c_1 \rho_1,&\text{если $T<T_f-\triangle$;}\\ c_2 \rho_2,&\text{если $T>T_f+\triangle$;}\\ c_1 \rho_1 + \frac{m \rho_w q_w}{2\triangle},&\text{если $T \in \left\lbrace T_f-\triangle;T_f+\triangle\right\rbrace$.} \end{cases}\\ \begin{cases} \lambda_1,&\text{если $T<T_f$;}\\ \lambda_2,&\text{если $T \geqslant T_f$;} \end{cases} $

на эти

$\\ \begin{cases} c_1 \rho_1,&\text{если $T \leqslant T_f$;}\\ c_1 \rho_2,&\text{если $T > T_f$;} \end{cases}\\ \begin{cases} \lambda_1,&\text{если $T \leqslant T_f$;}\\ \lambda_2,&\text{если $T > T_f$;} \end{cases}$

вот и все
Осталось только это все растянуть как-нибудь на побольше, чтоб не стыдно было в диплом на целую главу вставить

[svv]

RAEman в сообщении #1022661 писал(а):

$\frac{\partial S}{\partial t} = -\frac{S}{T} \frac{\partial S}{\partial t} $

А почему слева и справа одна и та же производная? Вроде как отсюда следует, что либо она равна нулю, либо $S=-T$. В любом случае, это нельзя назвать нормальным уравнением. Или я чего-то не понял?

[RAEman]

svv в сообщении #1022775 писал(а):

А почему слева и справа одна и та же производная? Вроде как отсюда следует, что либо она равна нулю, либо $S=-T$. В любом случае, это нельзя назвать нормальным уравнением. Или я чего-то не понял?

Я напутал, там не $\frac{\partial S}{\partial t} = -\frac{S}{T} \frac{\partial S}{\partial t} $, а $\frac{\partial S}{\partial t} = -\frac{S}{T} \frac{\partial T}{\partial t} $

[svv]

Но это уравнение легко решается. (А зачем Вам лишнее дифференциальное уравнение?)
$\frac{\partial S}{\partial t} = -\frac{S}{T} \frac{\partial T}{\partial t}$
Домножаем обе части на $T$, переносим всё в левую часть.
$T\frac{\partial S}{\partial t} + S\frac{\partial T}{\partial t}=0$
Применяем формулу производной произведения.
$\frac{\partial(TS)}{\partial t}=0$
Решение: $TS$ — любая функция (координаты $x$), не зависящая от $t$.
Больше из этого уравнения ничего нельзя получить.

[RAEman]

svv
Дело в том, что мы не решаем это уравнение, а подставляем в это уравнение

$c\rho \frac{\partial T}{\partial t} = \frac{\partial}{\partial x}(\lambda \frac{\partial T}{\partial x}) - m \rho_w q_w \bf{\frac{\partial S}{\partial t}}\\$

Вместо

$\frac{\partial S}{\partial t}$

ставим

$\frac{\partial S}{\partial t} = -\frac{S}{T} \frac{\partial T}{\partial t}$

перебрасываем с отрицательным знаком в левую часть, получаем

$(c\rho + m \rho_w q_w \frac{S}{T}) \frac{\partial T}{\partial t} = \frac{\partial}{\partial x}(\lambda \frac{\partial T}{\partial x}) \\$

и где у нас уже есть

$\\S = \frac{S_H T_H \alpha}{T}\\T_H=-c_H \alpha$

и с помощью этого мы можем уже нормально вывести коэффициенты для прогонки

[svv]

А, понятно. Но так, для себя, можете иметь в виду: произведение $TS$ не зависит от $t$ (хотя от $x$ может зависеть). Возможно, при проверке пригодится.

[RAEman]

svv
Большое спасибо за помощь.