Coverage for panelaero/VLM.py: 77%
117 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-09-04 09:17 +0000
1#!/usr/bin/env python2
2# -*- coding: utf-8 -*-
3import copy
4import numpy as np
7def calc_induced_velocities(aerogrid, Ma):
8 #
9 # l_2
10 # 4 o---------o 3
11 # | |
12 # u --> b_1 | l k j | b_2
13 # | |
14 # 1 o---------o 2
15 # y l_1
16 # |
17 # z.--- x
19 # define downwash location (3/4 chord and half span of the aero panel)
20 P0 = aerogrid['offset_j']
21 # define vortex location points
22 P1 = aerogrid['offset_P1']
23 P3 = aerogrid['offset_P3']
24 # P2 = mid-point between P1 and P3, not used
25 # normal vector part in vertical direction
26 n_hat_w = np.array(aerogrid['N'][:, 2], ndmin=2).T.repeat(aerogrid['n'], axis=1)
27 # normal vector part in lateral direction
28 n_hat_wl = np.array(aerogrid['N'][:, 1], ndmin=2).T.repeat(aerogrid['n'], axis=1)
30 # divide x coordinates with beta
31 # See Hedman 1965.
32 # However, Hedman divides by beta^2 ... why??
33 beta = (1 - (Ma ** 2.0)) ** 0.5
34 P0[:, 0] = P0[:, 0] / beta
35 P1[:, 0] = P1[:, 0] / beta
36 P3[:, 0] = P3[:, 0] / beta
38 # See Katz & Plotkin, Chapter 10.4.5
39 # get r1,r2,r0
40 r1x = np.array(P0[:, 0], ndmin=2).T - np.array(P1[:, 0], ndmin=2)
41 r1y = np.array(P0[:, 1], ndmin=2).T - np.array(P1[:, 1], ndmin=2)
42 r1z = np.array(P0[:, 2], ndmin=2).T - np.array(P1[:, 2], ndmin=2)
44 r2x = np.array(P0[:, 0], ndmin=2).T - np.array(P3[:, 0], ndmin=2)
45 r2y = np.array(P0[:, 1], ndmin=2).T - np.array(P3[:, 1], ndmin=2)
46 r2z = np.array(P0[:, 2], ndmin=2).T - np.array(P3[:, 2], ndmin=2)
48 # Step 1
49 r1Xr2_x = r1y * r2z - r1z * r2y
50 r1Xr2_y = -r1x * r2z + r1z * r2x # Plus-Zeichen Abweichung zu Katz & Plotkin ??
51 r1Xr2_z = r1x * r2y - r1y * r2x
52 mod_r1Xr2 = (r1Xr2_x ** 2.0 + r1Xr2_y ** 2.0 + r1Xr2_z ** 2.0) ** 0.5
53 # Step 2
54 r1 = (r1x ** 2.0 + r1y ** 2.0 + r1z ** 2.0) ** 0.5
55 r2 = (r2x ** 2.0 + r2y ** 2.0 + r2z ** 2.0) ** 0.5
56 # Step 4
57 r0r1 = (P3[:, 0] - P1[:, 0]) * r1x + (P3[:, 1] - P1[:, 1]) * r1y + (P3[:, 2] - P1[:, 2]) * r1z
58 r0r2 = (P3[:, 0] - P1[:, 0]) * r2x + (P3[:, 1] - P1[:, 1]) * r2y + (P3[:, 2] - P1[:, 2]) * r2z
59 # Step 5
60 gamma = np.ones((aerogrid['n'], aerogrid['n']))
61 D1_base = gamma / 4.0 / np.pi / mod_r1Xr2 ** 2.0 * (r0r1 / r1 - r0r2 / r2)
62 D1_v = r1Xr2_y * D1_base
63 D1_w = r1Xr2_z * D1_base
64 # Step 3
65 epsilon = 10e-6
66 ind = np.where(r1 < epsilon)[0]
67 D1_v[ind] = 0.0
68 D1_w[ind] = 0.0
70 ind = np.where(r2 < epsilon)[0]
71 D1_v[ind] = 0.0
72 D1_w[ind] = 0.0
74 ind = np.where(mod_r1Xr2 < epsilon)
75 D1_v[ind] = 0.0
76 D1_w[ind] = 0.0
78 # get final D1 matrix
79 # D1 matrix contains the perpendicular component of induced velocities at all panels.
80 # For wing panels, it's the z component of induced velocities (D1_w) while for
81 # winglets, it's the y component of induced velocities (D1_v)
82 D1 = D1_w * n_hat_w + D1_v * n_hat_wl
84 # See Katz & Plotkin, Chapter 10.4.7
85 # induced velocity due to inner semi-infinite vortex line
86 d2 = (r1y ** 2.0 + r1z ** 2.0) ** 0.5
87 cosBB1 = 1.0
88 cosBB2 = -r1x / r1
89 cosGamma = r1y / d2
90 sinGamma = -r1z / d2
92 D2_base = -(1.0 / (4.0 * np.pi)) * (cosBB1 - cosBB2) / d2
93 D2_v = sinGamma * D2_base
94 D2_w = cosGamma * D2_base
96 ind = np.where((r1 < epsilon) + (d2 < epsilon))[0]
97 D2_v[ind] = 0.0
98 D2_w[ind] = 0.0
100 D2 = D2_w * n_hat_w + D2_v * n_hat_wl
102 # induced velocity due to outer semi-infinite vortex line
103 d3 = (r2y ** 2.0 + r2z ** 2.0) ** 0.5
104 cosBB1 = r2x / r2
105 cosBB2 = -1.0
106 cosGamma = -r2y / d3
107 sinGamma = r2z / d3
109 D3_base = -(1.0 / (4.0 * np.pi)) * (cosBB1 - cosBB2) / d3
110 D3_v = sinGamma * D3_base
111 D3_w = cosGamma * D3_base
113 ind = np.where((r2 < epsilon) + (d3 < epsilon))[0]
114 D3_v[ind] = 0.0
115 D3_w[ind] = 0.0
117 D3 = D3_w * n_hat_w + D3_v * n_hat_wl
119 return D1, D2, D3
122def mirror_aerogrid_xz(aerogrid):
123 tmp = copy.deepcopy(aerogrid)
124 # mirror y-coord
125 tmp['offset_j'][:, 1] = -tmp['offset_j'][:, 1]
126 tmp['offset_k'][:, 1] = -tmp['offset_k'][:, 1]
127 tmp['offset_l'][:, 1] = -tmp['offset_l'][:, 1]
128 tmp['offset_P1'][:, 1] = -tmp['offset_P1'][:, 1]
129 tmp['offset_P3'][:, 1] = -tmp['offset_P3'][:, 1]
130 tmp['N'][:, 1] = -aerogrid['N'][:, 1]
131 tmp['N'][:, 2] = -aerogrid['N'][:, 2]
132 # assemble both grids
133 aerogrid_xzsym = {'offset_j': np.vstack((aerogrid['offset_j'], tmp['offset_j'])),
134 'offset_k': np.vstack((aerogrid['offset_k'], tmp['offset_l'])),
135 'offset_l': np.vstack((aerogrid['offset_k'], tmp['offset_l'])),
136 'offset_P1': np.vstack((aerogrid['offset_P1'], tmp['offset_P1'])),
137 'offset_P3': np.vstack((aerogrid['offset_P3'], tmp['offset_P3'])),
138 'N': np.vstack((aerogrid['N'], tmp['N'])),
139 'A': np.hstack((aerogrid['A'], tmp['A'])),
140 'l': np.hstack((aerogrid['l'], tmp['l'])),
141 'n': aerogrid['n'] * 2,
142 }
143 return aerogrid_xzsym
146def calc_Ajj(aerogrid, Ma):
147 D1, D2, D3 = calc_induced_velocities(aerogrid, Ma)
148 # define area, chord length and spann of each panel
149 A = aerogrid['A']
150 chord = aerogrid['l']
151 span = A / chord
152 # total D
153 D = D1 + D2 + D3
154 Ajj = D * 0.5 * A / span
155 Bjj = (D2 + D3) * 0.5 * A / span
156 return Ajj, Bjj
159def calc_Qjj(aerogrid, Ma, xz_symmetry=False):
160 '''
161 Symmetry about xz-plane:
162 Only the right hand side is give. The (missing) left hand side is created virtually using mirror_aerogrid_xz().
163 The AIC matrix is calculated for the whole system and can be partitioned into terms:
164 AIC = |RR|LR|
165 |RL|LL|
166 with
167 RR - influence of right side onto right side,
168 LL - influence of left side onto left side,
169 RL and LR - influence from left side onto right side. Due to symmetry, both term are identical.
170 Then, for symmetric motions: AIC_sym = RR - LR
171 And, for asymmetric motions: AIC_asym = RR + LR ??? -> to be checked !!!
172 '''
173 if xz_symmetry:
174 n = aerogrid['n']
175 aerogrid = mirror_aerogrid_xz(aerogrid)
177 # The function calc_Ajj() is Mach number dependent, which involves a scaling of the aerogrid in x-direction.
178 # To make sure that the geometrical scaling has no effect on the following calculations, a 'fresh' a copy of the aerogrid,
179 # created with copy.deepcopy(), is handed over.
180 Ajj, Bjj = calc_Ajj(aerogrid=copy.deepcopy(aerogrid), Ma=Ma)
181 Qjj = -np.linalg.inv(Ajj)
182 if xz_symmetry:
183 return Qjj[0:n, 0:n] - Qjj[n:2 * n, 0:n], Bjj[0:n, 0:n] - Bjj[n:2 * n, 0:n]
184 return Qjj, Bjj
187def calc_Qjjs(aerogrid, Ma, xz_symmetry=False):
188 Qjj = np.zeros((len(Ma), aerogrid['n'], aerogrid['n'])) # dim: Ma,n,n
189 Bjj = np.zeros((len(Ma), aerogrid['n'], aerogrid['n'])) # dim: Ma,n,n
190 for i, i_Ma in enumerate(Ma):
191 Qjj[i, :, :], Bjj[i, :, :] = calc_Qjj(aerogrid, i_Ma, xz_symmetry)
192 return Qjj, Bjj
195def calc_Gamma(aerogrid, Ma, xz_symmetry=False):
196 if xz_symmetry:
197 n = aerogrid['n']
198 aerogrid = mirror_aerogrid_xz(aerogrid)
199 D1, D2, D3 = calc_induced_velocities(aerogrid, Ma)
200 # total D
201 Gamma = -np.linalg.inv((D1 + D2 + D3))
202 Q_ind = D2 + D3
204 if xz_symmetry:
205 return Gamma[0:n, 0:n] - Gamma[n:2 * n, 0:n], Q_ind[0:n, 0:n] - Q_ind[n:2 * n, 0:n]
206 return Gamma, Q_ind
209def calc_Gammas(aerogrid, Ma, xz_symmetry=False):
210 Gamma = np.zeros((len(Ma), aerogrid['n'], aerogrid['n'])) # dim: Ma,n,n
211 Q_ind = np.zeros((len(Ma), aerogrid['n'], aerogrid['n'])) # dim: Ma,n,n
212 for i, i_Ma in enumerate(Ma):
213 Gamma[i, :, :], Q_ind[i, :, :] = calc_Gamma(aerogrid, i_Ma, xz_symmetry)
214 return Gamma, Q_ind