Closed anshchaube closed 4 years ago
@dshaver-ANL @katyhuff this is what I am doing for Boussinesq in Nek's userf
function (explanation follows):
g = 9.8
dh = 0.015845536
uav = 0.219457
uavsq = uav*uav
rhoe = 2692.0 !Janz 1988
a = 0.562
alpha = a/(rhoe-(a*temp))
tref = 0.0
ffx = 0.0
ffy = 0.0
ffz = (g*dh/uavsq)*(1-alpha*(temp-tref))
The Boussinesq term ffz
doesn't have a rho value multiplied to it as Nek automatically multiplies ffz by rho (which is 1 for our non-dimensionalized case anyway).
For the non dimensionalized case, the Boussinesq term is:
ffz = (gL / Vav^2) [ 1 - alpha_exp (T - Tref)
where L is the characteristic length (hydraulic dia), Vav the characteristic velocity, and the coefficient of expansion alpha_exp
is
alpha = - (1/rho) d(rho)/dT
rho_{molten-salt} = rho_experimental - (constant)*T
Plugging this rho into the definition of alpha yields what I have coded up.
Does this look correct? Furthermore, is this appropriate for our case where the difference between the inlet and the outlet temperatures is about 30 K ? Just for reference, our dRho/Rho_ref over this range is ~ 0.7%.
This is blowing up, just like #10. This one has lx1=8, genmap interval = 0.1, and variable dT with target CFL of 0.2.
Step 422, t= 1.0117773E+00, DT= 3.2136566E-42, C=******* 3.9314E+03 9.3883E+00
Solving for Hmholtz scalars
422 Scalars done 1.0118E+00 4.9035E-02
Solving for fluid
422 Project PRES NaN NaN NaN 7 8
422 PRES gmres 200 NaN NaN 1.0000E-05 4.6727E+00 8.9877E+00 F
L1/L2 DIV(V) 3.4499+120 6.2886+125
L1/L2 QTL 0.0000E+00 0.0000E+00
L1/L2 DIV(V)-QTL 3.4499+120 6.2886+125
WARNING: DIV(V)-QTL too large!
422 Fluid done 1.0118E+00 9.2433E+00
Step 423, t= 1.0117773E+00, DT= 2.5709253E-42, C=******* 3.9409E+03 9.4704E+00
This issue can be closed when Boussinesq approximation has been implemented and tested in our flat plenum model.