MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
physicalProperties.f90
Go to the documentation of this file.
1 
8 module physicalproperties
9  use constants
10  use globals, only: solmodel,diffmodel
11  use fmodena
12  implicit none
13  !modena variables
14  integer(c_int) :: ret
15  type(c_ptr) :: rhopmodena = c_null_ptr
16  type(c_ptr) :: rhopinputs = c_null_ptr
17  type(c_ptr) :: rhopoutputs = c_null_ptr
18  integer(c_size_t) :: rhoptemppos
19  type(c_ptr) :: kfoammodena = c_null_ptr
20  type(c_ptr) :: kfoaminputs = c_null_ptr
21  type(c_ptr) :: kfoamoutputs = c_null_ptr
22  integer(c_size_t) :: kfoamepspos
23  integer(c_size_t) :: kfoamdcellpos
24  integer(c_size_t) :: kfoamfstrutpos
25  integer(c_size_t) :: kfoamtemppos
26  integer(c_size_t), allocatable :: kfoamxg(:)
27  type(c_ptr) :: kgasmodena = c_null_ptr
28  type(c_ptr) :: kgasinputs = c_null_ptr
29  type(c_ptr) :: kgasoutputs = c_null_ptr
30  integer(c_size_t) :: kgastemppos
31  integer(c_size_t), allocatable :: kgasxg(:)
32  type(c_ptr), dimension(:), allocatable :: sgmodena
33  type(c_ptr), dimension(:), allocatable :: sginputs
34  type(c_ptr), dimension(:), allocatable :: sgoutputs
35  integer(c_size_t), dimension(:), allocatable :: sgtemppos
36  integer(c_size_t), dimension(:), allocatable :: sgxl1pos
37  integer(c_size_t), dimension(:), allocatable :: sgxl2pos
38  type(c_ptr), dimension(:), allocatable :: dgmodena
39  type(c_ptr), dimension(:), allocatable :: dginputs
40  type(c_ptr), dimension(:), allocatable :: dgoutputs
41  integer(c_size_t), dimension(:), allocatable :: dgtemppos
42  type(c_ptr), dimension(:), allocatable :: kgmodena
43  type(c_ptr), dimension(:), allocatable :: kginputs
44  type(c_ptr), dimension(:), allocatable :: kgoutputs
45  integer(c_size_t), dimension(:), allocatable :: kgtemppos
46 contains
47 !********************************BEGINNING*************************************
51 subroutine createmodels(ngas)
52  integer :: ngas
53  integer :: i
54  kfoammodena = modena_model_new(c_char_"foamConductivity"//c_null_char);
55  if (modena_error_occurred()) then
56  call exit(modena_error())
57  endif
58  kfoaminputs = modena_inputs_new(kfoammodena);
59  kfoamoutputs = modena_outputs_new(kfoammodena);
60  kfoamepspos = modena_model_inputs_argpos(&
61  kfoammodena, c_char_"eps"//c_null_char);
62  kfoamdcellpos = modena_model_inputs_argpos(&
63  kfoammodena, c_char_"dcell"//c_null_char);
64  kfoamfstrutpos = modena_model_inputs_argpos(&
65  kfoammodena, c_char_"fstrut"//c_null_char);
66  kfoamtemppos = modena_model_inputs_argpos(&
67  kfoammodena, c_char_"T"//c_null_char);
68  do i=1,ngas
69  kfoamxg(i) = modena_model_inputs_argpos(kfoammodena, &
70  c_char_"x["//trim(adjustl(gasname(i)))//"]"//c_null_char);
71  enddo
72  call modena_model_argpos_check(kfoammodena)
73  kgasmodena = modena_model_new(&
74  c_char_"gasMixtureConductivity"//c_null_char)
75  if (modena_error_occurred()) then
76  call exit(modena_error())
77  endif
78  kgasinputs = modena_inputs_new(kgasmodena)
79  kgasoutputs = modena_outputs_new(kgasmodena)
80  kgastemppos = modena_model_inputs_argpos(&
81  kgasmodena, c_char_"T"//c_null_char)
82  do i=1,ngas
83  kgasxg(i) = modena_model_inputs_argpos(kgasmodena, &
84  c_char_"x["//trim(adjustl(gasname(i)))//"]"//c_null_char);
85  enddo
86  call modena_model_argpos_check(kgasmodena)
87  do i=1,ngas
88  kgmodena(i) = modena_model_new(c_char_&
89  "gas_thermal_conductivity[A="//trim(adjustl(gasname(i)))//"]"//&
90  c_null_char);
91  if (modena_error_occurred()) then
92  call exit(modena_error())
93  endif
94  kginputs(i) = modena_inputs_new(kgmodena(i));
95  kgoutputs(i) = modena_outputs_new(kgmodena(i));
96  kgtemppos(i) = modena_model_inputs_argpos(&
97  kgmodena(i), c_char_"T"//c_null_char);
98  call modena_model_argpos_check(kgmodena(i))
99  if (solmodel(i)==1) then
100  sgmodena(i) = modena_model_new(c_char_&
101  "Solubility[A="//trim(adjustl(gasname(i)))//",B=2]"//&
102  c_null_char);
103  if (modena_error_occurred()) then
104  call exit(modena_error())
105  endif
106  sginputs(i) = modena_inputs_new(sgmodena(i));
107  sgoutputs(i) = modena_outputs_new(sgmodena(i));
108  sgtemppos(i) = modena_model_inputs_argpos(&
109  sgmodena(i), c_char_"T"//c_null_char);
110  sgxl1pos(i) = modena_model_inputs_argpos(&
111  sgmodena(i), c_char_"xl1"//c_null_char);
112  sgxl2pos(i) = modena_model_inputs_argpos(&
113  sgmodena(i), c_char_"xl2"//c_null_char);
114  call modena_model_argpos_check(sgmodena(i))
115  endif
116  if (diffmodel(i)==1) then
117  dgmodena(i) = modena_model_new(c_char_&
118  "diffusivityPol[A="//trim(adjustl(gasname(i)))//"]"//&
119  c_null_char);
120  if (modena_error_occurred()) then
121  call exit(modena_error())
122  endif
123  dginputs(i) = modena_inputs_new(dgmodena(i));
124  dgoutputs(i) = modena_outputs_new(dgmodena(i));
125  dgtemppos(i) = modena_model_inputs_argpos(&
126  dgmodena(i), c_char_"T"//c_null_char);
127  call modena_model_argpos_check(dgmodena(i))
128  endif
129  enddo
130 end subroutine createmodels
131 !***********************************END****************************************
132 
133 
134 !********************************BEGINNING*************************************
138 subroutine destroymodels(ngas)
139  integer, intent(in) :: ngas
140  integer :: i
141  call modena_inputs_destroy (kfoaminputs);
142  call modena_outputs_destroy (kfoamoutputs);
143  call modena_model_destroy (kfoammodena);
144  call modena_inputs_destroy (kgasinputs);
145  call modena_outputs_destroy (kgasoutputs);
146  call modena_model_destroy (kgasmodena);
147  do i=1,ngas
148  call modena_inputs_destroy (kginputs(i));
149  call modena_outputs_destroy (kgoutputs(i));
150  call modena_model_destroy (kgmodena(i));
151  if (solmodel(i)==1) then
152  call modena_inputs_destroy (sginputs(i));
153  call modena_outputs_destroy (sgoutputs(i));
154  call modena_model_destroy (sgmodena(i));
155  endif
156  if (diffmodel(i)==1) then
157  call modena_inputs_destroy (dginputs(i));
158  call modena_outputs_destroy (dgoutputs(i));
159  call modena_model_destroy (dgmodena(i));
160  endif
161  enddo
162 end subroutine destroymodels
163 !***********************************END****************************************
164 
165 
166 !********************************BEGINNING*************************************
170 real(dp) function polymerdensity(temp)
171  real(dp), intent(in) :: temp
172  call modena_inputs_set(rhopinputs, rhoptemppos, temp)
173  ret = modena_model_call(rhopmodena, rhopinputs, rhopoutputs)
174  if(ret /= 0) then
175  call exit(ret)
176  endif
177  polymerdensity=modena_outputs_get(rhopoutputs, 0_c_size_t)
178 end function polymerdensity
179 !***********************************END****************************************
180 
181 
182 !********************************BEGINNING*************************************
186 real(dp) function gasmixtureconductivity(temp,xg,ngas)
187  integer, intent(in) :: ngas
188  real(dp), intent(in) :: temp
189  real(dp), intent(in) :: xg(:)
190  integer :: i
191  call modena_inputs_set(kgasinputs, kgastemppos, temp)
192  do i=1,ngas
193  call modena_inputs_set(kgasinputs, kgasxg(i), xg(i))
194  enddo
195  ret = modena_model_call(kgasmodena, kgasinputs, kgasoutputs)
196  if(ret /= 0) then
197  call exit(ret)
198  endif
199  gasmixtureconductivity=modena_outputs_get(kgasoutputs, 0_c_size_t)
200 end function gasmixtureconductivity
201 !***********************************END****************************************
202 
203 
204 !********************************BEGINNING*************************************
208 real(dp) function gasconductivity(temp,index)
209  integer, intent(in) :: index
210  real(dp), intent(in) :: temp
211  call modena_inputs_set(kginputs(index), kgtemppos(index), temp)
212  ret = modena_model_call(kgmodena(index), kginputs(index), kgoutputs(index))
213  if(ret /= 0) then
214  call exit(ret)
215  endif
216  gasconductivity=modena_outputs_get(kgoutputs(index), 0_c_size_t)
217 end function gasconductivity
218 !***********************************END****************************************
219 
220 
221 !********************************BEGINNING*************************************
225 real(dp) function solubility(temp,index)
226  integer, intent(in) :: index
227  real(dp), intent(in) :: temp
228  real(dp) :: xl1,xl2
229  xl1=1.0e-3_dp
230  xl2=1-xl1
231  call modena_inputs_set(sginputs(index), sgtemppos(index), temp)
232  call modena_inputs_set(sginputs(index), sgxl1pos(index), xl1)
233  call modena_inputs_set(sginputs(index), sgxl2pos(index), xl2)
234  ret = modena_model_call(sgmodena(index), sginputs(index), sgoutputs(index))
235  if(ret /= 0) then
236  call exit(ret)
237  endif
238  solubility=modena_outputs_get(sgoutputs(index), 0_c_size_t)
239 end function solubility
240 !***********************************END****************************************
241 
242 
243 !********************************BEGINNING*************************************
247 real(dp) function diffusivity(temp,index)
248  integer, intent(in) :: index
249  real(dp), intent(in) :: temp
250  call modena_inputs_set(dginputs(index), dgtemppos(index), temp)
251  ret = modena_model_call(dgmodena(index), dginputs(index), dgoutputs(index))
252  if(ret /= 0) then
253  call exit(ret)
254  endif
255  diffusivity=modena_outputs_get(dgoutputs(index), 0_c_size_t)
256 end function diffusivity
257 !***********************************END****************************************
258 
259 
260 !********************************BEGINNING*************************************
264 real(dp) function gasdiffusivity(temp)
265  use globals, only: pressure
266  real(dp), intent(in) :: temp
267  real(dp) :: pcA, pcB, pcApcB, TcA, TcB, TcATcB
268  real(dp) :: MA, MB, Mterm ,a, b, aToverTcsb
269  pca = 33.5e0_dp ! N2
270  pcb = 72.9e0_dp ! CO2
271  pcapcb = (pca*pcb)**(1.0e0_dp/3.0e0_dp) ! CO2, N2, B-1 p. 744
272  tca = 126.2e0_dp ! N2
273  tcb = 304.2e0_dp ! CO2
274  tcatcb = (tca*tcb)**(5.0e0_dp/12.0e0_dp)
275  ma = 28.02e0_dp
276  mb = 44.01e0_dp
277  mterm = dsqrt(1/ma + 1/mb)
278  a = 2.7450e-4_dp ! non-polar pairs
279  b = 1.823e0_dp
280  atovertcsb = a*(temp/dsqrt(tca*tcb))**b
281  ! pressure in atmospheres, cm2/s
282  gasdiffusivity = (atovertcsb*pcapcb*tcatcb*mterm)*1.0e5_dp/pressure
283  gasdiffusivity = gasdiffusivity * 1.0e-4_dp ! m2/s
284 end function gasdiffusivity
285 !***********************************END****************************************
286 
287 
288 !********************************BEGINNING*************************************
292 real(dp) function cdheatcapacity(temp)
293  real(dp), intent(in) :: temp
294  real(dp) :: &
295  t,&
296  a=24.99735_dp,&
297  b=55.18696_dp,&
298  c=-33.69137_dp,&
299  d=7.948387_dp,&
300  e=-0.136638_dp
301  t=temp/1e3
302  cdheatcapacity=a+b*t+c*t**2+d*t**3+e/t**2
303 end function cdheatcapacity
304 !***********************************END****************************************
305 
306 
307 !********************************BEGINNING*************************************
311 real(dp) function cypheatcapacity(temp)
312  real(dp), intent(in) :: temp
313  real(dp) :: &
314  t,&
315  a=-25.6132057_dp,&
316  b=226.4176882_dp,&
317  c=574.2688767_dp,&
318  d=-670.5517907_dp,&
319  e=0.6765321_dp
320  t=temp/1e3
321  cypheatcapacity=a+b*t+c*t**2+d*t**3+e/t**2
322 end function cypheatcapacity
323 !***********************************END****************************************
324 
325 
326 !********************************BEGINNING*************************************
330 real(dp) function airheatcapacity(temp)
331  real(dp), intent(in) :: temp
332  airheatcapacity=0.21_dp*oxyheatcapacity(temp)+0.79_dp*nitrheatcapacity(temp)
333 end function airheatcapacity
334 !***********************************END****************************************
335 
336 
337 !********************************BEGINNING*************************************
341 real(dp) function nitrheatcapacity(temp)
342  real(dp), intent(in) :: temp
343  real(dp) :: &
344  t,&
345  a=28.98641_dp,&
346  b=1.853978_dp,&
347  c=-9.647459_dp,&
348  d=16.63537_dp,&
349  e=0.000117_dp
350  t=temp/1e3
351  nitrheatcapacity=a+b*t+c*t**2+d*t**3+e/t**2
352 end function nitrheatcapacity
353 !***********************************END****************************************
354 
355 
356 !********************************BEGINNING*************************************
360 real(dp) function oxyheatcapacity(temp)
361  real(dp), intent(in) :: temp
362  real(dp) :: &
363  t,&
364  a=31.32234_dp,&
365  b=-20.23531_dp,&
366  c=57.86644_dp,&
367  d=-36.50624_dp,&
368  e=-0.007374_dp
369  t=temp/1e3
370  oxyheatcapacity=a+b*t+c*t**2+d*t**3+e/t**2
371 end function oxyheatcapacity
372 !***********************************END****************************************
373 
374 
375 !********************************BEGINNING*************************************
379 subroutine strutcontent(strut_content,foam_density)
380  real(dp), intent(out) :: strut_content
381  real(dp), intent(in) :: foam_density
382  !modena variables
383  integer(c_size_t) :: fspos
384  integer(c_size_t) :: rhopos
385  type(c_ptr) :: fsModena = c_null_ptr
386  type(c_ptr) :: fsInputs = c_null_ptr
387  type(c_ptr) :: fsOutputs = c_null_ptr
388  fsmodena = modena_model_new(c_char_"strutContent"//c_null_char)
389  if (modena_error_occurred()) then
390  call exit(modena_error())
391  endif
392  fsinputs = modena_inputs_new(fsmodena)
393  fsoutputs = modena_outputs_new(fsmodena)
394  rhopos = modena_model_inputs_argpos(&
395  fsmodena, c_char_"rho"//c_null_char)
396  call modena_model_argpos_check(fsmodena)
397  call modena_inputs_set(fsinputs, rhopos, foam_density)
398  ret = modena_model_call(fsmodena, fsinputs, fsoutputs)
399  if(ret /= 0) then
400  call exit(ret)
401  endif
402  strut_content=modena_outputs_get(fsoutputs, 0_c_size_t)
403  call modena_inputs_destroy (fsinputs)
404  call modena_outputs_destroy (fsoutputs)
405  call modena_model_destroy (fsmodena)
406 end subroutine strutcontent
407 !***********************************END****************************************
408 end module physicalproperties
globals::solmodel
integer, dimension(:), allocatable solmodel
solubility model
Definition: globals.f90:16
globals::diffmodel
integer, dimension(:), allocatable diffmodel
diffusivity model
Definition: globals.f90:17
globals::pressure
real(dp), dimension(:), allocatable pressure
partial pressure(Pa)
Definition: globals.f90:106
globals
namespace with global variables
Definition: globals.f90:8