38 use json_module,
only : json_file
55 real(kind=
rp) :: kappa = 0.41_rp
57 real(kind=
rp) :: b = 5.2_rp
62 procedure, pass(this) :: init_from_components => &
69 procedure,
private, pass(this) ::
solve
82 type(
coef_t),
intent(in) :: coef
83 integer,
intent(in) :: msk(:)
84 integer,
intent(in) :: facet(:)
85 real(kind=
rp),
intent(in) :: nu
86 integer,
intent(in) :: h_index
87 type(json_file),
intent(inout) :: json
88 real(kind=
rp) :: kappa, b
93 call this%init_from_components(coef, msk, facet, nu, h_index, kappa, b)
107 type(
coef_t),
intent(in) :: coef
108 integer,
intent(in) :: msk(:)
109 integer,
intent(in) :: facet(:)
110 integer,
intent(in) :: h_index
111 real(kind=
rp),
intent(in) :: nu
112 real(kind=
rp),
intent(in) :: kappa
113 real(kind=
rp),
intent(in) :: b
116 call neko_error(
"Spalding's law is only available on the CPU backend.")
119 call this%init_base(coef, msk, facet, nu, h_index)
130 call this%free_base()
139 real(kind=
rp),
intent(in) :: t
140 integer,
intent(in) :: tstep
145 real(kind=
rp) :: ui, vi, wi, magu, utau, normu, guess
153 ui = u%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), this%ind_e(i))
154 vi = v%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), this%ind_e(i))
155 wi = w%x(this%ind_r(i), this%ind_s(i), this%ind_t(i), this%ind_e(i))
158 normu = ui * this%n_x%x(i) + vi * this%n_y%x(i) + wi * this%n_z%x(i)
160 ui = ui - normu * this%n_x%x(i)
161 vi = vi - normu * this%n_y%x(i)
162 wi = wi - normu * this%n_z%x(i)
164 magu = sqrt(ui**2 + vi**2 + wi**2)
167 if (tstep .eq. 1)
then
168 guess = sqrt(magu * this%nu / this%h%x(i))
170 guess = this%tau_x(i)**2 + this%tau_y(i)**2 + this%tau_z(i)**2
171 guess = sqrt(sqrt(guess))
174 utau = this%solve(magu, this%h%x(i), guess)
177 this%tau_x(i) = -utau**2 * ui / magu
178 this%tau_y(i) = -utau**2 * vi / magu
179 this%tau_z(i) = -utau**2 * wi / magu
188 function solve(this, u, y, guess)
result(utau)
190 real(kind=
rp),
intent(in) :: u
191 real(kind=
rp),
intent(in) :: y
192 real(kind=
rp),
intent(in) :: guess
193 real(kind=
rp) :: yp, up, kappa, b, utau
194 real(kind=
rp) :: error, f, df, old
195 integer :: niter, k, maxiter
205 yp = y * utau / this%nu
210 f = (up + exp(-kappa*b)* &
211 (exp(kappa*up) - 1.0_rp - kappa*up - 0.5_rp*(kappa*up)**2 - &
212 1.0_rp/6*(kappa*up)**3) - yp)
214 df = (-y / this%nu - u/utau**2 - kappa*up/utau*exp(-kappa*b) * &
215 (exp(kappa*up) - 1 - kappa*up - 0.5*(kappa*up)**2))
220 error = abs((old - utau)/old)
222 if (error < 1e-3)
then
229 write(*,*)
"Newton not converged", error, f, utau, old, guess
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Defines a mapping of the degrees of freedom.
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Utilities for retrieving parameters from the case files.
integer, parameter, public neko_log_debug
Debug log level.
type(log_t), public neko_log
Global log stream.
integer, parameter neko_bcknd_device
integer, parameter, public rp
Global precision used in computations.
real(kind=rp) function solve(this, u, y, guess)
Newton solver for the algebraic equation defined by the law.
subroutine spalding_compute(this, t, tstep)
Compute the wall shear stress.
subroutine spalding_init(this, coef, msk, facet, nu, h_index, json)
Constructor from JSON.
subroutine spalding_init_from_components(this, coef, msk, facet, nu, h_index, kappa, B)
Constructor from components.
subroutine spalding_free(this)
Destructor for the spalding_t (base) class.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Wall model based on Spalding's law of the wall. Reference: http://dx.doi.org/10.1115/1....
Base abstract type for wall-stress models for wall-modelled LES.