/home/runner/work/amr-wind/amr-wind/amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H Source File

AMR-Wind API: /home/runner/work/amr-wind/amr-wind/amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H Source File
AMR-Wind API v0.1.0
CFD solver for wind plant simulations
Loading...
Searching...
No Matches
stokes_waves_K.H
Go to the documentation of this file.
1#ifndef STOKES_WAVES_H_
2#define STOKES_WAVES_H_
3
4#include <AMReX_FArrayBox.H>
5#include <cmath>
6
7// Stokes waves theory adapted from
8// Fenton, J., Fifth Order Stokes Theory for Steady Waves
9// Journal of Waterway, Port, Coastal and Ocean Engineering, 1985, 111, 216-234
10
11// Updated Table 1 from Fenton 1985 found in
12// Fenton, J., Nonlinear Wave Theories
13// The Sea Vol. 9 Ocean Engineering Science, 1990
14// https://johndfenton.com/Papers/Fenton90b-Nonlinear-wave-theories.pdf
15
16// Relevant results are summarized in
17// Kinnas S. A., Notes on fifth order gravity wave theory
18// https://www.caee.utexas.edu/prof/kinnas/ce358/oenotes/kinnas_stokes11.pdf
19
20namespace amr_wind::ocean_waves::relaxation_zones {
21
22// Compute wavelength as a function of wave period, water depth, and g
23AMREX_FORCE_INLINE amrex::Real stokes_wave_length(
24 const amrex::Real T,
25 const amrex::Real d,
26 const amrex::Real H,
27 const int order,
28 const amrex::Real g,
29 const amrex::Real tol,
30 const int iter_max)
31{
32 // Calculate constants and derivatives that do not change with iteration
33 const amrex::Real omega = 2.0 * M_PI / T;
34 const amrex::Real depsdk = H / 2.0;
35
36 // Begin Newton-Raphson iterations
37 int iter = 0;
38 // First guess is first-order
39 amrex::Real k = (omega * omega) / g;
40 // Cannot skip loop
41 amrex::Real f = tol + 1.0;
42 while (std::abs(f) > tol && iter < iter_max) {
43 // Calculate current constants
44
45 // Exponential definition of S = sech(2kd)
46 const amrex::Real S =
47 2. * std::exp(2. * k * d) / (std::exp(4. * k * d) + 1.);
48 const amrex::Real C = 1.0 - S;
49 const amrex::Real eps = k * H / 2.0;
50 const amrex::Real C0 = std::sqrt(std::tanh(k * d));
51 const amrex::Real C2 = C0 * (2.0 + 7.0 * S * S) / (4.0 * C * C);
52 const amrex::Real numterm_C4 =
53 (4.0 + 32.0 * S - 116.0 * std::pow(S, 2) - 400.0 * std::pow(S, 3) -
54 71.0 * std::pow(S, 4) + 146.0 * std::pow(S, 5));
55 const amrex::Real C4 = C0 * numterm_C4 / (32.0 * std::pow(C, 5));
56 // Calculate pure derivates
57 const amrex::Real dSdk = -2.0 * d * std::sinh(2.0 * k * d) /
58 std::pow(std::cosh(2.0 * k * d), 2);
59 const amrex::Real dCdk = -dSdk;
60 const amrex::Real dC0dk =
61 d / (2.0 * C0 * std::pow(std::cosh(k * d), 2));
62 // Calculate derivatives with products
63 const amrex::Real dC2dk =
64 (4 * std::pow(C, 2) *
65 (dC0dk * (2.0 + 7.0 * std::pow(S, 2)) + C0 * 14.0 * S * dSdk) -
66 C0 * (2.0 + 7.0 * std::pow(S, 2)) * 8.0 * C * dCdk) /
67 (16.0 * std::pow(C, 4));
68 const amrex::Real dC4dk =
69 (32.0 * std::pow(C, 5) *
70 (dC0dk * numterm_C4 +
71 C0 * (32.0 * dSdk - 232.0 * S * dSdk -
72 1200 * std::pow(S, 2) * dSdk -
73 284.0 * std::pow(S, 3) * dSdk + 730 * std::pow(S, 4))) -
74 C0 * numterm_C4 * 160.0 * std::pow(C, 4) * dCdk) /
75 (1024.0 * std::pow(C, 10));
76
77 // Calculate derivative for loop convergence
78 amrex::Real dfdk = g * std::pow(C0, 2) + g * k * 2.0 * C0 * dC0dk;
79 // Add additional terms depending on order
80 if (order >= 2) {
81 dfdk +=
82 g * (2.0 * C0 * std::pow(eps, 2) * C2 +
83 std::pow(eps, 4) * std::pow(C2, 2)) +
84 g * k *
85 (2.0 * dC0dk * std::pow(eps, 2) * C2 +
86 2.0 * C0 *
87 (2.0 * eps * depsdk * C2 + std::pow(eps, 2) * dC2dk) +
88 4.0 * std::pow(eps, 3) * std::pow(C2, 2) * depsdk +
89 std::pow(eps, 4) * 2.0 * C2 * dC2dk);
90 }
91 if (order >= 4) {
92 dfdk += g * (2.0 * std::pow(eps, 4) * C0 * C4 +
93 2.0 * std::pow(eps, 6) * C2 * C4 +
94 std::pow(eps, 8) * std::pow(C4, 2)) +
95 g * k *
96 (8.0 * std::pow(eps, 3) * depsdk * C0 * C4 +
97 2.0 * std::pow(eps, 4) * (C0 * dC4dk + C4 * dC0dk) +
98 12 * std::pow(eps, 5) * depsdk * C2 * C4 +
99 2.0 * std::pow(eps, 6) * (C2 * dC4dk + C4 * dC2dk) +
100 8.0 * std::pow(eps, 7) * depsdk * std::pow(C4, 2) +
101 std::pow(eps, 8) * 2.0 * C4 * dC4dk);
102 }
103 k = k - f / dfdk;
104
105 iter += 1;
106 f = g * k * std::pow(C0, 2);
107 // Add additional terms depending on order
108 if (order >= 2) {
109 f += g * k *
110 (2.0 * C0 * std::pow(eps, 2) * C2 +
111 std::pow(eps, 4) * std::pow(C2, 2));
112 }
113 if (order >= 4) {
114 f += g * k *
115 (2.0 * std::pow(eps, 4) * C0 * C4 +
116 2.0 * std::pow(eps, 6) * C2 * C4 +
117 std::pow(eps, 8) * std::pow(C4, 2));
118 }
119 // Subtract omega^2 to measure convergence
120 f -= omega * omega;
121 }
122
123 if (k < tol) {
124 // Return negative wavelength if faulty result
125 return -1;
126 }
127 if (std::isnan(k)) {
128 return -2;
129 }
130 if (iter == iter_max) {
131 return -3;
132 }
133 // Return wavelength calculated from wavenumber
134 return 2.0 * M_PI / k;
135}
136
137AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(
138 int stokes_order,
139 amrex::Real wavenumber,
140 amrex::Real water_depth,
141 amrex::Real& c0,
142 amrex::Real& a11,
143 amrex::Real& a22,
144 amrex::Real& b22,
145 amrex::Real& c2,
146 amrex::Real& a31,
147 amrex::Real& a33,
148 amrex::Real& b31,
149 amrex::Real& a42,
150 amrex::Real& a44,
151 amrex::Real& b42,
152 amrex::Real& b44,
153 amrex::Real& c4,
154 amrex::Real& a51,
155 amrex::Real& a53,
156 amrex::Real& a55,
157 amrex::Real& b53,
158 amrex::Real& b55)
159{
160
161 amrex::Real kd = wavenumber * water_depth;
162 if (kd > 50. * M_PI) {
163 kd = 50 * M_PI;
164 }
165
166 // Exponential definition of S = sech(2kd)
167 amrex::Real S = 2. * std::exp(2. * kd) / (std::exp(4. * kd) + 1.);
168 amrex::Real C = 1.0 - S;
169 amrex::Real Sh = std::sinh(kd);
170 amrex::Real Th = std::tanh(kd);
171 // Exponential definition of coth(kd)
172 amrex::Real CTh = (1. + std::exp(-2 * kd)) / (1. - std::exp(-2 * kd));
173
174 c0 = std::sqrt(Th);
175 a11 = 1. / std::sinh(kd); // Hyperbolic cosecant
176 // Second order coefficients
177 a22 = 3. * std::pow(S, 2) / (2 * std::pow(C, 2));
178 b22 = CTh * (1 + 2. * S) / (2 * C);
179 c2 = std::sqrt(Th) * (2 + 7 * std::pow(S, 2)) / (4 * std::pow(C, 2));
180 if (stokes_order == 2) {
181 return;
182 }
183
184 // Third order coefficients
185 a31 = (-4 - 20 * S + 10 * std::pow(S, 2) - 13 * std::pow(S, 3)) /
186 (8 * Sh * std::pow(C, 3));
187 a33 =
188 (-2 * std::pow(S, 2) + 11 * std::pow(S, 3)) / (8 * Sh * std::pow(C, 3));
189 b31 = -3 * (1 + 3 * S + 3 * std::pow(S, 2) + 2 * std::pow(S, 3)) /
190 (8 * std::pow(C, 3));
191 if (stokes_order == 3) {
192 return;
193 }
194
195 // Fourth order coefficients
196 a42 = (12 * S - 14 * std::pow(S, 2) - 264 * std::pow(S, 3) -
197 45 * std::pow(S, 4) - 13 * std::pow(S, 5)) /
198 (24 * std::pow(C, 5));
199 a44 = (10 * std::pow(S, 3) - 174 * std::pow(S, 4) + 291 * std::pow(S, 5) +
200 278 * std::pow(S, 6)) /
201 (48 * (3 + 2 * S) * std::pow(C, 5));
202 b42 = CTh *
203 (6 - 26 * S - 182 * std::pow(S, 2) - 204 * std::pow(S, 3) -
204 25 * std::pow(S, 4) + 26 * std::pow(S, 5)) /
205 (6 * (3 + 2 * S) * std::pow(C, 4));
206 b44 = CTh *
207 (24 + 92 * S + 122 * std::pow(S, 2) + 66 * std::pow(S, 3) +
208 67 * std::pow(S, 4) + 34 * std::pow(S, 5)) /
209 (24 * (3 + 2 * S) * std::pow(C, 4));
210 c4 = std::sqrt(Th) *
211 (4 + 32 * S - 116 * std::pow(S, 2) - 400 * std::pow(S, 3) -
212 71 * std::pow(S, 4) + 146 * std::pow(S, 5)) /
213 (32 * std::pow(C, 5));
214 if (stokes_order == 4) {
215 return;
216 }
217
218 // Fifth order coefficients
219 a51 =
220 (-1184 + 32 * S + 13232 * std::pow(S, 2) + 21712 * std::pow(S, 3) +
221 20940 * std::pow(S, 4) + 12554 * std::pow(S, 5) -
222 500 * std::pow(S, 6) - 3341 * std::pow(S, 7) - 670 * std::pow(S, 8)) /
223 (64 * Sh * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
224 a53 = (4 * S + 105 * pow(S, 2) + 198 * std::pow(S, 3) -
225 1376 * std::pow(S, 4) - 1302 * std::pow(S, 5) -
226 117 * std::pow(S, 6) + 58 * std::pow(S, 7)) /
227 (32 * Sh * (3 + 2 * S) * std::pow(C, 6));
228 a55 =
229 (-6 * std::pow(S, 3) + 272 * std::pow(S, 4) - 1552 * std::pow(S, 5) +
230 852 * std::pow(S, 6) + 2029 * std::pow(S, 7) + 430 * std::pow(S, 8)) /
231 (64 * Sh * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
232 b53 = 9 *
233 (132 + 17 * S - 2216 * std::pow(S, 2) - 5897 * std::pow(S, 3) -
234 6292 * std::pow(S, 4) - 2687 * std::pow(S, 5) +
235 194 * std::pow(S, 6) + 467 * std::pow(S, 7) + 82 * std::pow(S, 8)) /
236 (128 * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
237 b55 =
238 5 *
239 (300 + 1579 * S + 3176 * std::pow(S, 2) + 2949 * std::pow(S, 3) +
240 1188 * std::pow(S, 4) + 675 * std::pow(S, 5) + 1326 * std::pow(S, 6) +
241 827 * std::pow(S, 7) + 130 * std::pow(S, 8)) /
242 (384 * (3 + 2 * S) * (4 + S) * std::pow(C, 6));
243 if (stokes_order == 5) {
244 return;
245 }
246
247 if (stokes_order > 5 || stokes_order < 2) {
248 amrex::Abort(
249 "invalid stokes order specified. It should be between 2,3,4 or 5 ");
250 }
251}
252
253// Based on Fenton 1985
254AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_waves(
255 int stokes_order,
256 amrex::Real wavelength,
257 amrex::Real water_depth,
258 amrex::Real wave_height,
259 amrex::Real zsl,
260 amrex::Real g,
261 amrex::Real x,
262 amrex::Real z,
263 amrex::Real time,
264 amrex::Real phase_offset,
265 amrex::Real& eta,
266 amrex::Real& u_w,
267 amrex::Real& v_w,
268 amrex::Real& w_w)
269{
270 const amrex::Real wavenumber = 2.0 * M_PI / wavelength;
271
272 // some parameters
273 amrex::Real c0{0.0};
274 amrex::Real a11{0.0}, a22{0.0}, b22{0.0}, c2{0.0};
275 amrex::Real a31{0.0}, a33{0.0}, b31{0.0}, a42{0.0}, a44{0.0};
276 amrex::Real b42{0.0}, b44{0.0}, c4{0.0};
277 amrex::Real a51{0.0}, a53{0.0}, a55{0.0}, b53{0.0}, b55{0.0};
278
279 stokes_coefficients(
280 stokes_order, wavenumber, water_depth, c0, a11, a22, b22, c2, a31, a33,
281 b31, a42, a44, b42, b44, c4, a51, a53, a55, b53, b55);
282
283 const amrex::Real eps = wavenumber * wave_height / 2.; // Steepness (ka)
284 const amrex::Real c = (c0 + std::pow(eps, 2) * c2 + std::pow(eps, 4) * c4) *
285 std::sqrt(g / wavenumber);
286
287 const amrex::Real omega = c * wavenumber;
288 const amrex::Real phase = wavenumber * x - omega * time - phase_offset;
289
290 eta = (eps * std::cos(phase) // first order term
291 + std::pow(eps, 2) * b22 * std::cos(2. * phase) // second order term
292 + std::pow(eps, 3) * b31 * (std::cos(phase) - std::cos(3. * phase)) +
293 std::pow(eps, 4) *
294 (b42 * std::cos(2. * phase) + b44 * std::cos(4. * phase)) +
295 std::pow(eps, 5) *
296 (-(b53 + b55) * std::cos(phase) + b53 * std::cos(3 * phase) +
297 b55 * std::cos(5 * phase))) /
298 wavenumber +
299 zsl;
300
301 // Compute velocities components using Eq.(21) Eq.(23) from Kinnas
302 // https://www.caee.utexas.edu/prof/kinnas/ce358/oenotes/kinnas_stokes11.pdf
303 // Define coefficients using Eq.(19)
304
305 const int MAX_ORDER = 5;
306 amrex::GpuArray<amrex::Real, MAX_ORDER> a;
307 if (stokes_order == 2) {
308 a[0] = a11;
309 a[1] = a22;
310 }
311 if (stokes_order == 3) {
312 a[0] = a11 + (eps * eps) * a31;
313 a[1] = a22;
314 a[2] = a33;
315 }
316 if (stokes_order == 4) {
317 a[0] = a11 + (eps * eps) * a31;
318 a[1] = a22 + (eps * eps) * a42;
319 a[2] = a33;
320 a[3] = a44;
321 }
322 if (stokes_order == 5) {
323 a[0] = a11 + (eps * eps) * a31 + std::pow(eps, 4) * a51;
324 a[1] = a22 + (eps * eps) * a42;
325 a[2] = a33 + (eps * eps) * a53;
326 a[3] = a44;
327 a[4] = a55;
328 }
329
330 u_w = 0.;
331 v_w = 0.;
332 w_w = 0.;
333 for (int n = 1; n <= stokes_order; ++n) {
334 // Upper bound for deep water case
335 // This ensure finite values of velocity for large kd's
336 amrex::Real nkdz = n * wavenumber * (water_depth + (z - zsl));
337 if (nkdz > 50. * M_PI) {
338 nkdz = 50. * M_PI;
339 }
340 u_w += std::pow(eps, n) * n * a[n - 1] * std::cosh(nkdz) *
341 std::cos(n * phase);
342 w_w += std::pow(eps, n) * n * a[n - 1] * std::sinh(nkdz) *
343 std::sin(n * phase);
344 }
345 u_w *= (c0 * std::sqrt(g / wavenumber));
346 w_w *= (c0 * std::sqrt(g / wavenumber));
347}
348
349} // namespace amr_wind::ocean_waves::relaxation_zones
350#endif
amr_wind::ocean_waves::relaxation_zones
Definition relaxation_zones_ops.cpp:13
amr_wind::ocean_waves::relaxation_zones::stokes_coefficients
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_coefficients(int stokes_order, amrex::Real wavenumber, amrex::Real water_depth, amrex::Real &c0, amrex::Real &a11, amrex::Real &a22, amrex::Real &b22, amrex::Real &c2, amrex::Real &a31, amrex::Real &a33, amrex::Real &b31, amrex::Real &a42, amrex::Real &a44, amrex::Real &b42, amrex::Real &b44, amrex::Real &c4, amrex::Real &a51, amrex::Real &a53, amrex::Real &a55, amrex::Real &b53, amrex::Real &b55)
Definition stokes_waves_K.H:137
amr_wind::ocean_waves::relaxation_zones::stokes_waves
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void stokes_waves(int stokes_order, amrex::Real wavelength, amrex::Real water_depth, amrex::Real wave_height, amrex::Real zsl, amrex::Real g, amrex::Real x, amrex::Real z, amrex::Real time, amrex::Real phase_offset, amrex::Real &eta, amrex::Real &u_w, amrex::Real &v_w, amrex::Real &w_w)
Definition stokes_waves_K.H:254
amr_wind::ocean_waves::relaxation_zones::stokes_wave_length
AMREX_FORCE_INLINE amrex::Real stokes_wave_length(const amrex::Real T, const amrex::Real d, const amrex::Real H, const int order, const amrex::Real g, const amrex::Real tol, const int iter_max)
Definition stokes_waves_K.H:23
Generated by doxygen 1.12.0