EXPERIMENTS WITH Ri-DEPENDENT SGS FORMULATIONS: (Run in: sultan:~/nek5c/2d07) ********************************************************************************* * smari41a: Ra=312500, dt=4.e-6 .. cs=0.05; Ri_c=1.0 Let's try to rejuvenate the numerically easy case via Ri-dependent formulation (if it works, of course)? - factor_Ri activated for u only: crashed and burned at 3 steps. - factor_Ri activated for T only: crashed and burned at 3 steps; look for bugs. - Parameter Ra is specified as Ra=5.e6 in smag routine (it is not passed from userf); change manually: crashed and burned at 3 steps. Look for other bugs: Calculation of Ri seems OK. Calculation of factor seems OK. Change if statement to IF((du/dz^2+dv/dy^2).gt.0.) then... since they are typically less than 1.e-6. - Flip to abs(Ri): treats Ri<0, unstable overturning, same as normal stuff, but keeps factor<=1 (for homog Smag, factor=1). Otherwise factor can exceed 1, presumably trying to model mixing due to convective instability. factor_Ri still activated for T only: it runs fine, but not a whole lot difference from smari40b. - Activate factor_Ri for both u and T: runs but not too different. Maybe Ri_c=1 is too inclusive for this smooth run, things behaving much like homog. Smag (even at Ri=0.75, factor=0.5). - Crank down Ri_c=0.25, theoretical number (differs from lab.): very slight difference from homog. Smag. Some problem at the very bottom as expected (stable place, r=0). ********************************************************************************* * smari41b: Ra=5.e6, dt=1.e-6 .. cs=0.05; Ri_c=0.25 There is naturally more structure in this case, giving more opportunity for Ri-dependent formulation to do some work. However, low Ri_c may create problems... It finished, but only slightly different from smari40a. It seems that no more information is provided in addition to what strain tensor can supply. CONCLUSION: Since deformation tensor in Smag. automatically takes care of anisotropy (actually better), efforts to cut off vertical mixing coefficient via Ri dependence is not going to make a big difference. It should apply to the entire SGS formulation. ********************************************************************************* * smari42a: Ra=5.e6, dt=1.e-6 .. cs=0.05; Ri_c=0.25 Define turbulent Prandtl number as a function of Ri in order to use more info. In u, factor=1 : leave it alone for the time being. Any impact could also come from forcing, as the status of temperature is immediately felt by the velocity field.Flow field Will simply react via homog. Smag. sgs. In T, if(Ri.lt.Ri_c) then factor = 1./(Ri/Ri_c*(Pr_lam-1.)+1.) else factor_Ri = 1./Pr_lam endif which is a function that changes the value of Pr from 1 at Ri=0 to Pr_lam at Ri=Ri_c and higher, where Pr is "laminar" value (10 for T, 1000 for S). I.e., if there is no turbulent overturning, entire eddy diffusivity applied to temperature fields will be reduced so that interfaces will be kept sharp. This is supported by various lab data, listed in Fig. 6 of Rohr et al. (JFM, 195, p. 88). We can worry about the shape of the curve later. Let's start with Pr_lam=10. step (x1000) | 1 | 2 | 3 | 4 | 5 | 6 | 7 ---------------------------------------------------------------------- Tmax | 1.027 | 1.242 | 1.349 | 1.420 | 1.385 | 1.235 | 1.138 Sharper features, but overshoots at the bottom, expected, stable stratification, low eddy viscosity. ********************************************************************************* * smari42b: Ra=80.e6, dt=0.25e-6; cs=0.05; Ri_c=0.25l Pr_lam=10 Will it work for higher Ra? It works. Generally sharper features, but not a drastic improvement over smari40c. ********************************************************************************* * smari43a: Let's take it a step further: apply Ri-dependence to all sgs terms both in u and T: if(Ri.lt.Ri_c) then factor = sqrt(1-Ri/Ri_c) ! use sqrt function, but not important else factor = 0.0 endif do j=1,ndim c if (j.eq.2) then c factor = factor_Ri ! all directions affected by Ri now. c else ! this is not needed anymore c factor = 1.0 c endif This implies that when dT/dz=0 (Ri=0), homog case is recovered, but sgs terms are eliminated when stratification acts to suppress mixing (Ri>Ri_c). run it with: Ra=5.e6, dt=1.e-6; cs=0.05; Ri_c=0.25, etc. looks OK. But may need to look into higher forcing to see more difference. ********************************************************************************* * smari43b: Ra=80.e6, dt=0.25e-6; cs=0.05; Ri_c=0.25 This is pretty wild (it looks "sick") and differs significantly from smari40c. step (x1000) | 1 | 2 | 3 | 4 | 5 | 6 | 7 ---------------------------------------------------------------------- Tmax | 1.055 | 1.436 | 1.645 | 1.497 | 1.507 | 1.484 | 1.036 step (x1000) | 8 | 9 | 10 | 11 | 12 ------------------------------------------------------ Tmax | 1.068 | 1.048 | 1.066 | | The overshoots are not wide spread, but are limited to only 1 location, just underneath the nose. Could be related to unstable case when Ri<0, which is badly handled here (factor=0). This is especially prominent under the nose, which overruns light water, leading to very high -dS/dz, and Ri<0. But the present model takes abs and applies factor=0. ********************************************************************************* * smari43c: Ra=80.e6, dt=0.25e-6; cs=0.05; Ri_c=0.25 Scheme to allow max. eddy diffusion (factor=1) in unstable regions (Ri<0): Ri = -Ra*w_T(kx,ky,kz,1,2)... ! remove abs(), correct sign mistake ... if(Ri.ge.Ri_c) then factor = 0.0 elseif(Ri.le.0.) then factor = 1.0 else factor = sqrt(1-Ri/Ri_c) endif step (x1000) | 1 | 2 | 3 | 4 | 5 | 6 | 7 ---------------------------------------------------------------------- Tmax | 1.020 | 1.084 | 1.044 | 1.122 | 1.056 | 1.067 | 1.016 step (x1000) | 8 | 9 | 10 | 11 | 12 ------------------------------------------------------ Tmax | 1.037 | 1.066 | 1.053 | 1.059 | 1.035 Compare to: smari40c (homog. Smag.) and smari43b (previous version with Ri): The results are quite close smari40c (homog. Smag.) but with sharper features. (and hence very different from smari43b). But continuous action leads to different evolutions by the time it reaches the bottom (nonlinearity). It looks quite good to me because the changes are 2nd order to the main Smag. model, but seem to be in the right direction. Some unharmful overshoots remain, but I cannot pinpoint the source/reason. Could be internal model response to zero SGS diffusivity here and there and some "background" value in the conventional Laplacian via vert_diff_ratio (now zero, maybe something like 1.e-4) may help. Behavior needs to be tested for (a) higher Ra, (b)higher resolution and (c) in 3d to see if it creates serious problems. * vert_diff_ratio = 1.e-4 in subroutine reset_vert_diff Add small background diffusion to help reduce numerical problems for regions with Ri>Ri_c. step (x1000) | 1 | 2 | 3 | 4 | 5 | 6 | 7 ---------------------------------------------------------------------- Tmax | 1.020 | 1.084 | 1.044 | 1.122 | 1.052 | 1.067 | 1.014 No difference, but let's keep it there. * Try c_mild=500 step (x1000) | 1 | 2 | 3 | 4 | 5 | 6 -------------------------------------------------------------- Tmax | 1.008 | 1.034 | 1.026 | 1.086 | 1.022 | 1.067 Sligthly better overshoots (nothing visible), leave it at that. ********************************************************************************* * smari43d: Ra=320.e6, dt=0.125e-6; cs=0.05; Ri_c=0.25; Will it differ from the case with Ra=80.e6 (homog. Smag did not)? It is different from the case with Ra=80.e6, but not much more scales. In any case, it is not identical as in the case of homog. Smag. model. ********************************************************************************* * smari44a: clean up a little bit, Ra has to be specified 3 times in userf, compute_sma and compute_sma_T, and Ri_c has to be specified twice. Remove Ra from compute_sma and compute_sma_T (userf only) and move Ri_c to userf_set, where sma_constant and other sgs variables are set. Clean up some other stuff. OK, seems to give identical results for Ra=80.e6. Reduce to Ra=5.e6