forked from tomflint22/MHD_Analytical
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmain.cpp
157 lines (88 loc) · 4.54 KB
/
main.cpp
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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
#include <iostream>
#include <math.h>
#include <fstream>
using namespace std;
double alpha_k(double k, double l);
double r1k(double k, double l, double B, double L, double sigma, double density, double k_viscosity);
double r2k(double k, double l, double B, double L, double sigma, double density, double k_viscosity);
double N(double k, double l, double B, double L, double sigma, double density, double k_viscosity);
double V2(double k, double l, double B, double L, double sigma, double density, double k_viscosity, double xi, double eta);
double V3(double k, double l, double B, double L, double sigma, double density, double k_viscosity, double xi, double eta);
double past_sum(double k, double l, double B, double L, double sigma, double density, double k_viscosity, double xi, double eta);
double sum(double l, double B, double L, double sigma, double density, double k_viscosity, double xi, double eta);
int main()
{
ofstream myfile;
myfile.open ("MHD_Analytical.txt");
double B =10.0;//0.00052552761542017795;
double L =1.0;
double mu =0.01;//1.2566156992899439E-06;//Permeability
double sigma =100.0;//4.55e6;//Electrical conductivity
double pressure_gradient = 3768.4/20.0;
double kinematic_viscosity = 1.0;//1.322401481089659E-07;
double density = 1.0;//11343;//
double asqr=1.0*1.0;
double Ha = B*L*sqrt(sigma/(density*kinematic_viscosity));
cout<<"Hartmann Number = "<<Ha<<endl;
for(double dist=0.0;dist<=1.01;dist+=0.01){
// cout<<dist<<"\t"<<(1.0/(kinematic_viscosity*density))*pressure_gradient*asqr*sum(1.0,B,L,sigma,mu,dist,0.0)<<endl;
myfile<<dist<<"\t"<<(1.0/(density*kinematic_viscosity))*pressure_gradient*asqr*sum(1.0,B,L,sigma,density,kinematic_viscosity,dist,0.0)<<endl;
}
myfile.close();
return 0;
}
double alpha_k(double k, double l){
double alpha_k_out = (k+0.5)*(M_PI/l);
return alpha_k_out;
}
double r1k(double k, double l, double B, double L, double sigma, double density, double k_viscosity){
double Ha = B*L*sqrt(sigma/(density*k_viscosity));
//double r1k_out = 0.5*(Ha+((Ha*Ha)+(4.0*pow(alpha_k(k,l),2.0))));
double r1k_out = 0.5*(Ha+N(k,l,B,L,sigma,density,k_viscosity));
return r1k_out;
}
double r2k(double k, double l, double B, double L, double sigma, double density, double k_viscosity){
double Ha = B*L*sqrt(sigma/(density*k_viscosity));
//double r2k_out = 0.5*(-Ha+((Ha*Ha)+(4.0*pow(alpha_k(k,l),2.0))));
double r2k_out = 0.5*(-Ha+N(k,l,B,L,sigma,density,k_viscosity));
return r2k_out;
}
double N(double k, double l, double B, double L, double sigma, double density, double k_viscosity){
double Ha = B*L*sqrt(sigma/(density*k_viscosity));
double N_out = sqrt((Ha*Ha)+(4.0*(pow(alpha_k(k,l),2.0))));//0.5*(-Ha+((Ha*Ha)+(4.0*alpha_k(k,l))));
return N_out;
}
double V2(double k, double l, double B, double L, double sigma, double density, double k_viscosity, double xi, double eta){
double db=1e99;
double N_=N(k,l,B,L,sigma,density,k_viscosity);
double r1k_=r1k(k,l,B,L,sigma,density,k_viscosity);
double r2k_=r2k(k,l,B,L,sigma,density,k_viscosity);
double numerator=((db*r2k_)+((1.0-exp(-2.0*r2k_))/(1.0+exp(-2.0*r2k_))))*(0.5*(exp(-r1k_*(1.0-eta))+exp(-r1k_*(1.0+eta))));
double denominator=(((1.0+exp(-2.0*r1k_))/2.0)*(db*N_))+((1.0+exp(-2.0*(r1k_+r2k_)))/(1.0+exp(-2.0*r2k_)));
double V2_out=(numerator)/(denominator);
return V2_out;
}
double V3(double k, double l, double B, double L, double sigma, double density, double k_viscosity, double xi, double eta){
double db=1e99;
double N_=N(k,l,B,L,sigma,density,k_viscosity);
double r1k_=r1k(k,l,B,L,sigma,density,k_viscosity);
double r2k_=r2k(k,l,B,L,sigma,density,k_viscosity);
double numerator=((db*r1k_)+((1.0-exp(-2.0*r1k_))/(1.0+exp(-2.0*r1k_))))*(0.5*(exp(-r2k_*(1.0-eta))+exp(-r2k_*(1.0+eta))));
double denominator=(((1.0+exp(-2.0*r2k_))/2.0)*(db*N_))+((1.0+exp(-2.0*(r1k_+r2k_)))/(1.0+exp(-2.0*r1k_)));
double V3_out=(numerator)/(denominator);
return V3_out;
}
double past_sum(double k, double l, double B, double L, double sigma, double density, double k_viscosity, double xi, double eta){
double alpha_k_=alpha_k(k,l);
double V2_ = V2(k,l,B,L,sigma,density, k_viscosity,xi,eta);
double V3_ = V3(k,l,B,L,sigma,density, k_viscosity,xi,eta);
double value = ((2.0*(pow(-1.0,k))*cos(alpha_k_*xi))/(l*pow(alpha_k_,3.0)))*(1.0-V2_-V3_);
return value;
}
double sum(double l, double B, double L, double sigma, double density, double k_viscosity, double xi, double eta){
double valuesum = 0.0;
for(double k=0.0;k<50000.0;k+=1.0){
valuesum+=past_sum(k,l,B,L,sigma,density,k_viscosity,xi,eta);
}
return valuesum;
}