flaps/analytical_solution.py at main · bartroot/flaps · 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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
###############################################################################

from numba import jit
import numpy as np
import math as math # <-- remove ?

###############################################################################

@jit(nopython=True)
def velocity_x(x,y,R1,R2,exp):
    if exp==0:
       k=4
       r=np.sqrt(x*x+y*y)
       theta=math.atan2(y,x)
       A=2.*(math.log(R1)-math.log(R2))/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       B=(R2**2-R1**2)/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       fr=A*r+B/r
       fpr=A-B/r**2
       gr=A/2.*r + B/r*math.log(r) - 1./r
       hr=(2*gr-fr)/r
       vr=k *gr * math.sin (k * theta)
       vtheta = fr *math.cos(k* theta)
       val=vr*math.cos(theta)-vtheta*math.sin(theta)
    else:
       val=0
    return val

###############################################################################

@jit(nopython=True)
def velocity_y(x,y,R1,R2,exp):
    if exp==0:
       k=4
       r=np.sqrt(x*x+y*y)
       theta=math.atan2(y,x)
       A=2.*(math.log(R1)-math.log(R2))/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       B=(R2**2-R1**2)/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       fr=A*r+B/r
       fpr=A-B/r**2
       gr=A/2.*r + B/r*math.log(r) - 1./r
       hr=(2*gr-fr)/r
       vr=k *gr * math.sin (k * theta)
       vtheta = fr *math.cos(k* theta)
       val=vr*math.sin(theta)+vtheta*math.cos(theta)
    else:
       val=0
    return val

###############################################################################

@jit(nopython=True)
def pressure(x,y,R1,R2,rho0,g0,exp):
    r=np.sqrt(x*x+y*y)
    if exp==0:
       k=4
       theta=math.atan2(y,x)
       A=2.*(math.log(R1)-math.log(R2))/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       B=(R2**2-R1**2)/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       fr=A*r+B/r
       gr=A/2.*r + B/r*math.log(r) - 1./r
       hr=(2*gr-fr)/r
       val=k*hr*math.sin(k*theta) + rho0*g0*(r-R2)
    else:
       val=rho0*g0*(R2-r)
    return val

###############################################################################

@jit(nopython=True)
def sr_xx(x,y,R1,R2,exp):
    if exp==0:
       k=4
       r=np.sqrt(x*x+y*y)
       theta=math.atan2(y,x)
       A=2.*(math.log(R1)-math.log(R2))/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       B=(R2**2-R1**2)/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       gr=A/2.*r + B/r*math.log(r) - 1./r
       gpr=A/2 + B*((1-math.log(r)) / r**2 ) +1./r**2
       fr=A*r+B/r
       fpr=A-B/r**2
       err=gpr*k*math.sin(k*theta)
       ert=0.5*(k**2/r*gr+fpr-fr/r)*math.cos(k*theta)
       ett=(gr-fr)/r*k*math.sin(k*theta)
       val=err*(math.cos(theta))**2\
          +ett*(math.sin(theta))**2\
          -2*ert*math.sin(theta)*math.cos(theta)
    else:
       val=0
    return val

###############################################################################

@jit(nopython=True)
def sr_yy(x,y,R1,R2,exp):
    if exp==0:
       k=4
       r=np.sqrt(x*x+y*y)
       theta=math.atan2(y,x)
       A=2.*(math.log(R1)-math.log(R2))/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       B=(R2**2-R1**2)/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       gr=A/2.*r + B/r*math.log(r) - 1./r
       gpr=A/2 + B*((1-math.log(r)) / r**2 ) +1./r**2
       fr=A*r+B/r
       fpr=A-B/r**2
       err=gpr*k*math.sin(k*theta)
       ert=0.5*(k**2/r*gr+fpr-fr/r)*math.cos(k*theta)
       ett=(gr-fr)/r*k*math.sin(k*theta)
       val=err*(math.sin(theta))**2\
          +ett*(math.cos(theta))**2\
          +2*ert*math.sin(theta)*math.cos(theta)
    else:
       val=0
    return val

###############################################################################

@jit(nopython=True)
def sr_xy(x,y,R1,R2,exp):
    if exp==0:
       k=4
       r=np.sqrt(x*x+y*y)
       theta=math.atan2(y,x)
       A=2.*(math.log(R1)-math.log(R2))/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       B=(R2**2-R1**2)/(R2**2*math.log(R1)-R1**2*math.log(R2) )
       gr=A/2.*r + B/r*math.log(r) - 1./r
       gpr=A/2 + B*((1-math.log(r)) / r**2 ) +1./r**2
       fr=A*r+B/r
       fpr=A-B/r**2
       err=gpr*k*math.sin(k*theta)
       ert=0.5*(k**2/r*gr+fpr-fr/r)*math.cos(k*theta)
       ett=(gr-fr)/r*k*math.sin(k*theta)
       val=ert*(math.cos(theta)**2-math.sin(theta)**2)\
          +(err-ett)*math.cos(theta)*math.sin(theta)
    else:
       val=0
    return val

###############################################################################