sagd/filtration.cpp at master · CompMathGroup/sagd · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include "const.h"
#include <math.h>

double theta_l(const double psi)
{
	return ( psi / (1 + psi) );

}

double theta_s(const double psi)
{
	return ( 1 / (1 + psi) );
}

void K_permeability(double* psi, double* K_perm) //Проницаемость
{
	int m;
	for (m=0;m<M;m++)
	{
		K_perm[m] = 4 * pow( psi[m], 3) / pow(1 + psi[m], 2) / (2 + psi[m]) ;
		//K_perm[m] = pow( theta_l( psi[m] ), 2);   //Задаю проницаемость как квадрат насыщенности.
	}															//(На самом деле в тетта лежит отношение насыщеносттей, поэтому такая формула)
}


void viscosity(double * eta, const double* T)
{
	int m;
	for (m = 1; m < M; m++)
	{
		eta[m] = 0.02 * exp( alpha * (Tcrit - T[m]) );
		//eta[m] = 1e-3;
	}
}

void filtr_satur(const double *K_perm, const double *psi, double *psi_new, double *W_fil, const double tau, double* a_max, const double* eta)
{
	int m;
	double psi_temporary;
	*a_max = 0;
	for (m = 1; m < M; m++)
	{
		W_fil[m]=K_abs*K_perm[m] / eta[m] * (ro_s - ro_l) * gravity;
		psi_temporary = psi[m] - tau * (W_fil[m] - W_fil[m-1])/h;

		if (2 * psi_temporary / pow((1+psi_temporary), 3) / eta[m] > *a_max ) //Вычисляем максимальное a для условия Куранта
		{
			*a_max = 2 * psi_temporary / pow((1+psi_temporary) , 3) / eta[m];
		}

		if (psi_temporary >= psi_crit)
		{
			psi_new[m] = psi_temporary;
			continue;
		}
		if (psi_temporary < psi_crit)
		{
			psi_new[m] = psi_crit;
			W_fil[m] = W_fil[m-1] + (psi[m] - psi_new[m]) * h / tau;
		}

		if (2 * psi_new[m] / pow((1+psi_new[m]), 3) / eta[m]  > *a_max ) 	//Вычисляем максимальное a для условия Куранта
		{
			*a_max = 2 * psi_new[m] / pow((1+psi_new[m]) , 3) / eta[m];
		}
	}

}