MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
conductivity.f90
Go to the documentation of this file.
1 
7 module conductivity
8  use constants
9  use fmodena
10  use physicalproperties
11  implicit none
12  private
13  public equcond
14 contains
15 !********************************BEGINNING*************************************
20 subroutine equcond(keq,ystate,ngas,nfv,mor,eps,dcell,fstrut,temp)
21  real(dp), intent(out) :: keq
22  real(dp), dimension(:), intent(in) :: ystate
23  integer, intent(in) :: ngas
24  integer, intent(in) :: nfv
25  integer, intent(in) :: mor(:)
26  real(dp), intent(in) :: temp
27  real(dp), intent(in) :: eps
28  real(dp), intent(in) :: dcell
29  real(dp), intent(in) :: fstrut
30  real(dp) :: kgas
31  real(dp), dimension(ngas) :: kg,yg,cc
32  integer :: i,j,k
33  integer :: gasModel=2
34  !calculate average concentrations
35  cpg(1)=oxyheatcapacity(temp)
36  cpg(2)=nitrheatcapacity(temp)
37  cpg(3)=cdheatcapacity(temp)
38  cpg(4)=cypheatcapacity(temp)
39  cc=0
40  j=0
41  do i=1,nfv
42  if (mor(i)==1) then
43  do k=1,ngas
44  cc(k)=cc(k)+ystate(ngas*(i-1)+k)
45  enddo
46  j=j+1
47  endif
48  enddo
49  cc=cc/j
50  yg=cc/sum(cc)
51  do i=1,ngas
52  kg(i)=gasconductivity(temp,i)
53  if (abs(yg(i))<1e-6_dp) then
54  yg(i)=0
55  endif
56  enddo
57  select case(gasmodel) ! determine conductivity of gas mixture
58  case (1)
59  kgas = weightedaverage(kg,yg)
60  case (2)
61  kgas = extwassiljewa(kg,yg,tc,pc,mg,temp,1._dp)
62  case (3)
63  kgas = lindsaybromley(kg,yg,tb,cpg,mg,temp)
64  case (4)
65  kgas = pandeyprajapati(kg,yg,tb,mg,temp)
66  case default
67  stop 'unknown gas conductivity model'
68  end select
69  call modena_inputs_set(kfoaminputs, kfoamepspos, eps)
70  call modena_inputs_set(kfoaminputs, kfoamdcellpos, dcell)
71  call modena_inputs_set(kfoaminputs, kfoamfstrutpos, fstrut)
72  call modena_inputs_set(kfoaminputs, kfoamtemppos, temp)
73  do i=1,ngas
74  call modena_inputs_set(kfoaminputs, kfoamxg(i), yg(i))
75  enddo
76  ret = modena_model_call(kfoammodena, kfoaminputs, kfoamoutputs)
77  if (modena_error_occurred()) then
78  call exit(modena_error())
79  endif
80  keq = modena_outputs_get(kfoamoutputs, 0_c_size_t); !fetch results
81 end subroutine equcond
82 !***********************************END****************************************
83 
84 
85 !********************************BEGINNING*************************************
89 real(dp) function weightedaverage(k,yin) result(kmix)
90  real(dp), dimension(:), intent(in) :: k !thermal conductivities
91  real(dp), dimension(:), intent(in) :: yin !molar fractions
92  integer :: n,i,j
93  real(dp), dimension(:), allocatable :: y
94  n=size(k)
95  allocate(y(n))
96  y=yin
97  if (minval(y)<0) then
98  print*, 'Input molar fractions to weightedAverage cannot be negative.'
99  print*, y
100  stop
101  endif
102  y=y/sum(y)
103  kmix=0
104  do i=1,size(k)
105  kmix=kmix+yin(i)*k(i)
106  enddo
107 end function weightedaverage
108 !***********************************END****************************************
109 
110 
111 !********************************BEGINNING*************************************
117 real(dp) function extwassiljewa(k,yin,Tc,pc,M,T,eps) result(kmix)
118  real(dp), dimension(:), intent(in) :: k !thermal conductivities
119  real(dp), dimension(:), intent(in) :: yin !molar fractions
120  real(dp), dimension(:), intent(in) :: Tc !critical temperatures
121  real(dp), dimension(:), intent(in) :: pc !critical pressures
122  real(dp), dimension(:), intent(in) :: M !molar masses
123  real(dp), intent(in) :: T !temperature
124  real(dp), intent(in) :: eps !parameter close to one
125  integer :: n,i,j
126  real(dp) :: x
127  real(dp), dimension(:), allocatable :: gam,y
128  real(dp), dimension(:,:), allocatable :: ktr,A
129  n=size(k)
130  allocate(y(n),gam(n),ktr(n,n),a(n,n))
131  y=yin
132  if (minval(y)<0) then
133  print*, 'Input molar fractions to extWassiljewa cannot be negative.'
134  print*, y
135  stop
136  endif
137  y=y/sum(y)
138  gam=210*(tc*m**3/pc**4)**(1/6._dp)
139  do i=1,n
140  do j=1,n
141  ktr(i,j)=gam(j)*(exp(0.0464_dp*t/tc(i))-exp(-0.2412_dp*t/tc(i)))/&
142  gam(i)/(exp(0.0464_dp*t/tc(j))-exp(-0.2412_dp*t/tc(j)))
143  a(i,j)=eps*(1+sqrt(ktr(i,j))*(m(i)/m(j))**0.25_dp)**2/&
144  sqrt(8*(1+m(i)/m(j)))
145  enddo
146  enddo
147  kmix=0
148  do i=1,n
149  x=0
150  do j=1,n
151  x=x+y(j)*a(i,j)
152  enddo
153  kmix=kmix+y(i)*k(i)/x
154  enddo
155 end function extwassiljewa
156 !***********************************END****************************************
157 
158 
159 !********************************BEGINNING*************************************
164 real(dp) function lindsaybromley(k,yin,Tb,cp,M,T) result(kmix)
165  real(dp), dimension(:), intent(in) :: &
166  k,& !thermal conductivities
167  yin,& !molar fractions
168  Tb,& !boiling point temperatures
169  cp,& !thermal capacities at constant pressure
170  M !molar masses
171  real(dp), intent(in) :: &
172  T !temperature
173  integer :: n,i,j
174  real(dp) :: x
175  real(dp), dimension(:), allocatable :: &
176  y,& !molar fractions
177  S,& !Sutherland constants
178  gam,& !heat capacity ratio
179  cv !thermal capacities at constant volume
180  real(dp), dimension(:,:), allocatable :: A
181  n=size(k)
182  allocate(y(n),cv(n),s(n),gam(n),a(n,n))
183  y=yin
184  if (minval(y)<0) then
185  print*, 'Input molar fractions to lindsayBromley cannot be negative.'
186  print*, y
187  stop
188  endif
189  y=y/sum(y)
190  do i=1,n
191  s(i)=1.5_dp*tb(i)
192  enddo
193  cv=cp-rg !ideal gas assumed
194  gam=cp/cv
195  do i=1,n
196  do j=1,n
197  x=k(i)/k(j)*cp(j)/cp(i)*(9-5/gam(j))/(9-5/gam(i))
198  a(i,j)=&
199  0.25_dp*(1+sqrt(x*(m(j)/m(i))**0.75_dp*(t+s(i))/(t+s(j))))**2*&
200  (t+sqrt(s(i)*s(j)))/(t+s(j))
201  enddo
202  enddo
203  kmix=0
204  do i=1,n
205  x=0
206  do j=1,n
207  x=x+y(j)*a(i,j)
208  enddo
209  kmix=kmix+y(i)*k(i)/x
210  enddo
211 end function lindsaybromley
212 !***********************************END****************************************
213 
214 
215 !********************************BEGINNING*************************************
220 real(dp) function pandeyprajapati(k,yin,Tb,M,T) result(kmix)
221  real(dp), dimension(:), intent(in) :: &
222  k,& !thermal conductivities
223  yin,& !molar fractions
224  Tb,& !boiling point temperatures
225  M !molar masses
226  real(dp), intent(in) :: &
227  T !temperature
228  integer :: n,i,j
229  real(dp) :: x
230  real(dp), dimension(:), allocatable :: &
231  y,& !molar fractions
232  S !Sutherland constants
233  real(dp), dimension(:,:), allocatable :: A
234  n=size(k)
235  allocate(y(n),s(n),a(n,n))
236  y=yin
237  if (minval(y)<0) then
238  print*, 'Input molar fractions to pandeyPrajapati cannot be negative.'
239  print*, y
240  stop
241  endif
242  y=y/sum(y)
243  do i=1,n
244  s(i)=1.5_dp*tb(i)
245  enddo
246  do i=1,n
247  do j=1,n
248  a(i,j)=0.25_dp*(1+sqrt(k(i)/k(j)*(m(i)/m(j))**0.25_dp*&
249  (t+s(i))/(t+s(j))))**2*(t+sqrt(s(i)*s(j)))/(t+s(j))
250  enddo
251  enddo
252  kmix=0
253  do i=1,n
254  x=0
255  do j=1,n
256  x=x+y(j)*a(i,j)
257  enddo
258  kmix=kmix+y(i)*k(i)/x
259  enddo
260 end function pandeyprajapati
261 !***********************************END****************************************
262 end module conductivity