# sandbox/abockmann/Stokes5.h

 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 ``` ``````#ifndef Stokes5_H #define Stokes5_H #define eta Stokes5_eta // rename functions #define u Stokes5_u #define v Stokes5_v // #include // #include #ifndef PI #define PI 3.14159265358979323846 #endif typedef struct{ /* Coefficients depending on wave and environmental properties */ double A11, A22, A31, A33, A42, A44, A51, A53, A55, B22, B31, B42, B44, B53, B55, C0, C2, C4, D2, D4, E2, E4, S, eps, wave_length, wave_height, current, gravity, depth, k, kd, kH, c, cs, u_mean, Q, R, x0, y0; } Stokes5; /* Stokes 5th order waves based on Fenton's paper */ void set_stokes5_properties(Stokes5 *stokes5, double wave_length, double wave_height, double current, double depth, double gravity, double x0, double y0) /* Compute coefficients based on wave properties and depth */ { double eps, S, k, kd; // compute basic properties, dimensional numbers, and set attributes of stokes5 k = 2.*PI/wave_length; kd = k*depth; eps = k*wave_height/2.; S = 1./cosh(2.*kd); stokes5->k = k; stokes5->kd = kd; stokes5->wave_height = wave_height; stokes5->wave_length = wave_length; stokes5->current = current; stokes5->depth = depth; stokes5->gravity = gravity; stokes5->kH = k*wave_height; stokes5->eps = eps; stokes5->S = S; stokes5->x0 = x0; stokes5->y0 = y0; // compute coefficients stokes5->A11 = 1/sinh(kd); stokes5->A22 = 3*pow(S,2)/(2*pow(1-S,2)); stokes5->A31 = (-4 -20*S + 10*pow(S,2) - 13*pow(S,3))/(8*sinh(kd)*pow(1-S,3)); stokes5->A33 = (-2*pow(S,2) + 11*pow(S,3))/(8*sinh(kd)*pow(1-S,3)); stokes5->A42 = (12*S - 14*pow(S,2) - 264*pow(S,3) - 45*pow(S,4) - 13*pow(S,5))/(24*pow(1-S,5)); stokes5->A44 = (10*pow(S,3) - 174*pow(S,4) + 291*pow(S,5) + 278*pow(S,6))/(48*(3+2*S)*pow(1-S,5)); stokes5->A51 = (-1184 + 32*S + 13232*pow(S,2) + 21712*pow(S,3) + 20940*pow(S,4) + 12554*pow(S,5) - 500*pow(S,6) - 3341*pow(S,7) - 670*pow(S,8))/(64*sinh(kd)*(3+2*S)*(4+S)*pow(1-S,6)); stokes5->A53 = (4*S + 105*pow(S,2) + 198*pow(S,3) - 1376*pow(S,4) - 1302*pow(S,5) - 117*pow(S,6) + 58*pow(S,7))/(32*sinh(kd)*(3+2*S)*(pow(1-S,6))); stokes5->A55 = (-6*pow(S,3) + 272*pow(S,4) - 1552*pow(S,5) + 852*pow(S,6) + 2029*pow(S,7) + 430*pow(S,8))/(64*sinh(kd)*(3+2*S)*(4+S)*pow(1-S,6)); stokes5->B22 = (cosh(kd)/sinh(kd))*(1+2*S)/(2*(1-S)); stokes5->B31 = -3*(1 + 3*S + 3*pow(S,2) + 2*pow(S,3))/(8*pow(1-S,3)); stokes5->B42 = (cosh(kd)/sinh(kd))*(6 - 26*S - 182*pow(S,2) - 204*pow(S,3) - 25*pow(S,4) + 26*pow(S,5))/(6*(3+2*S)*pow(1-S,4)); stokes5->B44 = (cosh(kd)/sinh(kd))*(24 + 92*S + 122*pow(S,2) + 66*pow(S,3) + 67*pow(S,4) + 34*pow(S,5))/(24*(3+2*S)*pow(1-S,4)); stokes5->B53 = 9*(132 + 17*S - 2216*pow(S,2) - 5897*pow(S,3) - 6292*pow(S,4) - 2687*pow(S,5) + 194*pow(S,6) + 467*pow(S,7) + 82*pow(S,8))/(128*(3+2*S)*(4+S)*pow(1-S,6)); stokes5->B55 = 5*(300 + 1579*S + 3176*pow(S,2) + 2949*pow(S,3) + 1188*pow(S,4) + 675*pow(S,5) + 1326*pow(S,6) + 827*pow(S,7) + 130*pow(S,8))/(384*(3+2*S)*(4+S)*pow(1-S,6)); stokes5->C0 = pow(tanh(kd),0.5); stokes5->C2 = pow(tanh(kd),0.5)*(2+7*pow(S,2))/(4*pow(1-S,2)); stokes5->C4 = pow(tanh(kd),0.5)*(4 + 32*S - 116*pow(S,2) - 400*pow(S,3) - 71*pow(S,4) + 146*pow(S,5))/(32*pow(1-S,5)); stokes5->D2 = -pow(cosh(kd)/sinh(kd),0.5)/2; stokes5->D4 = pow(cosh(kd)/sinh(kd),0.5)*(2 + 4*S + pow(S,2) + 2*pow(S,3))/(8*pow(1-S, 3)); stokes5->E2 = tanh(kd)*(2 + 2*S + 5*pow(S,2))/(4*pow(1-S,2)); stokes5->E4 = tanh(kd)*(8 + 12*S - 152*pow(S,2) - 308*pow(S,3) - 42*pow(S,4) + 77*pow(S,5))/(32*pow(1-S,5)); stokes5->u_mean = (stokes5->C0 + pow(eps,2)*stokes5->C2 + pow(eps,4)*stokes5->C4)/sqrt(k/gravity); stokes5->c = stokes5->u_mean + stokes5->current; // wave celerity stokes5->Q = (stokes5->u_mean)*depth + (stokes5->D2*pow(eps,2) + stokes5->D4*pow(eps,4))/sqrt(pow(k,3)/gravity); stokes5->cs = stokes5->c - (stokes5->Q)/depth; // cs = Stokes drift (set the current to Q/d (Stokes drift=0) for closed wave tanks stokes5->R = 0.5*(gravity/k)*pow(stokes5->C0,2) + kd + (stokes5->E2)*pow(eps,2) + (stokes5->E4)*pow(eps,4); // Bernoulli constant } double eta(Stokes5 *stokes5, double t, double X) /* Given an initialized stokes wave (stokes5), compute the surface elevation eta at time t and coordinate X */ { double keta, kx, eps; eps = stokes5->eps; kx = (stokes5->k)*(X - (stokes5->c)*t - (stokes5->x0)); keta = stokes5->kd + eps*cos(kx) + pow(eps,2)*(stokes5->B22)*cos(2*kx) + pow(eps,3)*(stokes5->B31)*(cos(kx) - cos(3*kx)) + pow(eps,4)*((stokes5->B42)*cos(2*kx) + (stokes5->B44)*cos(4*kx)) + pow(eps,5)*(-((stokes5->B53) + (stokes5->B55))*cos(kx) + (stokes5->B53)*cos(3*kx) + (stokes5->B55)*cos(5*kx)); return stokes5->y0 + keta/(stokes5->k) - stokes5->depth; } double u(Stokes5 *stokes5, double t, double X, double Y) /* Given an initialized stokes wave (stokes5), compute the horizontal velocity at time t and coordinate X, Y. Values above the surface are computed uncritically. */ { double k, kx, ky, eps; eps = stokes5->eps; k = stokes5->k; kx = (stokes5->k)*(X - (stokes5->c)*t - (stokes5->x0)); ky = (stokes5->k)*(Y + stokes5->depth - (stokes5->y0)); return stokes5->c - stokes5->u_mean + (stokes5->C0)*sqrt((stokes5->gravity)/pow(k,3))*(pow(eps,1)*(1*k*(stokes5->A11)*cosh(1*ky)*cos(1*kx)) + pow(eps,2)*(2*k*(stokes5->A22)*cosh(2*ky)*cos(2*kx)) + pow(eps,3)*(1*k*(stokes5->A31)*cosh(1*ky)*cos(1*kx) + 3*k*(stokes5->A33)*cosh(3*ky)*cos(3*kx)) + pow(eps,4)*(2*k*(stokes5->A42)*cosh(2*ky)*cos(2*kx) + 4*k*(stokes5->A44)*cosh(4*ky)*cos(4*kx)) + pow(eps,5)*(1*k*(stokes5->A51)*cosh(1*ky)*cos(1*kx) + 3*k*(stokes5->A53)*cosh(3*ky)*cos(3*kx) + 5*k*(stokes5->A55)*cosh(5*ky)*cos(5*kx))); } double v(Stokes5 *stokes5, double t, double X, double Y) /* Given an initialized stokes wave (stokes5), compute the vertical velocity at time t and coordinate X, Y. Values above the surface are computed, but should later be set to zero. */ { double k, kx, ky, eps; eps = stokes5->eps; k = stokes5->k; kx = (stokes5->k)*(X - (stokes5->c)*t - (stokes5->x0)); ky = (stokes5->k)*(Y + stokes5->depth - (stokes5->y0)); return (stokes5->C0)*sqrt((stokes5->gravity)/pow(k,3))*(pow(eps,1)*(1*k*(stokes5->A11)*sinh(1*ky)*sin(1*kx)) + pow(eps,2)*(2*k*(stokes5->A22)*sinh(2*ky)*sin(2*kx)) + pow(eps,3)*(1*k*(stokes5->A31)*sinh(1*ky)*sin(1*kx) + 3*k*(stokes5->A33)*sinh(3*ky)*sin(3*kx)) + pow(eps,4)*(2*k*(stokes5->A42)*sinh(2*ky)*sin(2*kx) + 4*k*(stokes5->A44)*sinh(4*ky)*sin(4*kx)) + pow(eps,5)*(1*k*(stokes5->A51)*sinh(1*ky)*sin(1*kx) + 3*k*(stokes5->A53)*sinh(3*ky)*sin(3*kx) + 5*k*(stokes5->A55)*sinh(5*ky)*sin(5*kx))); } #endif``````