2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


В этом разделе нельзя создавать новые темы.

Если Вы хотите задать новый вопрос, то не дописывайте его в существующую тему, а создайте новую в корневом разделе "Помогите решить/разобраться (М)".

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

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

Обязательно просмотрите тему Правила данного раздела, иначе Ваша тема может быть удалена или перемещена в Карантин, а Вы так и не узнаете, почему.



Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Промерзание грунта
Сообщение13.05.2015, 18:14 


13/05/15
17
г. Якутск
Здравствуйте ув.Форумчане.
Помогите разобраться в задаче Стефана:
$
c_1 \rho_1\frac{\partial T}{\partial t} = \lambda_1 \frac{\partial^2 T}{\partial x^2}  , 0<x<S(t)\\
c_2 \rho_2\frac{\partial T}{\partial t} = \lambda_2 \frac{\partial^2 T}{\partial x^2}  , S(t)<x<l\\
\lambda_1 \frac{\partial T}{\partial x} - \lambda_2 \frac{\partial T}{\partial x} = q_w \rho_w \frac{\partial S}{\partial t}\\
\\
l = 50 \text{ метров}\\
c_1 \rho_1 = 2 \cdot10^6, \lambda_1 = 2\\
c_2 \rho_2 = 2 \cdot 10^6, \lambda_2 = 0,5\\
\\
x = S(t)\\
T_- = T_{\text{Ф}}=T_+\\
T(x,0) = T^0(x), t=0, 0 \leq x \leq l\\
T(0,t) = T_1(t) = \sin(\pi \cdot a \cdot x)\\
\begin{cases}
T(l,t) = T_H &\\
&\text{ $t>0,x=l$;}\\
\frac{\partial T}{\partial x} = q &
\end{cases}
Есть набросок программы которую я не совсем понимаю
код: [ скачать ] [ спрятать ]
Используется синтаксис C#
#include "stdafx.h"
#include <math.h>
 
#define Pi 3.14159265358979
#define q_w 3.34*1000000
#define ro_w 1000
#define n 10
#define N 10
#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
 
int 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;
}
int _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;
}
int Tgran(double a,double x){
    return sin(Pi*a*x);
}
int main(){
    int k;
    double tau = 1, t = 1, maxy = 1, alpha[n+1], beta[n+1], lambda, cro, y[n+1], yj[n+1],a , b, c, d;
    for (int j = 0; j <= N; j++){
        y[j] = Tnach;
        yj[j] = y[j];
    }
    for (int 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,1)/*Синус пока от балды написал*/;
            for (int i = 1; i < n; i++ ){
                if (y[i] < Tfaza) lambda = _lambda(1.0,m);
                else lambda = _lambda(0.0, m);
                if (y[i] < (Tfaza - delta)) cro = c_ro(1.0,m);
                else if (y[i] < (Tfaza + delta)) cro = c_ro(0.0,m);
                else cro = c_ro(1.0,m) + m*q_w*ro_w / (2 * delta);
                a = lambda / (h * h);
                b = a;
                c = cro / tau;
                d = c*yj[i];
                c = a + b + c;
                alpha[i] = b / (c - a*alpha[i - 1]);
                beta[i] = (d + a*beta[i - 1]) / (c - a*alpha[i - 1]);
            }
            y[n] = Tnach;
            y[n] = beta[n - 1] / (1 - alpha[n - 1]);
            maxy = y[n];
            for (int i = n - 1; i >= 0; i--){
                a = y[i];
                y[i] = y[i + 1] * alpha[i] + beta[i];
                if (maxy < abs(y[i] - a)) maxy = abs(y[i] - a);
            }
            for (int i = 1; i <= n; i++){
                printf("%f", y[i]);
            }
            printf("\n");
            yj[j] = y[j];
           
        }
        getchar();
    }
}
 

Я понимаю, что из постановку задачи мы аппроксимируем, получаем сетку, (и потом вроде как сглаживание(? не знаю что это) проводим) и ее прогонкой решаем.
1)А что эта сетка означает и как она мне помогает построить график(скажем в той же проге) я не понимаю.
2)Зачем было введено maxy.
3)И как мне по времени разбить результат

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение15.05.2015, 01:03 
Заслуженный участник


23/07/08
10609
Crna Gora
1. Понятен ли Вам физический смысл всех уравнений, которые Вы записали?
2. Понятно ли, зачем в численных методах решения ДУЧП вообще вводятся сетки?
3. Смогли бы Вы написать расчётные формулы, по которым должна работать программа (не та, что Вы привели, а Ваша, ещё не написанная)?

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение15.05.2015, 07:00 


13/05/15
17
г. Якутск
svv
Цитата:
1. Понятен ли Вам физический смысл всех уравнений, которые Вы записали?

Нет, что я и хочу понять....
Цитата:
2. Понятно ли, зачем в численных методах решения ДУЧП вообще вводятся сетки?

Не очень, видимо для численной реализации...
Цитата:
3. Смогли бы Вы написать расчётные формулы, по которым должна работать программа (не та, что Вы привели, а Ваша, ещё не написанная)?

Аппроксимацию, по которой и строится сетка, я вроде бы усвоил и подогнать для программы могу

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение15.05.2015, 10:38 
Заслуженный участник


23/07/08
10609
Crna Gora
Вообще я в этой теме не специалист, но смысл уравнений, вроде, понятен.
Здесь, похоже, используется такая модель. Есть две среды, «одномерные», с одной координатой $x$ — глубиной (от других координат никакие величины не зависят и потому они в задаче не фигурируют). Среда 1 — промёрзший грунт, среда 2 — непромёрзший. Вся «фишка» задачи в том, что глубина границы раздела $S(t)$, то есть глубина промерзания, зависит от времени. Грубо говоря, если становится теплее, граница поднимается, а если холоднее — опускается. Каждая среда характеризуется набором известных постоянных параметров:
$c$удельная теплоёмкость (массовая), $\frac{\text{Дж}}{\text{кг}\cdot\text{K}}$;
$\rho$ — плотность, $\frac{\text{кг}}{\text{м}^3}$;
$\lambda$удельная теплопроводность, $\frac{\text{Вт}}{\text{м}\cdot\text{K}}$, она же коэффициент теплопроводности.
Кроме того, для описания перехода воды из твердого в жидкое состояние и обратно нужна ещё плотность воды $\rho_w$ и её удельная теплота плавления $q_w$, $\frac{\text{Дж}}{\text{кг}}$. Похоже, в программе ошибка, потому что по таблицам последняя величина равна $3.3477\cdot 10^5 \frac{\text{Дж}}{\text{кг}}$, а в программе q_w=3.34*1000000, то есть в 10 раз больше.
В задаче рассматриваются глубины не до ядра Земли :D , а только до некоторой глубины $\ell=50\text{м}$. Вероятно, ниже этой глубины температура практически постоянна, независимо от времени года.

Помимо констант, в задаче имеется неизвестная функция — температура $T(x, t)$, зависящая не только от координаты $x$, но и от времени $t$. Собственно, в задаче и требуется найти зависимость $T(x, t)$. Как известно, изменение распределения температуры в среде со временем подчиняется уравнению теплопроводности. В Ваших обозначениях оно имеет вид $c\rho\frac{\partial T}{\partial t} = \lambda \frac{\partial^2 T}{\partial x^2}$.
Уравнение, приведённое в Википедии, отличается от Вашего тем, что там 1) константы собраны в одну, 2) имеются дополнительные внутренние источники тепла и 3) рассматривается трехмерный случай.

$c_1 \rho_1\frac{\partial T}{\partial t} = \lambda_1 \frac{\partial^2 T}{\partial x^2}, \quad 0<x<S(t)$
Это уравнение теплопроводности для среды 1 (промёрзший грунт). Среда 1 имеет координаты от $0$ (поверхность земли) до $S(t)$ (глубина промерзания в момент $t$).

$c_2 \rho_2\frac{\partial T}{\partial t} = \lambda_2 \frac{\partial^2 T}{\partial x^2} , \quad S(t)<x<\ell$
Это уравнение теплопроводности для среды 2 (непромёрзший грунт). Среда 2 имеет координаты от $S(t)$ до $\ell$ (дальше не рассматриваем).

$\lambda_1 \frac{\partial T}{\partial x} - \lambda_2 \frac{\partial T}{\partial x} = q_w \rho_w \frac{\partial S}{\partial t}$
Это уравнение описывает смещение границы раздела сред со временем, то есть процесс промерзания и оттаивания грунта. В левой части стоит разность потоков тепла в обеих средах у границы раздела. Эта разница идёт на замораживание или размораживание грунта на границе раздела, что и приводит к смещению границы. По сути, это уравнение баланса энергии на границе. Тепловая энергия, необходимая для перевода 1 килограмма льда в жидкое состояние, известна — это $q_w$, а $\rho_w$ позволяет пересчитать массу в объём.

$\ell = 50 \text{м}$
Выбрали такую предельную глубину, до которой рассматриваются процессы теплопередачи.

$c_1 \rho_1 = 2 \cdot10^6, \lambda_1 = 2$
$c_2 \rho_2 = 2 \cdot 10^6, \lambda_2 = 0,5$
Для каждой среды можно было вообще все константы перенести в правую часть и заменить их одним коэффициентом температуропроводности: $\alpha=\frac{\lambda}{c\rho}$ (как сделано в Википедии), тогда бы уравнения теплопроводности приняли вид $\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}$, где $\alpha$ в каждой из сред своё.

$T_- = T_{\text{Ф}}=T_+, \quad\quad x=S(t)$
На границе раздела сред $x = S(t)$ температура сверху от границы $T_-$ равна температуре снизу от границы $T_+$, и обе совпадают с температурой фазового перехода лёд-вода $T_{\text{Ф}}$.

В задаче заданы начальные и граничные условия.

$T(x,0) = T^0(x), \quad\quad 0\leqslant x \leqslant l,\;\;t=0$
Начальная температура в момент времени $t=0$ в каждой точке из рассматриваемого диапазона координат равна заданной функции $T^0(x)$.

$T(0,t) = T_1(t) = \sin(\pi a t),\quad\quad x=0,\;\;t>0$
Температура в точке $x=0$, т.е. на поверхности земли, в любой момент времени $t$ равна заданной функции $T_1(t)=\sin(\pi at)$. Эта функция, вероятно, описывает годовые колебания температуры: один период синусоиды — один год. У Вас в функции было $\pi a x$, я исправил на $\pi a t$, здесь нужна зависимость от времени.

$T(\ell,t) = T_H,  \;\;\frac{\partial T}{\partial x}(\ell,t) = q, \quad\quad x=\ell,\;\;t>0$
Температура и её производная по координате на нижней границе среды 2 в любой момент времени постоянны и заданы (на такой глубине погода уже не сказывается :-) ).

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение20.05.2015, 10:31 


13/05/15
17
г. Якутск
svv
Спасибо за столь подробное разьяснение постановки :roll:
Вы можете мне сказать зачем вообще было введено ограничение (maxy>1e-4) и что оно означает?

-- 20.05.2015, 16:50 --

Я забыл указать в постановке еще одно уравнение
$(\tilde{c \rho} + \delta(T_\text{Ф} - T) \cdot q_w \cdot \rho_w)\frac{\partial T}{\partial t} = \tilde{\lambda}\frac{\partial T}{\partial x}$
И мне кажется видимо первые два (или три уравнения) судя по программе
Цитата:
$c_1 \rho_1\frac{\partial T}{\partial t} = \lambda_1 \frac{\partial^2 T}{\partial x^2}  , 0<x<S(t)\\
c_2 \rho_2\frac{\partial T}{\partial t} = \lambda_2 \frac{\partial^2 T}{\partial x^2}  , S(t)<x<l\\
\lambda_1 \frac{\partial T}{\partial x} - \lambda_2 \frac{\partial T}{\partial x} = q_w \rho_w \frac{\partial S}{\partial t}\\$

были заменены им для выставления условий:
Код:
if (y[i] < Tfaza) lambda = _lambda(1.0,m);
else lambda = _lambda(0.0, m);
if (y[i] < (Tfaza - delta)) cro = c_ro(1.0,m);
else if (y[i] < (Tfaza + delta)) cro = c_ro(0.0,m);
else cro = c_ro(1.0,m) + m*q_w*ro_w / (2 * delta);


-- 20.05.2015, 16:52 --

либо я что-то путаю

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение20.05.2015, 10:58 
Заслуженный участник


23/07/08
10609
Crna Gora
Программу я ещё не смотрел, чуть позже посмотрю и постараюсь ответить.

По-хорошему, между дифференциальными уравнениями и программой должен быть ещё один этап — расчётные формулы, или запись алгоритма «почти как в будущей программе», но в математических обозначениях. Чтобы было понятно — вот пример. (Даже всю ту тему бегло просмотрите.)

Если этот этап выполнен, программу очень легко будет и самому написать. Если уже есть чужая программа, её будет легче анализировать, если в такой вид перевести.

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение20.05.2015, 11:20 


13/05/15
17
г. Якутск
svv в сообщении #1017719 писал(а):
По-хорошему, между дифференциальными уравнениями и программой должен быть ещё один этап — расчётные формулы

Эти формулы у меня уже есть, программа на них вроде нормально работает, но с какими-то ошибками(вероятнее всего в вычислении maxy)
Используется неявная схема(она мне больше нравится ибо там не надо проверять на устойчивость)
$c \rho \frac{y_i}{\tilde{y_i}}{\tau} = \lambda \frac{y_{i+1}+2y_{i}+y_{i-1}}{h^2}\\
y_i(\frac{\tilde{c \rho}}{\tau} + \frac{2 \tilde{\lambda}}{h^2}) = (\frac{c \rho \tilde{y_i}}{\tau} + y_{i+1}\frac{\lambda}{h^2}+y_i-1\frac{\lambda}{h^2})
Используется метод прогонки, т.е. $a_i {y__{i-1}} - c_i y_i + b_i y_{i+1} +d_i$
Вы про этот этап имели ввиду?

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение20.05.2015, 11:31 
Заслуженный участник


23/07/08
10609
Crna Gora
Да, про этот. :-)
Пожалуйста, подождите день-два, пока я посмотрю Вашу программу.

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение20.05.2015, 11:36 


13/05/15
17
г. Якутск
svv
Хорошо, Спасибо.

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение21.05.2015, 18:59 
Заслуженный участник


23/07/08
10609
Crna Gora
Предварительно: я нашел в программе несколько ошибок, из-за которых она вряд ли будет работать. Вы можете аргументированно возражать :-) . Для удобства я их пронумерую.

1. Согласно таблицам, для воды (льда) удельная теплота плавления воды равна $3.35\cdot 10^5 \frac{\text{Дж}}{\text{кг}}$, а в программе в 10 раз больше.

2. Дальше обратите внимание на формулы для $c\rho$ и $\lambda$. Пусть индекс $s$ означает грунт сам по себе, $w$ — вода в жидком состоянии, $i$ — лёд. Пусть $M$ — пористость (т.е. доля воды в любом состоянии в общем количестве вещества), $S_w$ — водонасыщенность (доля воды в жидком состоянии в общем количестве воды). Тогда, действительно, считая, что каждая фракция вносит вклад пропорционально своей доле, получаем формулы вроде такой:
$\lambda=(1 - M)\lambda_s +M S_w \lambda_w  + M (1-S_w) \lambda_i$
Именно так в программе. Здесь всё правильно и логично. Первое слагаемое — доля грунта в веществе, т.е. $1-M$. Грунт не делится на твердый и жидкий. Доля воды $M$. Доля жидкой воды в веществе $MS_w$ (второе слагаемое), а доля твердой воды $M(1-S_w)$ (третье слагаемое).
Чтобы лучше почувствовать написанную формулу для $\lambda$, посмотрите, что она даёт в трёх случаях:
$\bullet$ $M=0,\; S_w$ любое
$\bullet$ $M=1,\; S_w=1$
$\bullet$ $M=1,\; S_w=0$

Понятно, что $M$ от температуры не зависит: если в грунте есть вода, изменение температуры может только перевести её в другое агрегатное состояние, но не может изменить долю воды. А вот $S_w$ от температуры зависит: когда очень тепло, вся вода жидкая ($S_w=1$), когда очень холодно, вся вода твёрдая ($S_w=0$). Так должно быть.

Но теперь посмотрите, как в программе вычисляется «суммарная» $\lambda$ грунта:
if (y[i] < Tfaza) lambda = _lambda(1.0,m);
else lambda = _lambda(0.0, m);

и аналогично вычисляется «суммарное» $c\rho$. Всё наоборот. Когда температура ниже фазового перехода, используется водонасыщенность $1$, иначе $0$.

3. Функции _lambda и c_ro вычисляют число типа double (с плавающей точкой, двойной точности), это правильно. Но поскольку тип возвращаемого значения у них обеих int, заботливо вычисленная дробная часть отбрасывается и остаётся только целая часть. :shock: :shock:

4. Посмотрите, как вычисляется $c\rho$:
if (y[i] < (Tfaza - delta)) cro = c_ro(1.0,m);
else if (y[i] < (Tfaza + delta)) cro = c_ro(0.0,m);
else cro = c_ro(1.0,m) + m*q_w*ro_w / (2 * delta);

Первая альтернатива должна сработать, когда температура определённо ниже фазового перехода ($T<T_{\text{Ф}}-\delta$). Про неправильную водонасыщенность я уже сказал в п.2, но знак < правильный.
Вторая альтернатива должна сработать, когда температура определённо выше фазового перехода ($T>T_{\text{Ф}}+\delta$). Во второй строке кода знак < неправильный.
Третья альтернатива должна сработать в остальных случаях, когда температура с точностью до небольшого зазора $\pm\delta$ равна температуре фазового перехода $T_{\text{Ф}}$. Здесь теплоёмкость записана с дополнительным слагаемым, которое представляет собой удельную теплоту плавления, «размазанную» по зазору $T_{\text{Ф}}-\delta<T<T_{\text{Ф}}+\delta$. Это я могу понять.

Исправьте все эти ошибки и покажите, что получилось. Скоро дойдём и до maxy.

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение23.05.2015, 06:28 


13/05/15
17
г. Якутск
svv
Я не знаю чему изначально будет равна $\tau$ и $t$ и поэтому приравнял от балды, как и $T_{gran}$
Вот исправленная программа

код: [ скачать ] [ спрятать ]
Используется синтаксис 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 100
#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;
        for (j = 0; j <= N; j++)
        {
                y[j] = Tnach;
                yj[j] = y[j];
        }
        for (j = 1; j < 100; 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];
                                c = a + b + c;
                                alpha[i] = b / (c - a*alpha[i - 1]);
                                beta[i] = (d + a*beta[i - 1]) / (c - 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 < abs(y[i] - a)) maxy = abs(y[i] - a);
                        }
//                      printf("%d::",k);
                        for (int i = 1; i <= n; i++) printf("%.5f |",y[i]);
                        printf("\n\n");
                        yj[j] = y[j];
                }
        }
}

и результат практически тот же самый
https://yadi.sk/i/m6XXXzOOgpHj6

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение23.05.2015, 09:43 


13/05/15
17
г. Якутск
svv
И еще одно замечание по вашим исправлениям:
Цитата:
Понятно, что $M$ от температуры не зависит: если в грунте есть вода, изменение температуры может только перевести её в другое агрегатное состояние, но не может изменить долю воды.

Пористость у меня по условию равна 0,4.

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение26.05.2015, 16:16 
Заслуженный участник


23/07/08
10609
Crna Gora
В своих попытках понять логику работы программы (ох, легче свою с нуля написать!) я привёл программу к такому виду:
код: [ скачать ] [ спрятать ]
Используется синтаксис C++
const int
   n=100,
   N=100;
const double
   Pi=M_PI,
   q_w=3.34e5,
   delta=0.01,
   h=0.1,
   c_s=700,
   ro_s=4000,
   c_w=4190,
   ro_w=1000,
   c_ice=2000,
   ro_ice=900,
   lambda_s=900,
   lambda_w=900,
   lambda_ice=900,
   m=0.4,
   Tfaza=0,
   Tnach=-3;

double get_lambda(double T)
{
   double S_w = T<Tfaza? 0: 1;
   return (1-m)*lambda_s + lambda_w*S_w*m + lambda_ice*(1-S_w)*m;
}

double get_cro(double T)
{
   double
      S_w = (T <= Tfaza+delta)? 0: 1,
      c_ro = (1-m)*c_s*ro_s + c_w*ro_w*S_w*m + c_ice*ro_ice*(1-S_w)*m;
   if (fabs(T - Tfaza)<=delta) c_ro+=m*q_w*ro_w / (2 * delta);
   return c_ro;
}

double Tgran(double a, double t)
{
   return sin(Pi*a*t);
}

double max(double a, double b)
{
   return a>b? a: b;
}

double min(double a, double b)
{
   return a<b? a: b;
}

void Calc()
{
   double t = 0, maxy = 1, y[n+1], yj[n+1];
   for (int j = 0; j <= N; j++)
   {
      y [j] = Tnach;
      yj[j] = Tnach;
   }
   for (int j = 1; j < 100; j++)
   {
      double tau = min(pow(1.1, j), 1e6);
      t = t + tau;
      for (int k = 1; k <= 20 && maxy > 1e-4; k++)
      {
         double alpha[n+1], beta[n+1];

         alpha[0] = 0;
         beta [0] = Tgran(1,t);
         for (int i = 1; i < n; i++)
         {
            double
               lambda = get_lambda(y[i]),
               cro = get_cro(y[i]),
               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] = beta[n - 1] / (1 - alpha[n - 1]);
         double maxy = y[n];
         for (int i = n - 1; i >= 0; i--)
         {
            double a = y[i];
            y[i] = y[i + 1] * alpha[i] + beta[i];
            maxy = max(maxy, abs(y[i] - a));
         }
         yj[j] = y[j];
      }
   }
}
 

Вам не нужно делать так, как у меня. Это только для пояснения хода мыслей. Здесь язык C++, а у Вас — C или совместимое с ним подмножество C++. Я постарался сделать область видимости переменных как можно более узкой, чтобы легче было анализировать их поведение.

В строке
if (maxy < abs(y[i] - a)) maxy = abs(y[i] - a);
вместо abs точно надо поставить fabs, иначе будет преобразование к целому.

Много непонятного. Самый внешний цикл, с управляющей переменной j, очевидно, по дискретному времени. Почему с каждым следующим j шаг по времени tau растёт, вплоть до астрономической величины $10^6$?
Группа операторов, вычисляющих a, b, c, d и затем alpha[i], beta[i], малопонятна. Не стоит одну и ту же переменную (c) использовать в двух разных смыслах, это источник ошибок, лучше завести дополнительную (у меня e). Всю эту группу надо тщательно проверять.
Совершенно непонятен оператор yj[j] = y[j];. Почему надо копировать в yj только одно значение, почему именно это, и почему это надо делать при каждом k? Можно ведь сделать это один раз, после цикла по $k$, когда элемент y[j] примет окончательное для данного j, «устаканившееся» значение.

Что касается maxy. Это максимальное значение модуля производной от y[i] (т.е. температуры). Вот только пространственной или временной производной — Вам виднее. Когда этот максимальный модуль производной становится меньше $10^{-4}$, итерации (т.е. цикл по k) можно завершить, не дожидаясь k=20 — предельного количества итераций.

Похоже, дальше сможете продвинуться только Вы сами. :P

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение26.05.2015, 18:18 


13/05/15
17
г. Якутск
svv
Спасибо большое, за столь подробный и полный ответ. Мне многое наконец понятно Благодаря Вам.
svv писал(а):
Много непонятного.
Совершенно непонятен оператор yj[j] = y[j];. Почему надо копировать в yj только одно значение, почему именно это, и почему это надо делать при каждом k? Можно ведь сделать это один раз, после цикла по $k$, когда элемент y[j] примет окончательное для данного j, «устаканившееся» значение.

Видимо это надо для того чтобы из каждого вычисленного слоя сделать "исходный" слой, от которого можно вычислить последующий.
Насчет fabc спасибо теперь у меня есть хоть какие-то нормальные результаты. Однако при изменении константы n у меня выводит ошибку сегментации :(
И Спасибо за разъяснение про maxy, но меня волнует не столь его значение, сколько вопрос "для чего он тут? и для чего сделан цикл с maxy и k?".

 Профиль  
                  
 
 Re: Промерзание грунта
Сообщение26.05.2015, 21:25 
Заслуженный участник


23/07/08
10609
Crna Gora
По всей видимости, в этой программе уравнение решается методом последовательных приближений. Процесс последовательных приближений должен быстро сходиться к точному решению уравнения. Для этого в программе имеется цикл по k — это номер последовательного приближения (итерации). Этот цикл не должен быть бесконечным. Из него существует два выхода:
1) выход «по-хорошему» — достигли необходимой точности, а maxy — индикатор этой точности: когда он становится достаточно малым (1e-4), процесс можно прекращать; почему именно maxy является индикатором — другой вопрос;
2) выход «по-плохому», когда совершили уже 20 итераций, а процесс почему-то не сходится; в такой ситуации ясно, что «пора кончать», несмотря на неудачу — и для этого проверяется условие k<=20.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 29 ]  На страницу 1, 2  След.

Модераторы: Модераторы Математики, Супермодераторы



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

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


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group