MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
BASources.H
Go to the documentation of this file.
1 
7 if (PBESwitch)
8 {
9  scalar prsr, tmpt, alpha_f, g1_BA_val,
10  cc1_value, tmptrOld, timeStep, dTdt,
11  dXWdt, g1_CO2_val, XOH_val, rhoPolyS;
12  scalar henry_Coefficient, mZero_val, mOne_val,
13  bubble_radius, partialPressure_CO2, wCO2_g_val,
14  wBA_g_val, wCO2_Max,wBL_Max, partialPressure_BA;
15 
16  // call surrogate model for density reaction mixture
17  size_t T_pos =
18  modena_model_inputs_argPos(density_reaction_mixturemodel, "T");
19  size_t XOH_pos =
20  modena_model_inputs_argPos(density_reaction_mixturemodel, "XOH");
21  modena_model_argPos_check(density_reaction_mixturemodel);
22 
23  forAll(mesh.C(), celli)
24  {
25  tmpt = TS[celli];
26  tmptrOld = TS.oldTime()[celli];
27  // tmptrOldOld = TS.oldTime().oldTime()[celli];
28  timeStep = runTime.deltaTValue();
29  prsr = p[celli];
30  // rho_f = rho_foam[celli];
31  alpha_f = alpha2[celli];
32  mZero_val = mZero[celli];
33  mOne_val = mOne[celli];
34  wCO2_g_val = wCO2_g[celli];
35  wBA_g_val = wBA_g[celli];
36  cc1_value = cc1[celli];
37  g1_BA_val = g1_BA[celli];
38  XOH_val = XOH[celli];
39 
40  if (liquidMixtureDensitySurrogate)
41  {
42  // set input vector
43  modena_inputs_set(inputs_den, T_pos, tmpt);
44  modena_inputs_set(inputs_den, XOH_pos, XOH_val);
45  // call the model
46  int ret_den =
47  modena_model_call
48  (
49  density_reaction_mixturemodel,
50  inputs_den,
51  outputs_den
52  );
53  if (ret_den != 0)
54  {
55  modena_inputs_destroy (inputs_den);
56  modena_outputs_destroy (outputs_den);
57  modena_model_destroy (density_reaction_mixturemodel);
58  exit(ret_den);
59  }
60  rhoPolyS = modena_outputs_get(outputs_den, 0);
61  }
62  else
63  {
64  rhoPolyS = rhoPoly;
65  }
66 
67  bubble_radius = bubbleRadius(mZero_val, mOne_val);
68 
69  partialPressure_BA = partialPressureBA
70  (
71  M_B, M_CO2, surfaceTension, wBA_g_val,
72  wCO2_g_val, prsr, bubble_radius
73  );
74  if (blowingAgent == "n-pentane")
75  {
76  if (alpha_f > 0.5)
77  {
78  // if (LliqMax(tmpt) < L0)
79  // {
80  // dTdt = dTdtFirstOrder(tmpt, tmptrOld, timeStep);
81  // wBA_gSource[celli] = (-cc1_value*ddT_LliqMax(tmpt)*dTdt);
82 
83  wBA_lSource[celli] =
84  (
85  (M_B/1000.0)*(1.0/rhoPolyS)
86  *(partialPressure_BA/(max((RR*tmpt),ROOTVSMALL)))*(-g1_BA_val)
87  );
88  wBA_gSource[celli] = -wBA_lSource[celli];
89  // }
90  // else
91  // {
92  // wBA_gSource[celli] = 0.0;
93  // wBA_lSource[celli] = 0.0;
94  // }
95  }
96  else
97  {
98  wBA_lSource[celli] = 0.0;
99  wBA_gSource[celli] = 0.0;
100  }
101  }
102 
103  if (blowingAgent == "R-11")
104  {
105  scalar xBL_val;
106  xBL_val = xBL(tmpt, dxdT);
107  wBL_Max = wBL_D(xBL_val, M_B, M_NCO, L0);
108  if (wBL_Max > 1.0e-4 && alpha_f > 0.5)
109  {
110  dTdt = dTdtFirstOrder(tmpt, tmptrOld, timeStep);
111  wBA_gSource[celli] =
112  (
113  cc1_value*(-(M_B/M_NCO)
114  *(1.0/(Foam::pow((1.0 - xBL_val),2)))*(dxdT))*dTdt
115  );
116  wBA_lSource[celli] =
117  (
118  -cc1_value*(-(M_B/M_NCO)
119  *(1.0/(Foam::pow((1.0 - xBL_val),2)))*(dxdT))*dTdt
120  );
121  }
122  else
123  {
124  wBA_gSource[celli] = 0.0;
125  wBA_lSource[celli] = 0.0;
126  }
127  }
128  if (alpha_f > 0.5)
129  {
130  g1_CO2_val = g1_CO2[celli];
131  dXWdt =
132  (
133  (mag(XW[celli] - XW.oldTime()[celli]))
134  /(max(runTime.deltaTValue(),ROOTVSMALL))
135  );
136  partialPressure_CO2 =
137  partialPressureCO2
138  (
139  M_CO2, M_B, surfaceTension,
140  wCO2_g_val, wBA_g_val, prsr, bubble_radius
141  );
142  if (CW_0 == 0.0)
143  {
144  wCO2_lSource[celli] = 0.0;
145  }
146  else
147  {
148  if (XOH[celli] < XOH_Gel)
149  {
150  wCO2_lSource[celli] =
151  (
152  (CW_0*dXWdt*(M_CO2/1000.0)*(1.0/rhoPolyS)
153  - g1_CO2_val*(partialPressure_CO2/(RR*tmpt))*(M_CO2/1000.0)*(1.0/rhoPolyS))
154  );
155  }
156  }
157 
158  henry_Coefficient = henryCoefficient(tmpt);
159 
160  wCO2_Max = wCO2Max(M_CO2, M_liq, partialPressure_CO2, henry_Coefficient);
161 
162  if (wCO2_l[celli] < wCO2_Max || CW_0 == 0.0 || XOH[celli] > XOH_Gel)
163  {
164  wCO2_gSource[celli] = 0.0;
165 
166  }
167  else
168  {
169  wCO2_gSource[celli] =
170  (
171  g1_CO2_val
172  *(partialPressure_CO2/(RR*tmpt))*(M_CO2/1000.0)*(1.0/rhoPolyS)
173  );
174  }
175  }
176  else
177  {
178  wCO2_gSource[celli] = 0.0;
179  wCO2_lSource[celli] = 0.0;
180  }
181  }
182 }
183 else
184 {
185  scalar prsr, tmpt, alpha_f,
186  QBA, tmptrOld, timeStep, dTdt;
187 
188  scalar henry_Coefficient,
189  bubble_radius, partialPressure_CO2, wCO2_g_val,
190  wBA_g_val, wCO2_Max, wBL_Max, XOH_val, rhoPolyS;
191 
192  scalar m0_val, m1_val;
193  m0_val = m1_val = 0.0;
194  scalar xBL_val;
195 
196  // call surrogate model for density reaction mixture
197  size_t T_pos =
198  modena_model_inputs_argPos(density_reaction_mixturemodel, "T");
199  size_t XOH_pos =
200  modena_model_inputs_argPos(density_reaction_mixturemodel, "XOH");
201  modena_model_argPos_check(density_reaction_mixturemodel);
202 
203  forAll(mesh.C(), celli)
204  {
205  tmpt = TS[celli];
206  tmptrOld = TS.oldTime()[celli];
207  // tmptOldOld = TS.oldTime().oldTime()[celli];
208  timeStep = runTime.deltaTValue();
209  prsr = p[celli];
210  // rho_f = rho_foam[celli];
211  alpha_f = alpha2[celli];
212  wCO2_g_val = wCO2_g[celli];
213  wBA_g_val = wBA_g[celli];
214  // wBA_l_val = wBA_l[celli];
215  XOH_val = XOH[celli];
216 
217  if (liquidMixtureDensitySurrogate)
218  {
219  // set input vector
220  modena_inputs_set(inputs_den, T_pos, tmpt);
221  modena_inputs_set(inputs_den, XOH_pos, XOH_val);
222  // call the model
223  int ret_den =
224  modena_model_call
225  (
226  density_reaction_mixturemodel,
227  inputs_den,
228  outputs_den
229  );
230  if (ret_den != 0)
231  {
232  modena_inputs_destroy (inputs_den);
233  modena_outputs_destroy (outputs_den);
234  modena_model_destroy (density_reaction_mixturemodel);
235  exit(ret_den);
236  }
237  rhoPolyS = modena_outputs_get(outputs_den, 0);
238  }
239  else
240  {
241  rhoPolyS = rhoPoly;
242  }
243 
244  if (blowingAgent == "n-pentane")
245  {
246  if (alpha_f > 0.5)
247  {
248  dTdt = dTdtFirstOrder(tmpt, tmptrOld, timeStep);
249 
250  if (L0 < LliqMax(tmpt))
251  {
252  QBA = -ddT_LliqMax(tmpt)*dTdt;
253  wBA_gSource[celli] = QBA;
254  wBA_lSource[celli] = ddT_LliqMax(tmpt)*dTdt;
255  }
256  else
257  {
258  wBA_gSource[celli] = 0.0;
259  wBA_lSource[celli] = 0.0;
260  }
261  }
262  else
263  {
264  wBA_lSource[celli] = 0.0;
265  wBA_gSource[celli] = 0.0;
266  }
267  }
268 
269  if (blowingAgent == "R-11")
270  {
271  xBL_val = xBL(tmpt, dxdT);
272  wBL_Max = wBL_D(xBL_val, M_B, M_NCO, L0);
273  if (wBL_Max > 1.0e-4 && alpha_f > 0.5)
274  {
275  dTdt = dTdtFirstOrder(tmpt, tmptrOld, timeStep);
276  wBA_gSource[celli] =
277  (
278  (-(M_B/M_NCO)
279  *(1.0/(Foam::pow((1.0 - xBL_val),2)))*(dxdT))*dTdt
280  );
281  wBA_lSource[celli] =
282  (
283  ((M_B/M_NCO)
284  *(1.0/(Foam::pow((1.0 - xBL_val),2)))*(dxdT))*dTdt
285  );
286  }
287  else
288  {
289  wBA_lSource[celli] = 0.0;
290  wBA_gSource[celli] = 0.0;
291  }
292  }
293 
294  wCO2_lSource[celli] = scalar(0.0);
295 
296  henry_Coefficient = henryCoefficient(tmpt);
297  bubble_radius = bubbleRadius(m0_val, m1_val);
298  partialPressure_CO2 =
299  partialPressureCO2
300  (
301  M_CO2, M_B, surfaceTension,
302  wCO2_g_val, wBA_g_val, prsr, bubble_radius
303  );
304  wCO2_Max = wCO2Max(M_CO2, M_liq, partialPressure_CO2, henry_Coefficient);
305 
306  if (alpha_f > 0.5 && XOH[celli] <= XOH_Gel && CW_0 != 0.0)
307  {
308  wCO2_gSource[celli] =
309  (
310  ((CW_0*XW[celli]*(M_CO2/1000.0)*(1.0/rhoPolyS)) - wCO2_Max)
311  );
312  }
313  else
314  {
315  wCO2_gSource[celli] = 0.0;
316  }
317 
318  if (wCO2_gSource[celli] < 0.0)
319  {
320  wCO2_gSource[celli] = scalar(0.0);
321  }
322  }
323 }
modena_model_inputs_argPos
size_t modena_model_inputs_argPos(const modena_model_t *self, const char *name)
Function determining position of an argument in the input vector.
Definition: model.c:366
henryCoefficient
double henryCoefficient(double &T)
Henry coefficient for CO2.
Definition: KineticsFunctions.H:185
partialPressureCO2
double partialPressureCO2(const state_type &y)
partial pressure of CO2
Definition: partialPressure.h:89
surfaceTension
double surfaceTension
required for the computation of partial pressure
Definition: experimentalInputs.h:136
M_CO2
double M_CO2
Molecular mass of carbon dioxide, kg/kmol.
Definition: experimentalInputs.h:129
M_NCO
double M_NCO
Molecular weight of NCO, kg/kmol.
Definition: experimentalInputs.h:131
modena_model_call
int modena_model_call(modena_model_t *self, modena_inputs_t *inputs, modena_outputs_t *outputs)
Function calling the surrogate model and checking for errors.
Definition: model.c:553
RR
double RR
ideal gas constant, J/mol K
Definition: experimentalInputs.h:108
ddT_LliqMax
double ddT_LliqMax(double &)
derivative of LliqMax with respect to temperature
Definition: KineticsFunctions.H:156
bubbleRadius
double bubbleRadius(const double m0, const double m1)
radius of bubbles based on the moments
Definition: bubbleRadius.h:28
LliqMax
double LliqMax(double &)
maximum allowable amount of liquid blowing agent (n-pentane) in liquid
Definition: KineticsFunctions.H:139
wBL_D
double wBL_D(double &xBL, double &M_B, double &M_NCO, double &L0)
weight fraction of maximum allowable blowing agent (R-11) in liquid
Definition: KineticsFunctions.H:294
rhoPoly
double rhoPoly
density of the liquid polymer, kg/m3
Definition: experimentalInputs.h:109
modena_model_destroy
void modena_model_destroy(modena_model_t *self)
Function deallocating the memory allocated for the surrogate model.
Definition: model.c:665
M_B
double M_B
Molecular mass of blowing agent, kg/kmol.
Definition: experimentalInputs.h:130
xBL
double xBL(double &T, double &dxdT)
mole fraction of blowing agent (R-11) in liquid polymer
Definition: KineticsFunctions.H:275
partialPressureBA
double partialPressureBA(const state_type &y)
partial pressure of the physical blowing agent
Definition: partialPressure.h:65
L0
double L0
Initial weight fraction of blowing agent in the liquid, -.
Definition: experimentalInputs.h:133
wCO2Max
double wCO2Max(double &M_CO2, double &M_liq, double &pCO2, double &henryCoeff)
dissolved amount of CO2 in liquid
Definition: KineticsFunctions.H:260
modena_model_argPos_check
void modena_model_argPos_check(const modena_model_t *self)
Function checking that the user has queried all input positions.
Definition: model.c:408
dTdtFirstOrder
double dTdtFirstOrder(double &Tc, double &To, double &timeStep)
first order accurate dTdt
Definition: KineticsFunctions.H:328