WindWaveShape-DeepWater-Code/NumericPressMagChecker.mpl at master · TZdyrski/WindWaveShape-DeepWater-Code · GitHub

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
#!/home/tzdyrski/bin/maple -q

### Deep water
print(`Deep Water`):

print(`Section 2.4 calculations`):
epsilon := 0.1; # ak
lambda := 2; # meters

rhoAonrhoW := 1.225*10^(-3):
psiP := 3*Pi/4:

#P := 1;
P := epsilon^2;

k := 2*Pi/lambda:
a := epsilon/k:
H := 2*a;

g := 9.8: # meter per second^2
vonKarman := 0.4: # von Karman constant

gammaOnF := 2*Pi*P*sin(psiP):
gamma/f = evalf(gammaOnF);

uStarOnC0 := sqrt(gammaOnF/2/Pi/32.5/rhoAonrhoW):
'u_star'/c_0 = evalf(%);

z0 := H*1200*(H/lambda)^4.5;

U_10 = 'u_star'*ln(10/z0)/vonKarman;

U10OnC0 := uStarOnC0*ln(10/z0)/vonKarman:
U_10/c_0 = evalf(%);

C0 := sqrt(g/k):

U10 := U10OnC0*C0:
U_10 = evalf(%);

restart:

print(`Section 5.3 calculations`):
epsilon := 0.15; # ak
kh := 2.5:

rhoAonrhoW := 1.225*10^(-3):
psiP := 3*Pi/4:

uStarOnC0 := 0.5;
#uStarOnC0 := 1.5;

P := 32.5*rhoAonrhoW/sin(psiP)*uStarOnC0^2:
'P' = evalf(P);
'P' = 'epsilon'^(evalf(ln(P)/ln(epsilon)));

### Shallow water
restart:

print(`Shallow Water`):

print(`Parameters from Numerics`):
epsilon := 0.1; # a/h
mu := 6*epsilon; # (kh)^2

rhoAonrhoW := 1.225*10^(-3):
psiP := Pi/2:

lambda := 20; # wavelength in meters
z := lambda/2; # wind speed measurement height in meters

k := 2*Pi/lambda: # wavenumber in meters^(-1)
kh := sqrt(mu): # kh
h := kh/k; # depth in meters

a := epsilon*h: # amplitude in meters
ak := epsilon*sqrt(mu): # ak
H := 2*a; # wave height in meters

g := 9.8: # acceleration due to gravity in meter per second^2
c0 := sqrt(g*h); # unforced phase speed

# law of the wall height conversion from lambda/2 to z
Z := ln(z*lambda^(4.5)/(1200*H^(3.5)))/ln(lambda^(5.5)/(2400*H^(3.5)));

Pk__on__rhoW__g := 0.25*epsilon: # Pressure magnitude
`Pk/(rho_w*g)` = Pk__on__rhoW__g;

# Shape constant
s := int(diff(sech(x/2)^2,x)^2,x=-infinity..infinity)/int((sech(x/2)^2)^2,x=-infinity..infinity);

print(`Jeffreys Forcing: Donelan 1999 Parameterization`):
S__lambda__on__2 := 0.17; # Jeffreys sheltering coefficient from Donelan 1999

# Calculate wind speed
Uonz := c0*Z*(1+sqrt(s*Pk__on__rhoW__g*sin(psiP)/rhoAonrhoW/S__lambda__on__2)):
cat(`U_`,z) = Uonz;

print(`Jeffreys Forcing: Donelan 1999 Parameterization`):
# Jeffreys sheltering coefficient from Donelan 2006
b := 4.91:
q := 3.98:
if s*Pk__on__rhoW__g < b*rhoAonrhoW then
        print(`Non-separated`):
        print(cat(`s*Pk/(rho_w g) = `,s*Pk__on__rhoW__g,`< b*rho_a/rho_w = `,b*rhoAonrhoW)):

        S__lambda__on__2 := ak*b:
        `S_(lambda/2)` = S__lambda__on__2;
else:
        print(`Separated`):
        print(cat(`s*Pk/(rho_w g) = `,s*Pk__on__rhoW__g,`> b*rho_a/rho_w = `,b*rhoAonrhoW)):

        S__lambda__on__2 := ak*(b-q): # Jeffreys sheltering coefficient from Donelan 2006
        `S_(lambda/2)` = S__lambda__on__2;
end if:

# Calculate wind speed
Uonz := c0*Z*(1+sqrt(s*Pk__on__rhoW__g*sin(psiP)/rhoAonrhoW/S__lambda__on__2)):
cat(`U_`,z) = Uonz;