diff --git a/doc/jcl_template.py b/doc/jcl_template.py index 03b0eb5c..cdc27c47 100755 --- a/doc/jcl_template.py +++ b/doc/jcl_template.py @@ -14,7 +14,7 @@ class jcl: def __init__(self): # Give your aircraft a name and set some general parameters - self.general = {'aircraft': 'DLR F-19-S', + self.general = {'aircraft': 'MyAircraft', # Reference span width (from tip to tip) 'b_ref': 15.375, # Reference chord length @@ -30,7 +30,7 @@ def __init__(self): to be implemented as a python module. """ # Electronic flight control system - self.efcs = {'version': 'mephisto', # Name of the corresponding python module + self.efcs = {'version': 'MyEFCS', # Name of the corresponding python module # Path where to find the EFCS module 'path': '/path/to/EFCS', } @@ -70,10 +70,11 @@ def __init__(self): self.aero = {'method': 'mona_steady', # 'mona_steady' - steady trim and quasi-steady time domain simulations # 'mona_unsteady' - unsteady time domain simulation based on the RFA, e.g. for gust - # 'freq_dom' - frequency domain simulations, e.g. gust, continuous turbulence, flutter, etc + # 'mona_freq_dom' - frequency domain simulations, e.g. gust, continuous turbulence, flutter, etc # 'nonlin_steady' - steady trim and quasi-steady time domain simulations with some non-linearities # 'cfd_steady' - steady trim # 'cfd_unsteady' - unsteady time domain simulation, e.g. for gust + # 'cfd_freq_dom' - frequency domain simulations, e.g. flutter # # True or False, aerodynamic feedback of elastic structure on aerodynamics can be deactivated. # You will still see deformations, but there is no coupling. @@ -101,12 +102,14 @@ def __init__(self): # number of poles for rational function approximation (RFA) 'n_poles': 4, # Additional parameters for CFD - 'para_path': '/scratch/tau/', - 'para_file': 'para', + 'para_path': '/path/to/CFD-workspace', + 'para_file': 'my_baseline_para_file', # Currently implemented interfaces: 'tau' or 'su2' 'cfd_solver': 'tau', 'tau_solver': 'el', 'tau_cores': 16, + # Name of job that was used to precompute CFD-based GAFs + 'job_name_gafs': 'jcl_dc3_gafs', # --- Start of experimental section, only for special cases --- # Correction coefficient at CG, negativ = destabilizing 'Cn_beta_corr': [-0.012], @@ -126,7 +129,7 @@ def __init__(self): # General CFD surface mesh information self.meshdefo = {'surface': {'fileformat': 'netcdf', # implemented file formats: 'cgns', 'netcdf', 'su2' # file name of the CFD mesh - 'filename_grid': 'tau.grid', + 'filename_grid': 'CFD-mesh.grid', # list of markers [1, 2, ...] or ['upper', 'lower', ...] of surfaces to be included in # deformation 'markers': [1, 3], @@ -211,8 +214,8 @@ def __init__(self): 'n_blades': [2, 2], # Mach number for VLM4Prop 'Ma': [0.25], - # Input-file ('.yaml') for PyPropMAt and VLM4Prop - 'propeller_input_file': 'HAP_O6_PROP_pitch.yaml', + # Input-file ('.yaml') for PyPropMat and VLM4Prop + 'propeller_input_file': 'Prop.yaml', } # CFD-specific # In case a pressure inlet boundary is modeled in the CFD mesh, the boundary condition will be @@ -338,10 +341,17 @@ def __init__(self): # Flutter parameters for k and ke method 'flutter_para': {'method': 'k', 'k_red': np.linspace(2.0, 0.001, 1000)}, # Flutter parameters for pk method - # There are two implementations of the PK method: 'pk_schwochow', 'pk_rodden' - # Available mode tracking algortihms: 'MAC', 'MAC*PCC' (recommended), 'MAC*HDM' - # 'flutter_para': {'method': 'pk', 'Vtas': np.linspace(100.0, 500.0, 100), - # 'tracking': 'MAC*PCC'}, + # 'flutter_para': {'method': 'pk_schwochow', # Available implementations: 'pk_schwochow', 'pk_rodden' + # 'Vtas': np.linspace(100.0, 500.0, 100), + # 'tracking': 'MAC*PCC', # Available: 'MAC', 'MAC*PCC' (recommended), 'MAC*HDM' + # }, + # True or False, enables computation of CFD-based generalized aerodynamic forces (GAFs) + 'gaf': True, + 'gaf_para': {'df': 1.0, # The frequency resolution governs the time length. + # The max. frequency governs the dt of the time domain simulation. + # Rule of thumb: 20x the highest frequency of interest. + 'fmax': 500.0 + }, }, ] # End diff --git a/loadskernel/cfd_interfaces/meshdefo.py b/loadskernel/cfd_interfaces/meshdefo.py index 12e5fd8b..70cb18f6 100644 --- a/loadskernel/cfd_interfaces/meshdefo.py +++ b/loadskernel/cfd_interfaces/meshdefo.py @@ -10,26 +10,30 @@ class Meshdefo(): """ def Ux2(self, Ux2): - logging.info('Apply control surface deflections to cfdgrid.') - Ujx2 = np.zeros(self.aerogrid['n'] * 6) - if 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'y': - hingeline = 'y' - elif 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'z': - hingeline = 'z' - else: # default - hingeline = 'y' - for i_x2 in range(len(self.x2grid['key'])): - logging.debug('Apply deflection of {} for {:0.4f} [deg].'.format( - self.x2grid['key'][i_x2], Ux2[i_x2] / np.pi * 180.0)) - if hingeline == 'y': - Ujx2 += np.dot(self.Djx2[i_x2], [0, 0, 0, 0, Ux2[i_x2], 0]) - elif hingeline == 'z': - Ujx2 += np.dot(self.Djx2[i_x2], [0, 0, 0, 0, 0, Ux2[i_x2]]) - self.transfer_deformations(self.aerogrid, Ujx2, '_k', rbf_type='wendland2', surface_spline=False, support_radius=1.5) + if np.any(Ux2): + logging.debug('Apply control surface deflections to cfdgrid.') + Ujx2 = np.zeros(self.aerogrid['n'] * 6) + if 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'y': + hingeline = 'y' + elif 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'z': + hingeline = 'z' + else: # default + hingeline = 'y' + for i_x2 in range(len(self.x2grid['key'])): + logging.debug('Apply deflection of {} for {:0.4f} [deg].'.format( + self.x2grid['key'][i_x2], Ux2[i_x2] / np.pi * 180.0)) + if hingeline == 'y': + Ujx2 += np.dot(self.Djx2[i_x2], [0, 0, 0, 0, Ux2[i_x2], 0]) + elif hingeline == 'z': + Ujx2 += np.dot(self.Djx2[i_x2], [0, 0, 0, 0, 0, Ux2[i_x2]]) + self.transfer_deformations(self.aerogrid, Ujx2, '_k', rbf_type='wendland2', + surface_spline=False, support_radius=1.5) + else: + logging.debug('Apply NO control surface deflections to cfdgrid.') def Uf(self, Uf): - if 'flex' in self.jcl.aero and self.jcl.aero['flex']: - logging.info('Apply flexible deformations to cfdgrid.') + if 'flex' in self.jcl.aero and self.jcl.aero['flex'] and np.any(Uf): + logging.debug('Apply flexible deformations to cfdgrid.') # set-up spline grid if self.jcl.spline['splinegrid']: # make sure that there are no double points in the spline grid as this would cause a singularity of the @@ -45,4 +49,4 @@ def Uf(self, Uf): self.transfer_deformations(splinegrid, Ug_f_body, '', rbf_type='tps', surface_spline=False) else: - logging.info('Apply NO flexible deformations to cfdgrid.') + logging.debug('Apply NO flexible deformations to cfdgrid.') diff --git a/loadskernel/cfd_interfaces/su2_interface.py b/loadskernel/cfd_interfaces/su2_interface.py index f88d5701..05b38652 100644 --- a/loadskernel/cfd_interfaces/su2_interface.py +++ b/loadskernel/cfd_interfaces/su2_interface.py @@ -113,7 +113,7 @@ def set_deformations(self): Communicate the change of coordinates of the fluid interface to the fluid solver. Prepare the fluid solver for mesh deformation. """ - logging.info('Sending surface deformations to SU2.') + logging.debug('Sending surface deformations to SU2.') for x in range(self.local_mesh['n']): self.FluidSolver.SetMarkerCustomDisplacement(self.local_mesh['MarkerID'][x], self.local_mesh['VertexIndex'][x], @@ -199,7 +199,7 @@ def init_solver(self): self.get_local_mesh() def run_solver(self, i_timestep=0): - logging.info('Waiting until all processes are ready to perform a coordinated start...') + logging.debug('Waiting until all processes are ready to perform a coordinated start...') self.comm.barrier() logging.info('Launch SU2 for time step {}.'.format(i_timestep)) # start timer @@ -510,16 +510,23 @@ def update_timedom_para(self): os.symlink(filename_steady, filename_unsteady1) except FileExistsError: pass - - """ - In the time domain, simply rely on the the density residual to establish convergence. - This is because SU2 only needs a few inner iterations per time step, which are too few for a meaningful - cauchy convergence. - """ + # Set-up inner iterations + if 'gaf' in self.simcase and self.simcase['gaf']: + # In case of GAF computation, no convergence criterion can be used because the reference / zero solution + # (without a pulse) has to have exactly the same number of steps as the pulse solution. + # Possibly, the number of inner iterations depends on the aircraft, mesh size, etc. + config['INNER_ITER'] = 4 + if 'CONV_RESIDUAL_MINVAL' in config: + config.pop('CONV_RESIDUAL_MINVAL') + else: + # In the time domain, simply rely on the the density residual to establish convergence. + # This is because SU2 only needs a few inner iterations per time step, which are too few for a meaningful + # cauchy convergence. + config['INNER_ITER'] = 30 + config['CONV_RESIDUAL_MINVAL'] = -6 + # Remove other convergence criteria if 'CONV_FIELD' in config: config.pop('CONV_FIELD') - config['INNER_ITER'] = 30 - config['CONV_RESIDUAL_MINVAL'] = -6 # There is no need for restart solutions, they only take up storage space. Write plotting files only. config['OUTPUT_FILES'] = ['TECPLOT', 'SURFACE_TECPLOT'] @@ -561,7 +568,6 @@ def update_gust_para(self, Vtas, Vgust): # Note: In SU2 this is the full gust length, not the gust gradient H (half gust length). config['GUST_WAVELENGTH'] = 2.0 * self.simcase['gust_gradient'] config['GUST_PERIODS'] = 1.0 - config['GUST_AMPL'] = Vgust config['GUST_BEGIN_TIME'] = 0.0 config['GUST_BEGIN_LOC'] = -2.0 * self.simcase['gust_gradient'] - self.simcase['gust_para']['T1'] * Vtas diff --git a/loadskernel/equations/cfd_frequency_domain.py b/loadskernel/equations/cfd_frequency_domain.py new file mode 100644 index 00000000..cd67b204 --- /dev/null +++ b/loadskernel/equations/cfd_frequency_domain.py @@ -0,0 +1,129 @@ +import logging +import numpy as np + +from scipy.interpolate import interp1d +from scipy.fftpack import fft + +from loadskernel.interpolate import MatrixInterpolation +from loadskernel.equations.mona_frequency_domain import KMethod as MonaKMethod +from loadskernel.equations.mona_frequency_domain import KEMethod as MonaKEMethod +from loadskernel.equations.mona_frequency_domain import PKMethodRodden as MonaPKMethodRodden +from loadskernel.equations.mona_frequency_domain import GustExcitation as MonaGustExcitation + + +class GustExcitation(MonaGustExcitation): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + # Initialize additional attributes to avoid defining them outside __init__ + # Interpolators + self.Qhh_interp = None + self.Qhk_interp = None + self.Qk_gust_interp = None + + def build_AIC_interpolators(self): + # Similar as in the fultter solutions, but not scaled with the dynamic pressure q_dyn. + Qhh = [] + for i_k, _ in enumerate(self.GAFs['k_red']): + Qhh.append(self.PHIkh.T.dot(self.GAFs['Qhk'][:, :, i_k])) + self.Qhh_interp = MatrixInterpolation(self.GAFs['k_red'], Qhh) + # Reorder matrices (freq, aero panels, modes) because the interpolator works along the first axis. + Qhk = np.moveaxis(self.GAFs['Qhk'], -1, 0) + Qk_gust = np.moveaxis(self.GAFs['Qk_gust'], -1, 0) + self.Qhk_interp = MatrixInterpolation(self.GAFs['k_red'], Qhk) + self.Qk_gust_interp = MatrixInterpolation(self.GAFs['k_red'], Qk_gust) + + def transfer_function(self, f): + omega = 2.0 * np.pi * f + Qhh = self.Qhh_interp(self.f2k(f)) + TF = np.linalg.inv(-self.Mhh * omega ** 2 + complex(0, 1) * omega * self.Dhh + self.Khh - Qhh) + return TF + + def calc_gust_excitation(self, freqs, t): + gust_f = fft(self.one_m_cosine_gust(t)) + Ph_fourier = np.zeros((self.n_modes, len(freqs)), dtype='complex128') + Pk_fourier = np.zeros((self.aerogrid['n'] * 6, len(freqs)), dtype='complex128') + for i, f in enumerate(freqs): + Qk_gust = self.Qk_gust_interp(self.f2k(f)) + Pk_fourier[:, i] = Qk_gust.dot(gust_f[i]) + Ph_fourier[:, i] = self.PHIkh.T.dot(Pk_fourier[:, i]) + return Ph_fourier, Pk_fourier + + def one_m_cosine_gust(self, t): + # This is "only" the gust signal; unlike with panel methods, there is no relationship with the aicraft geometry here. + # The effect of the aircraft penetrating into the gust was already captured during the GAF computations. + tw = self.simcase['gust_gradient'] * 2.0 / self.Vtas + gust = self.Vtas * self.WG_TAS * 0.5 * (1 - np.cos(2.0 * np.pi * t / tw)) + gust[np.where(t > tw)] = 0.0 + return gust + + def calc_aero_response(self, freqs, Uh, dUh_dt): + # Because the motion is included in the GAF computations, matrix Qhk is multiplied "only" by the generalized + # deformations Uh. + Ph_fourier = np.zeros((self.n_modes, len(freqs)), dtype='complex128') + Pk_fourier = np.zeros((self.aerogrid['n'] * 6, len(freqs)), dtype='complex128') + for i, f in enumerate(freqs): + Qhk = self.Qhk_interp(self.f2k(f)) + Pk_fourier[:, i] = Qhk.dot(Uh[:, i]) + Ph_fourier[:, i] = self.PHIkh.T.dot(Pk_fourier[:, i]) + return Ph_fourier, Pk_fourier + + +class KMethod(MonaKMethod): + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + # initialize frequently used attributes to satisfy linters/static analyzers + self.n_freqs = None + self.n_modes = None + self.k_reds = None + + def build_AIC_interpolators(self): + Qhh = [] + for i_k, _ in enumerate(self.GAFs['k_red']): + Qhh.append(self.PHIkh.T.dot(self.GAFs['Qhk'][:, :, i_k]) / self.GAFs['q_dyn']) + self.Qhh_interp = interp1d(self.GAFs['k_red'], Qhh, kind='cubic', axis=0, fill_value="extrapolate") + + def setup_frequence_parameters(self): + self.n_modes = self.model['mass'][self.trimcase['mass']]['n_modes'][()] + 5 + self.k_reds = self.simcase['flutter_para']['k_red'] + self.n_freqs = len(self.k_reds) + + if self.k_reds.max() > np.max(self.GAFs['k_red']): + logging.warning('Required reduced frequency = %0.3f but GAFs given only up to %0.3f', + self.k_reds.max(), np.max(self.GAFs['k_red'])) + + +class KEMethod(MonaKEMethod, KMethod): + """ + The CFD-based KE-Method uses the combined formulations of the CFD-based K-Method (see above) + and the Mona-based KE-Method (imported as MonaKMethod). This is achieved by inheriting twice. + """ + + +class PKMethodRodden(MonaPKMethodRodden): + + def build_AIC_interpolators(self): + # Same formulation as in K-Method, but with custom, linear matrix interpolation + Qhh = [] + for i_k, _ in enumerate(self.GAFs['k_red']): + Qhh.append(self.PHIkh.T.dot(self.GAFs['Qhk'][:, :, i_k]) / self.GAFs['q_dyn']) + self.Qhh_interp = MatrixInterpolation(self.GAFs['k_red'], Qhh) + + def system(self, k_red): + rho = self.atmo['rho'] + # Make sure that k_red is not zero due to the division by k_red. If k_red=0.0, set to a small value. + # This line is the only difference to the mona-based PKMethodRodden, because GAFs from CFD are also + # calculated for k_red=0.0. + k_red = np.max([k_red, 0.001]) + + Qhh = self.Qhh_interp(k_red) + Mhh_inv = np.linalg.inv(self.Mhh) + + upper_part = np.concatenate((np.zeros((self.n_modes, self.n_modes)), + np.eye(self.n_modes)), axis=1) + lower_part = np.concatenate((-Mhh_inv.dot(self.Khh - rho / 2 * self.Vtas ** 2.0 * Qhh.real), + -Mhh_inv.dot(self.Dhh - rho / 4 * self.Vtas * self.macgrid['c_ref'] / k_red * Qhh.imag)), + axis=1) + A = np.concatenate((upper_part, lower_part)) + return A diff --git a/loadskernel/equations/cfd.py b/loadskernel/equations/cfd_time_domain.py similarity index 98% rename from loadskernel/equations/cfd.py rename to loadskernel/equations/cfd_time_domain.py index fc5f444e..62c0969a 100644 --- a/loadskernel/equations/cfd.py +++ b/loadskernel/equations/cfd_time_domain.py @@ -1,6 +1,6 @@ import numpy as np -from loadskernel.equations.steady import Steady +from loadskernel.equations.mona_time_domain import Steady from loadskernel.solution_tools import gravitation_on_earth @@ -127,7 +127,7 @@ def equations(self, X, t, modus): # aerodynamics self.cfd_interface.update_general_para() self.cfd_interface.update_timedom_para() - if self.simcase['gust']: + if 'gust' in self.simcase and self.simcase['gust']: self.cfd_interface.update_gust_para(Vtas, self.WG_TAS * Vtas) self.cfd_interface.init_solver() self.cfd_interface.set_euler_transformation(delta_XYZ, PhiThetaPsi) diff --git a/loadskernel/equations/common.py b/loadskernel/equations/common.py index 39bfbb50..1366b4c1 100644 --- a/loadskernel/equations/common.py +++ b/loadskernel/equations/common.py @@ -38,6 +38,7 @@ def __init__(self, solution, X0=None): self.defo_old = 0.0 # load data needed for subsequent simulation self.load_data() + self.load_GAFs() # set-up simulation parameters self.setup_hingeline() self.setup_efcs() @@ -93,6 +94,20 @@ def load_data(self): self.Qjj = self.aero['Qjj'] + def load_GAFs(self): + if self.jcl.aero['method'] in ['cfd_freq_dom']: + # Derive the paht/key from the trimcase description + key = '.'.join(self.trimcase['desc'].split('.')[:-1]) + if key in self.model['GAFs']: + self.GAFs = load_hdf5_dict(self.model['GAFs'][key]) + self.X0 = self.GAFs['X0'] + else: + logging.error('No GAFs found for "%s" in model!', key) + msg = """GAFs depend on the operational point, mass case, etc. They are organized by the trimcase description + 'desc' in the JCL, e.g. 'CC.M1.OVCFL000.xxx'. Please make sure that the GAFs for the requested + trimcase have bee computed and gathered/added to the model file.""" + logging.info(msg) + def setup_hingeline(self): # set hingeline for cs deflections if 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'z': @@ -207,10 +222,10 @@ def setup_controller(self): self.trimcond_Y[np.where(self.trimcond_Y[:, 0] == 'Nz')[0][0], 2])]), }) else: - logging.error('Unknown EFCS: {}'.format(self.jcl.efcs['version'])) + logging.error('Unknown EFCS: %s', self.jcl.efcs['version']) def setup_aero_matrices(self): - if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'freq_dom']: + if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'freq_dom', 'mona_freq_dom']: self.Djf_1 = self.aerogrid['Nmat'].dot(self.aerogrid['Rmat'].dot(self.PHIjf)) self.Djf_2 = self.aerogrid['Nmat'].dot(self.PHIjf) * -1.0 self.Djh_1 = self.aerogrid['Nmat'].dot(self.aerogrid['Rmat'].dot(self.PHIjh)) diff --git a/loadskernel/equations/landing.py b/loadskernel/equations/landing.py deleted file mode 100644 index 8dc50fd9..00000000 --- a/loadskernel/equations/landing.py +++ /dev/null @@ -1,153 +0,0 @@ -import numpy as np -from loadskernel.equations.common import Common -from loadskernel.solution_tools import gravitation_on_earth - - -class Landing(Common): - - def equations(self, X, t, modus): - self.counter += 1 - # recover states - Tgeo2body, Tbody2geo = self.geo2body(X) - dUcg_dt, Uf, dUf_dt = self.recover_states(X) - Vtas, q_dyn = self.recover_Vtas(X) - onflow = self.recover_onflow(X) - alpha, beta, gamma = self.windsensor(X, Vtas, Uf, dUf_dt) - Ux2 = self.get_Ux2(X) - - # aerodynamics - Pk_rbm, wj_rbm = self.rbm(onflow, q_dyn, Vtas) - Pk_cam, wj_cam = self.camber_twist(q_dyn) - Pk_cs, wj_cs = self.cs(X, Ux2, q_dyn) - Pk_f, wj_f = self.flexible(Uf, dUf_dt, onflow, q_dyn, Vtas) - Pk_gust, wj_gust = self.gust(X, q_dyn) - - wj = wj_rbm + wj_cam + wj_cs + wj_f + wj_gust - Pk_idrag = self.idrag(wj, q_dyn) - Pk_unsteady = Pk_rbm * 0.0 - - # correction coefficients - Pb_corr = self.correctioon_coefficients(alpha, beta, q_dyn) - - # landing gear - Pextra, _, dp2, ddp2, F1, F2 = self.landinggear(X, Tbody2geo) - - # summation of forces - Pk_aero = Pk_rbm + Pk_cam + Pk_cs + Pk_f + Pk_gust + Pk_idrag + Pk_unsteady - Pmac = np.dot(self.Dkx1.T, Pk_aero) - Pb = np.dot(self.PHImac_cg.T, Pmac) + Pb_corr + np.dot(self.PHIextra_cg.T, Pextra) - - g_cg = gravitation_on_earth(self.PHInorm_cg, Tgeo2body) - - # EoM - d2Ucg_dt2, Nxyz = self.rigid_EoM(dUcg_dt, Pb, g_cg, modus) - Pf = np.dot(self.PHIkf.T, Pk_aero) + self.Mfcg.dot( - np.hstack((d2Ucg_dt2[0:3] - g_cg, d2Ucg_dt2[3:6]))) + np.dot(self.PHIf_extra, Pextra) # viel schneller! - d2Uf_dt2 = self.flexible_EoM(dUf_dt, Uf, Pf) - - # CS derivatives - dcommand = self.get_command_derivatives(t, X, Vtas, gamma, alpha, beta, Nxyz, np.dot(Tbody2geo, X[6:12])[0:3]) - - # output - if modus in ['trim', 'trim_full_output']: - Y = np.hstack((np.dot(Tbody2geo, X[6:12]), - np.dot(self.PHIcg_norm, d2Ucg_dt2), - dUf_dt, - d2Uf_dt2, - dcommand, - Nxyz[2], - Vtas, - beta, - )) - elif modus in ['sim', 'sim_full_output']: - Y = np.hstack((np.dot(Tbody2geo, X[6:12]), - np.dot(self.PHIcg_norm, d2Ucg_dt2), - dUf_dt, - d2Uf_dt2, - dcommand, - np.hstack((dp2, ddp2)), - Nxyz[2], - Vtas, - beta, - )) - if modus in ['trim', 'sim']: - return Y - elif modus in ['trim_full_output', 'sim_full_output']: - # calculate translations, velocities and accelerations of some additional points - # (might also be used for sensors in a closed-loop system - # position LG attachment point over ground - p1 = -self.cggrid['offset'][:, 2] + self.extragrid['offset'][:, 2] \ - + self.PHIextra_cg.dot(np.dot(self.PHInorm_cg, X[0:6]))[self.extragrid['set'][:, 2]] \ - + self.PHIf_extra.T.dot(X[12:12 + self.n_modes])[self.extragrid['set'][:, 2]] - # velocity LG attachment point - dp1 = self.PHIextra_cg.dot(np.dot(self.PHInorm_cg, np.dot(Tbody2geo, X[6:12])))[self.extragrid['set'][:, 2]] \ - + self.PHIf_extra.T.dot(X[12 + self.n_modes:12 + self.n_modes * 2])[self.extragrid['set'][:, 2]] - # acceleration LG attachment point - ddp1 = self.PHIextra_cg.dot(np.dot(self.PHInorm_cg, np.dot(Tbody2geo, Y[6:12])))[self.extragrid['set'][:, 2]] \ - + self.PHIf_extra.T.dot(Y[12 + self.n_modes:12 + self.n_modes * 2])[self.extragrid['set'][:, 2]] - - response = {'X': X, - 'Y': Y, - 't': np.array([t]), - 'Pk_rbm': Pk_rbm, - 'Pk_cam': Pk_cam, - 'Pk_aero': Pk_aero, - 'Pk_cs': Pk_cs, - 'Pk_f': Pk_f, - 'Pk_gust': Pk_gust, - 'Pk_unsteady': Pk_unsteady, - 'Pk_idrag': Pk_idrag, - 'q_dyn': np.array([q_dyn]), - 'Pb': Pb, - 'Pmac': Pmac, - 'Pf': Pf, - 'alpha': np.array([alpha]), - 'beta': np.array([beta]), - 'Ux2': Ux2, - 'dUcg_dt': dUcg_dt, - 'd2Ucg_dt2': d2Ucg_dt2, - 'Uf': Uf, - 'dUf_dt': dUf_dt, - 'd2Uf_dt2': d2Uf_dt2, - 'Nxyz': Nxyz, - 'g_cg': g_cg, - 'Pextra': Pextra, - 'p1': p1, - 'dp1': dp1, - 'ddp1': ddp1, - 'F1': F1, - 'F2': F2, - } - return response - - def eval_equations(self, X_free, time, modus='trim_full_output'): - - if modus in ['trim', 'trim_full_output']: - # get inputs from trimcond and apply inputs from fsolve - X = np.array(self.trimcond_X[:, 2], dtype='float') - X[np.where((self.trimcond_X[:, 1] == 'free'))[0]] = X_free - elif modus in ['sim', 'sim_full_output']: - X = X_free - - # evaluate model equations - if modus == 'trim': - Y = self.equations(X, time, 'trim') - # get the current values from Y and substract tamlab.figure() - # fsolve only finds the roots; Y = 0 - Y_target_ist = Y[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] - Y_target_soll = np.array(self.trimcond_Y[:, 2], dtype='float')[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] - out = Y_target_ist - Y_target_soll - return out - - elif modus == 'sim': - Y = self.equations(X, time, 'sim') - return Y[self.solution.idx_state_derivatives + self.solution.idx_input_derivatives - + self.solution.idx_lg_derivatives] # Nz ist eine Rechengroesse und keine Simulationsgroesse! - - elif modus == 'sim_full_output': - response = self.equations(X, time, 'sim_full_output') - return response - - elif modus == 'trim_full_output': - response = self.equations(X, time, 'trim_full_output') - return response diff --git a/loadskernel/equations/frequency_domain.py b/loadskernel/equations/mona_frequency_domain.py similarity index 89% rename from loadskernel/equations/frequency_domain.py rename to loadskernel/equations/mona_frequency_domain.py index 05445162..3e347eb7 100644 --- a/loadskernel/equations/frequency_domain.py +++ b/loadskernel/equations/mona_frequency_domain.py @@ -15,50 +15,70 @@ class GustExcitation(Common): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + # Initialize additional attributes to avoid defining them outside __init__ + self.n_modes = None + self.Vtas = None + self.q_dyn = None + self.fmax = None + self.n_freqs = None + self.t = None + self.idx_output = None + self.freqs = None + self.fftomega = None + self.positiv_fftfreqs = None + self.positiv_fftomega = None + self.k_reds = None # For compatibility with KMethod + # Interpolators + self.Qjj_interp = None + self.Qhh_1_interp = None + self.Qhh_2_interp = None + def eval_equations(self): self.setup_frequence_parameters() - logging.info('building transfer functions') + logging.info('Building transfer functions') self.build_AIC_interpolators() # unsteady positiv_TFs = self.build_transfer_functions(self.positiv_fftfreqs) TFs = np.zeros((self.n_modes, self.n_modes, self.n_freqs), dtype='complex128') for i_mode in range(self.n_modes): TFs[:, i_mode, :] = self.mirror_fouriersamples_even(positiv_TFs[:, i_mode, :]) - logging.info('calculating gust excitation (in physical coordinates, this may take considerable time and memory)') + logging.info('Calculating gust excitation (in physical coordinates, this may take considerable time and memory)') Ph_gust_fourier, Pk_gust_fourier = self.calc_gust_excitation(self.positiv_fftfreqs, self.t) Ph_gust_fourier = self.mirror_fouriersamples_even(Ph_gust_fourier) Pk_gust_fourier = self.mirror_fouriersamples_even(Pk_gust_fourier) - logging.info('calculating responses') + logging.info('Calculating aircraft response') Uh_fourier = TFs * Ph_gust_fourier # [Antwort, Anregung, Frequenz] Uh = ifft(np.array((Uh_fourier) * (1j * self.fftomega) ** 0).sum(axis=1)) dUh_dt = ifft(np.array((Uh_fourier) * (1j * self.fftomega) ** 1).sum(axis=1)) d2Uh_dt2 = ifft(np.array((Uh_fourier) * (1j * self.fftomega) ** 2).sum(axis=1)) - logging.info('reconstructing aerodynamic forces (in physical coordinates, this may take considerable time and memory)') + logging.info('Reconstructing aerodynamic forces (in physical coordinates, this may take considerable time and memory)') Ph_aero_fourier, Pk_aero_fourier = self.calc_aero_response( self.positiv_fftfreqs, np.array((Uh_fourier) * (1j * self.fftomega) ** 0).sum(axis=1)[:, :self.n_freqs // 2 + 1], np.array((Uh_fourier) * (1j * self.fftomega) ** 1).sum(axis=1)[:, :self.n_freqs // 2 + 1],) Ph_aero_fourier = self.mirror_fouriersamples_even(Ph_aero_fourier) Pk_aero_fourier = self.mirror_fouriersamples_even(Pk_aero_fourier) - Pk_aero = np.real(ifft(Pk_gust_fourier) + ifft(Pk_aero_fourier))[:, self.t_output] - Pk_gust = np.real(ifft(Pk_gust_fourier))[:, self.t_output] + Pk_aero = np.real(ifft(Pk_gust_fourier) + ifft(Pk_aero_fourier))[:, self.idx_output] + Pk_gust = np.real(ifft(Pk_gust_fourier))[:, self.idx_output] # split h-set into b- and f-set # remember that the x-component was omitted - Ucg = np.concatenate((np.zeros((len(self.t_output), 1)), Uh[:5, self.t_output].T.real - Uh[:5, 0].real), axis=1) - dUcg_dt = np.concatenate((np.zeros((len(self.t_output), 1)), - dUh_dt[:5, self.t_output].T.real - dUh_dt[:5, 0].real), axis=1) - d2Ucg_dt2 = np.concatenate((np.zeros((len(self.t_output), 1)), - d2Uh_dt2[:5, self.t_output].T.real - d2Uh_dt2[:5, 0].real), axis=1) - Uf = Uh[5:, self.t_output].T.real - Uh[5:, 0].real - dUf_dt = dUh_dt[5:, self.t_output].T.real - dUh_dt[5:, 0].real - d2Uf_dt2 = d2Uh_dt2[5:, self.t_output].T.real - d2Uh_dt2[5:, 0].real - - g_cg = np.zeros((len(self.t_output), 3)) - commands = np.zeros((len(self.t_output), self.solution.n_inputs)) + Ucg = np.concatenate((np.zeros((len(self.idx_output), 1)), Uh[:5, self.idx_output].T.real - Uh[:5, 0].real), axis=1) + dUcg_dt = np.concatenate((np.zeros((len(self.idx_output), 1)), + dUh_dt[:5, self.idx_output].T.real - dUh_dt[:5, 0].real), axis=1) + d2Ucg_dt2 = np.concatenate((np.zeros((len(self.idx_output), 1)), + d2Uh_dt2[:5, self.idx_output].T.real - d2Uh_dt2[:5, 0].real), axis=1) + Uf = Uh[5:, self.idx_output].T.real - Uh[5:, 0].real + dUf_dt = dUh_dt[5:, self.idx_output].T.real - dUh_dt[5:, 0].real + d2Uf_dt2 = d2Uh_dt2[5:, self.idx_output].T.real - d2Uh_dt2[5:, 0].real + + g_cg = np.zeros((len(self.idx_output), 3)) + commands = np.zeros((len(self.idx_output), self.solution.n_inputs)) # x, y, z, Phi, Theta, Psi, u, v, w, p, q, r in DIN 9300 body fixed system for flight physics X = np.concatenate((Ucg * np.array([-1., 1., -1., -1., 1., -1.]), @@ -68,7 +88,7 @@ def eval_equations(self): commands, ), axis=1) response = {'X': X, - 't': np.array([self.t[self.t_output]]).T, + 't': np.array([self.t[self.idx_output]]).T, 'Pk_aero': Pk_aero.T, 'Pk_gust': Pk_gust.T, 'Pk_unsteady': Pk_aero.T * 0.0, @@ -96,8 +116,8 @@ def setup_frequence_parameters(self): self.n_freqs += 1 # make even # sample spacing self.t = np.linspace(0.0, self.n_freqs * dt, self.n_freqs) - # indices of time samples to returned for post-processing - self.t_output = np.where(self.t <= self.simcase['t_final'])[0] + # indices of time samples to be returned for post-processing + self.idx_output = np.where(self.t <= self.simcase['t_final'])[0] # samples only from zero up to the Nyquist frequency self.freqs = np.linspace(0.0, self.fmax / 2.0, self.n_freqs // 2) # whole frequency space including negative frequencies @@ -107,11 +127,8 @@ def setup_frequence_parameters(self): self.positiv_fftfreqs = np.abs(fftfreqs[:self.n_freqs // 2 + 1]) self.positiv_fftomega = 2.0 * np.pi * self.positiv_fftfreqs - logging.info('Frequency domain solution with tfinal = {}x{} s, nfreq = {}, fmax={} Hz and df = {} Hz'.format( - t_factor, self.simcase['t_final'], self.n_freqs // 2, self.fmax / 2.0, self.fmax / self.n_freqs)) - if self.f2k(self.freqs.max()) > np.max(self.aero['k_red']): - logging.warning('Required reduced frequency = {:0.3} but AICs given only up to {:0.3}'.format( - self.f2k(self.freqs.max()), np.max(self.aero['k_red']))) + logging.info('Frequency domain solution with t_final = %s x %s s, n_freq = %s, f_max = %s Hz and df = %s Hz', + t_factor, self.simcase['t_final'], self.n_freqs // 2, self.fmax / 2.0, self.fmax / self.n_freqs) def mirror_fouriersamples_even(self, fouriersamples): mirrored_fourier = np.zeros((fouriersamples.shape[0], self.n_freqs), dtype='complex128') @@ -218,7 +235,8 @@ def calc_psd_vonKarman(self, freqs): if freqs[0] == 0.0: psd_karman[0] = 0.0 # Calculate the RMS value for cross-checking. Exclude first frequency with f=0.0 from the integral. - logging.info('RMS of PSD input (should approach 1.0): {:.4f}'.format(np.trapezoid(psd_karman[1:], freqs[1:]) ** 0.5)) + logging.info('RMS of PSD input (should approach 1.0): %.4f', np.trapezoid(psd_karman[1:], freqs[1:]) ** 0.5) + return psd_karman def calc_gust_excitation(self, freqs, t): @@ -227,7 +245,7 @@ def calc_gust_excitation(self, freqs, t): # during the simulation time. sigma = self.calc_sigma(self.n_freqs) u_sigma = self.u_sigma / sigma # turbulence gust intensity [m/s] - logging.info("Using RMS turbulence intensity u_sigma = {:.4f} m/s, sigma = {:.4f}.".format(u_sigma, sigma)) + logging.info("Using RMS turbulence intensity u_sigma = %.4f m/s, sigma = %.4f.", u_sigma, sigma) psd_karman = self.calc_psd_vonKarman(freqs) # Apply a scaling in the frequency domain to achieve the correct amplitude in the time domain. @@ -374,6 +392,16 @@ def calc_white_noise_excitation(self, freqs): class KMethod(GustExcitation): + def __init__(self, *args, **kwargs): + # Initialize additional attributes to avoid defining them outside __init__ + super().__init__(*args, **kwargs) + self.k_reds = None + self.A = None + self.B = None + self.freqs = None + self.damping = None + self.Qhh_interp = None + def eval_equations(self): self.setup_frequence_parameters() @@ -395,8 +423,8 @@ def setup_frequence_parameters(self): self.n_freqs = len(self.k_reds) if self.k_reds.max() > np.max(self.aero['k_red']): - logging.warning('Required reduced frequency = {:0.3} but AICs given only up to {:0.3}'.format( - self.k_reds.max(), np.max(self.aero['k_red']))) + logging.warning('Required reduced frequency = %0.3f but AICs given only up to %0.3f', + self.k_reds.max(), np.max(self.aero['k_red'])) def build_AIC_interpolators(self): Qhh = [] @@ -463,6 +491,7 @@ def calc_eigenvalues(self): class KEMethod(KMethod): + # No new attributes needed, thus no init def system(self, k_red): rho = self.atmo['rho'] @@ -532,6 +561,14 @@ class PKMethodSchwochow(KMethod): Modellierungsunsicherheiten am Flugzeug zur”, Dissertation, Universität Kassel, Kassel, 2012. """ + def __init__(self, *args, **kwargs): + # Initialize additional attributes to avoid defining them outside __init__ + super().__init__(*args, **kwargs) + self.n_modes_rbm = None + self.n_modes_f = None + self.states = None + self.Vvec = None + def setup_frequence_parameters(self): self.n_modes_rbm = 5 self.n_modes_f = self.model['mass'][self.trimcase['mass']]['n_modes'][()] @@ -542,7 +579,7 @@ def setup_frequence_parameters(self): self.states += ['Uf' + str(i_mode)] self.states += ["v'", "w'", "p'", "q'", "r'"] for i_mode in range(1, self.n_modes_f + 1): - self.states += ['$\\mathrm{{ \\dot Uf{} }}$'.format(str(i_mode))] # noqa: W605 + self.states += [f'$\\mathrm{{ \\dot Uf{i_mode} }}$'] # noqa: W605 self.Vvec = self.simcase['flutter_para']['Vtas'] @@ -555,9 +592,9 @@ def eval_equations(self): # Compute initial guess at k_red=0.0 and first flight speed self.Vtas = self.Vvec[0] eigenvalue, eigenvector = linalg.eig(self.system(k_red=0.0)) - bandbreite = eigenvalue.__abs__().max() - eigenvalue.__abs__().min() + bandbreite = np.abs(eigenvalue).max() - np.abs(eigenvalue).min() # No zero eigenvalues - idx_pos = np.where(eigenvalue.__abs__() / bandbreite >= 1e-3)[0] + idx_pos = np.where(np.abs(eigenvalue) / bandbreite >= 1e-3)[0] # Sort initial results by eigenvalue idx_sort = np.argsort(np.abs(eigenvalue.imag[idx_pos])) eigenvalues0 = eigenvalue[idx_pos][idx_sort] @@ -571,7 +608,7 @@ def eval_equations(self): Vtas = [] # Loop over modes for i_mode in range(len(eigenvalues0)): - logging.debug('Mode {}'.format(i_mode + 1)) + logging.debug('Mode %d', i_mode + 1) eigenvalues_per_mode = [] eigenvectors_per_mode = [] k_old = copy.deepcopy(k0[i_mode]) @@ -593,6 +630,9 @@ def eval_equations(self): elif self.simcase['flutter_para']['method'] in ['pk_rodden']: # Allow only positive reduced frequencies in the implementation following Rodden. k_now = np.abs(eigenvalues_new[i_mode].imag) * self.macgrid['c_ref'] / 2.0 / self.Vtas + else: + # We should never reach this point + k_now = 0.0 # Use relaxation for improved convergence, which helps in some cases to avoid oscillations of the # iterative solution. k_new = k_old + 0.8 * (k_now - k_old) @@ -602,8 +642,8 @@ def eval_equations(self): # If no convergence is achieved, stop and issue a warning. Typically, the iteration converges in less than # ten loops, so 50 should be more than enough and prevents excessive calculation times. if n_iter > 50: - logging.warning('No convergence for mode {} at Vtas={:.2f} with k_red={:.5f} and e={:.5f}'.format( - i_mode + 1, V, k_new, e)) + logging.warning('No convergence for mode %d at Vtas=%.2f with k_red=%.5f and e=%.5f', + i_mode + 1, V, k_new, e) break eigenvalues_old = eigenvalues_new eigenvectors_old = eigenvectors_new @@ -632,7 +672,7 @@ def calc_eigenvalues(self, A, eigenvalues_old, eigenvectors_old): # To match the modes with the previous step, use a correlation cirterion as specified in the JCL. if 'tracking' not in self.simcase['flutter_para']: # Set a default. - tracking_method = 'MAC' + tracking_method = 'MAC*PCC' else: tracking_method = self.simcase['flutter_para']['tracking'] # Calculate the correlation bewteen the old and current modes. @@ -645,6 +685,9 @@ def calc_eigenvalues(self, A, eigenvalues_old, eigenvectors_old): elif tracking_method == 'MAC*HDM': # Combining MAC and hyperboloic distance metric (HDM) for improved handling of complex conjugate pairs. correlation = fem_helper.calc_MAC(eigenvectors_old, eigenvector) * fem_helper.calc_HDM(eigenvalues_old, eigenvalue) + else: + # Because a default was selected above, we should never reach this point. + correlation = None # Based on the correlation matrix, find the best match and apply to the modes. idx_pos = self.get_best_match(correlation) eigenvalues = eigenvalue[idx_pos] @@ -729,8 +772,8 @@ def system(self, k_red): Qhh = self.Qhh_interp(k_red) Mhh_inv = np.linalg.inv(self.Mhh) - upper_part = np.concatenate((np.zeros((self.n_modes, self.n_modes), dtype='complex128'), - np.eye(self.n_modes, dtype='complex128')), axis=1) + upper_part = np.concatenate((np.zeros((self.n_modes, self.n_modes)), + np.eye(self.n_modes)), axis=1) lower_part = np.concatenate((-Mhh_inv.dot(self.Khh - rho / 2 * self.Vtas ** 2.0 * Qhh.real), -Mhh_inv.dot(self.Dhh - rho / 4 * self.Vtas * self.macgrid['c_ref'] / k_red * Qhh.imag)), axis=1) diff --git a/loadskernel/equations/state_space.py b/loadskernel/equations/mona_state_space.py similarity index 98% rename from loadskernel/equations/state_space.py rename to loadskernel/equations/mona_state_space.py index 2e457a0e..ce633f83 100644 --- a/loadskernel/equations/state_space.py +++ b/loadskernel/equations/mona_state_space.py @@ -4,7 +4,7 @@ from scipy import linalg -from loadskernel.equations.frequency_domain import PKMethodSchwochow +from loadskernel.equations.mona_frequency_domain import PKMethodSchwochow class StateSpaceAnalysis(PKMethodSchwochow): diff --git a/loadskernel/equations/mona_time_domain.py b/loadskernel/equations/mona_time_domain.py new file mode 100644 index 00000000..1dc2c40e --- /dev/null +++ b/loadskernel/equations/mona_time_domain.py @@ -0,0 +1,604 @@ +import logging +import numpy as np +from scipy import linalg + +from loadskernel.equations.common import Common, ConvergenceError +from loadskernel.solution_tools import gravitation_on_earth + + +class Steady(Common): + + def equations(self, X, t, modus): + self.counter += 1 + # recover states + Tgeo2body, Tbody2geo = self.geo2body(X) + dUcg_dt, Uf, dUf_dt = self.recover_states(X) + Vtas, q_dyn = self.recover_Vtas(X) + onflow = self.recover_onflow(X) + alpha, beta, gamma = self.windsensor(X, Vtas, Uf, dUf_dt) + Ux2 = self.get_Ux2(X) + + # aerodynamics + Pk_rbm, wj_rbm = self.rbm(onflow, q_dyn, Vtas) + Pk_cam, wj_cam = self.camber_twist(q_dyn) + Pk_cs, wj_cs = self.cs(X, Ux2, q_dyn) + Pk_f, wj_f = self.flexible(Uf, dUf_dt, onflow, q_dyn, Vtas) + Pk_gust, wj_gust = self.gust(X, q_dyn) + + wj = wj_rbm + wj_cam + wj_cs + wj_f + wj_gust + Pk_idrag = self.idrag(wj, q_dyn) + Pk_unsteady = Pk_rbm * 0.0 + + Pextra, Pb_ext, Pf_ext = self.engine(X, Vtas, q_dyn, Uf, dUf_dt, t) + + # correction coefficients + Pb_corr = self.correctioon_coefficients(alpha, beta, q_dyn) + + # summation of forces + Pk_aero = Pk_rbm + Pk_cam + Pk_cs + Pk_f + Pk_gust + Pk_idrag + Pk_unsteady + Pmac = np.dot(self.Dkx1.T, Pk_aero) + Pb = np.dot(self.PHImac_cg.T, Pmac) + Pb_corr + Pb_ext + + g_cg = gravitation_on_earth(self.PHInorm_cg, Tgeo2body) + + # EoM + d2Ucg_dt2, Nxyz = self.rigid_EoM(dUcg_dt, Pb, g_cg, modus) + Pf = np.dot(self.PHIkf.T, Pk_aero) + self.Mfcg.dot(np.hstack((d2Ucg_dt2[0:3] - g_cg, d2Ucg_dt2[3:6]))) + Pf_ext + d2Uf_dt2 = self.flexible_EoM(dUf_dt, Uf, Pf) + + # CS derivatives + dcommand = self.get_command_derivatives(t, X, Vtas, gamma, alpha, beta, Nxyz, np.dot(Tbody2geo, X[6:12])[0:3]) + + # output + Y = np.hstack((np.dot(Tbody2geo, X[6:12]), + np.dot(self.PHIcg_norm, d2Ucg_dt2), + dUf_dt, + d2Uf_dt2, + dcommand, + Nxyz[2], + Vtas, + beta, + )) + + if modus in ['trim', 'sim']: + return Y + elif modus in ['trim_full_output']: + response = {'X': X, + 'Y': Y, + 't': np.array([t]), + 'Pk_rbm': Pk_rbm, + 'Pk_cam': Pk_cam, + 'Pk_aero': Pk_aero, + 'Pk_cs': Pk_cs, + 'Pk_f': Pk_f, + 'Pk_gust': Pk_gust, + 'Pk_unsteady': Pk_unsteady, + 'Pk_idrag': Pk_idrag, + 'q_dyn': np.array([q_dyn]), + 'Pb': Pb, + 'Pmac': Pmac, + 'Pf': Pf, + 'alpha': np.array([alpha]), + 'beta': np.array([beta]), + 'Ux2': Ux2, + 'dUcg_dt': dUcg_dt, + 'd2Ucg_dt2': d2Ucg_dt2, + 'Uf': Uf, + 'dUf_dt': dUf_dt, + 'd2Uf_dt2': d2Uf_dt2, + 'Nxyz': Nxyz, + 'g_cg': g_cg, + 'Pextra': Pextra, + } + return response + elif modus in ['sim_full_output']: + # For time domain simulations, typically not all results are required. To reduce the amount of data while + # maintaining compatibility with the trim, empty arrays are used. + response = {'X': X, + 'Y': Y, + 't': np.array([t]), + 'Pk_aero': Pk_aero, + 'q_dyn': np.array([q_dyn]), + 'Pb': Pb, + 'Pmac': Pmac, + 'alpha': np.array([alpha]), + 'beta': np.array([beta]), + 'Ux2': Ux2, + 'dUcg_dt': dUcg_dt, + 'd2Ucg_dt2': d2Ucg_dt2, + 'Uf': Uf, + 'dUf_dt': dUf_dt, + 'd2Uf_dt2': d2Uf_dt2, + 'Nxyz': Nxyz, + 'g_cg': g_cg, + 'Pextra': Pextra, + } + return response + + def eval_equations(self, X_free, time, modus='trim_full_output'): + # this is a wrapper for the model equations 'eqn_basic' + if modus in ['trim', 'trim_full_output']: + # get inputs from trimcond and apply inputs from fsolve + X = np.array(self.trimcond_X[:, 2], dtype='float') + X[np.where((self.trimcond_X[:, 1] == 'free'))[0]] = X_free + elif modus in ['sim', 'sim_full_output']: + X = X_free + + # evaluate model equations + if modus == 'trim': + Y = self.equations(X, time, 'trim') + # get the current values from Y and substract tamlab.figure() + # fsolve only finds the roots; Y = 0 + Y_target_ist = Y[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] + Y_target_soll = np.array(self.trimcond_Y[:, 2], dtype='float')[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] + out = Y_target_ist - Y_target_soll + return out + + elif modus == 'sim': + Y = self.equations(X, time, 'sim') + return Y[self.solution.idx_state_derivatives + self.solution.idx_input_derivatives] + + elif modus == 'sim_full_output': + response = self.equations(X, time, 'sim_full_output') + return response + + elif modus == 'trim_full_output': + response = self.equations(X, time, 'trim_full_output') + return response + + def eval_equations_iteratively(self, X_free, time, modus='trim_full_output'): + # this is a wrapper for the model equations + + # get inputs from trimcond and apply inputs from fsolve + X = np.array(self.trimcond_X[:, 2], dtype='float') + X[np.where((self.trimcond_X[:, 1] == 'free'))[0]] = X_free + logging.debug('X_free: {}'.format(X_free)) + converged = False + inner_loops = 0 + while not converged: + inner_loops += 1 + response = self.equations(X, time, 'trim_full_output') + logging.info('Inner iteration {}, calculate structural deformations.'.format(self.counter)) + Uf_new = linalg.solve(self.Kff, response['Pf']) + + # Add a relaxation factor between each loop to reduce overshoot and/or oscillations. + # In case the solution hasn't converged in a reasonable time (say 8 loops), increase the relaxation. + # A low relaxation factor is slower but more robust, f_relax = 1.0 implies no relaxation. + # After 16 inner loops, decrease the convergence criterion epsilon. + if inner_loops < 8: + f_relax = 0.8 + else: + f_relax = 0.5 + if inner_loops < 16: + epsilon = self.jcl.general['b_ref'] * 1.0e-5 + else: + epsilon = self.jcl.general['b_ref'] * 1.0e-4 + logging.info(' - Relaxation factor: {}'.format(f_relax)) + logging.info(' - Epsilon: {:0.6g}'.format(epsilon)) + + # recover Uf_old from last step and blend with Uf_now + Uf_old = [self.trimcond_X[np.where((self.trimcond_X[:, 0] == 'Uf' + str(i_mode)))[0][0], 2] for i_mode in + range(1, self.n_modes + 1)] + Uf_old = np.array(Uf_old, dtype='float') + Uf_new = Uf_new * f_relax + Uf_old * (1.0 - f_relax) + + # set new values for Uf in trimcond for next loop and store in response + for i_mode in range(self.n_modes): + self.trimcond_X[np.where((self.trimcond_X[:, 0] == 'Uf' + str(i_mode + 1)))[0][0], 2] = '{:g}'.format( + Uf_new[i_mode]) + response['X'][12 + i_mode] = Uf_new[i_mode] + + # convergence parameter for iterative evaluation + Ug_f_body = np.dot(self.PHIf_strc.T, Uf_new.T).T + defo_new = Ug_f_body[self.strcgrid['set'][:, :3]].max() # Groesste Verformung, meistens Fluegelspitze + ddefo = defo_new - self.defo_old + self.defo_old = np.copy(defo_new) + if np.abs(ddefo) < epsilon: + converged = True + logging.info(' - Max. deformation: {:0.6g}, delta: {:0.6g} smaller than epsilon, converged.'.format( + defo_new, ddefo)) + else: + logging.info(' - Max. deformation: {:0.6g}, delta: {:0.6g}'.format(defo_new, ddefo)) + if inner_loops > 20: + raise ConvergenceError('No convergence of structural deformation achieved after {} inner loops. \ + Check convergence of CFD solution and/or convergence criterion "delta".'.format(inner_loops)) + # get the current values from Y and substract tamlab.figure() + # fsolve only finds the roots; Y = 0 + Y_target_ist = response['Y'][np.where((self.trimcond_Y[:, 1] == 'target'))[0]] + Y_target_soll = np.array(self.trimcond_Y[:, 2], dtype='float')[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] + out = Y_target_ist - Y_target_soll + + if modus in ['trim']: + return out + elif modus == 'trim_full_output': + return response + + +class Unsteady(Common): + + def equations(self, X, t, modus): + self.counter += 1 + # recover states + Tgeo2body, Tbody2geo = self.geo2body(X) + dUcg_dt, Uf, dUf_dt = self.recover_states(X) + Vtas, q_dyn = self.recover_Vtas(self.X0) + onflow = self.recover_onflow(X) + alpha, beta, gamma = self.windsensor(X, Vtas, Uf, dUf_dt) + Ux2 = self.get_Ux2(X) + + # aerodynamics + Pk_rbm, wj_rbm = self.rbm(onflow, q_dyn, Vtas) + Pk_cam, wj_cam = self.camber_twist(q_dyn) + Pk_cs, wj_cs = self.cs(X, Ux2, q_dyn) + Pk_f, wj_f = self.flexible(Uf, dUf_dt, onflow, q_dyn, Vtas) + Pk_gust, wj_gust = self.gust(X, q_dyn) + + wj = wj_rbm + wj_cam + wj_cs + wj_f + wj_gust + Pk_idrag = self.idrag(wj, q_dyn) + Pk_unsteady, dlag_states_dt = self.unsteady(X, t, wj, Uf, dUf_dt, onflow, q_dyn, Vtas) + + Pextra, Pb_ext, Pf_ext = self.engine(X, Vtas, q_dyn, Uf, dUf_dt, t) + + # correction coefficients + Pb_corr = self.correctioon_coefficients(alpha, beta, q_dyn) + + # summation of forces + Pk_aero = Pk_rbm + Pk_cam + Pk_cs + Pk_f + Pk_gust + Pk_idrag + Pk_unsteady + Pmac = np.dot(self.Dkx1.T, Pk_aero) + Pb = np.dot(self.PHImac_cg.T, Pmac) + Pb_corr + Pb_ext + + g_cg = gravitation_on_earth(self.PHInorm_cg, Tgeo2body) + + # EoM + d2Ucg_dt2, Nxyz = self.rigid_EoM(dUcg_dt, Pb, g_cg, modus) + Pf = np.dot(self.PHIkf.T, Pk_aero) + self.Mfcg.dot(np.hstack((d2Ucg_dt2[0:3] - g_cg, d2Ucg_dt2[3:6]))) + Pf_ext + d2Uf_dt2 = self.flexible_EoM(dUf_dt, Uf, Pf) + + # CS derivatives + dcommand = self.get_command_derivatives(t, X, Vtas, gamma, alpha, beta, Nxyz, np.dot(Tbody2geo, X[6:12])[0:3]) + + # output + Y = np.hstack((np.dot(Tbody2geo, X[6:12]), + np.dot(self.PHIcg_norm, d2Ucg_dt2), + dUf_dt, + d2Uf_dt2, + dcommand, + dlag_states_dt, + Nxyz[2], + Vtas, + beta, + )) + + if modus in ['trim', 'sim']: + return Y + elif modus in ['trim_full_output']: + response = {'X': X, + 'Y': Y, + 't': np.array([t]), + 'Pk_rbm': Pk_rbm, + 'Pk_cam': Pk_cam, + 'Pk_aero': Pk_aero, + 'Pk_cs': Pk_cs, + 'Pk_f': Pk_f, + 'Pk_gust': Pk_gust, + 'Pk_unsteady': Pk_unsteady, + 'Pk_idrag': Pk_idrag, + 'q_dyn': np.array([q_dyn]), + 'Pb': Pb, + 'Pmac': Pmac, + 'Pf': Pf, + 'alpha': np.array([alpha]), + 'beta': np.array([beta]), + 'Ux2': Ux2, + 'dUcg_dt': dUcg_dt, + 'd2Ucg_dt2': d2Ucg_dt2, + 'Uf': Uf, + 'dUf_dt': dUf_dt, + 'd2Uf_dt2': d2Uf_dt2, + 'Nxyz': Nxyz, + 'g_cg': g_cg, + 'Pextra': Pextra, + } + return response + elif modus in ['sim_full_output']: + # For time domain simulations, typically not all results are required. To reduce the amount of data while + # maintaining compatibility with the trim, empty arrays are used. + response = {'X': X, + 'Y': Y, + 't': np.array([t]), + 'Pk_aero': Pk_aero, + 'Pk_gust': Pk_gust, + 'Pk_unsteady': Pk_unsteady, + 'q_dyn': np.array([q_dyn]), + 'Pb': Pb, + 'Pmac': Pmac, + 'alpha': np.array([alpha]), + 'beta': np.array([beta]), + 'Ux2': Ux2, + 'dUcg_dt': dUcg_dt, + 'd2Ucg_dt2': d2Ucg_dt2, + 'Uf': Uf, + 'dUf_dt': dUf_dt, + 'd2Uf_dt2': d2Uf_dt2, + 'Nxyz': Nxyz, + 'g_cg': g_cg, + 'Pextra': Pextra, + } + return response + + def eval_equations(self, X_free, time, modus='trim_full_output'): + if modus in ['sim', 'sim_full_output']: + X = X_free + + # evaluate model equations + if modus == 'sim': + Y = self.equations(X, time, 'sim') + return Y[self.solution.idx_state_derivatives + self.solution.idx_input_derivatives + + self.solution.idx_lag_derivatives] + + elif modus == 'sim_full_output': + response = self.equations(X, time, 'sim_full_output') + return response + + +class NonlinSteady(Steady): + + def equations(self, X, t, modus): + self.counter += 1 + # recover states + Tgeo2body, Tbody2geo = self.geo2body(X) + dUcg_dt, Uf, dUf_dt = self.recover_states(X) + Vtas, q_dyn = self.recover_Vtas(X) + onflow = self.recover_onflow(X) + alpha, beta, gamma = self.windsensor(X, Vtas, Uf, dUf_dt) + Ux2 = self.get_Ux2(X) + + # aerodynamics + Pk_rbm, wj_rbm = self.rbm_nonlin(onflow) + Pk_cam, wj_cam = self.camber_twist_nonlin(onflow, Vtas) + Pk_cs, wj_cs = self.cs_nonlin(onflow, X, Ux2, Vtas) + Pk_f, wj_f = self.flexible_nonlin(onflow, Uf, dUf_dt, Vtas) + Pk_gust = Pk_rbm * 0.0 + + wj = (wj_rbm + wj_cam + wj_cs + wj_f) / Vtas + Pk_idrag = self.idrag(wj, q_dyn) + + Pk_unsteady = Pk_rbm * 0.0 + + Pextra, Pb_ext, Pf_ext = self.engine(X, Vtas, q_dyn, Uf, dUf_dt, t) + + # correction coefficients + Pb_corr = self.correctioon_coefficients(alpha, beta, q_dyn) + Pmac_vdrag = self.vdrag(alpha, q_dyn) + + # summation of forces + Pk_aero = Pk_rbm + Pk_cam + Pk_cs + Pk_f + Pk_gust + Pk_idrag + Pk_unsteady + Pmac = np.dot(self.Dkx1.T, Pk_aero) + Pmac_vdrag + Pb = np.dot(self.PHImac_cg.T, Pmac) + Pb_corr + Pb_ext + + g_cg = gravitation_on_earth(self.PHInorm_cg, Tgeo2body) + + # EoM + d2Ucg_dt2, Nxyz = self.rigid_EoM(dUcg_dt, Pb, g_cg, modus) + Pf = np.dot(self.PHIkf.T, Pk_aero) + self.Mfcg.dot(np.hstack((d2Ucg_dt2[0:3] - g_cg, d2Ucg_dt2[3:6]))) + Pf_ext + d2Uf_dt2 = self.flexible_EoM(dUf_dt, Uf, Pf) + + # CS derivatives + dcommand = self.get_command_derivatives(t, X, Vtas, gamma, alpha, beta, Nxyz, np.dot(Tbody2geo, X[6:12])[0:3]) + + # output + Y = np.hstack((np.dot(Tbody2geo, X[6:12]), + np.dot(self.PHIcg_norm, d2Ucg_dt2), + dUf_dt, + d2Uf_dt2, + dcommand, + Nxyz[2], + Vtas, + beta, + )) + + if modus in ['trim', 'sim']: + return Y + elif modus in ['trim_full_output']: + response = {'X': X, + 'Y': Y, + 't': np.array([t]), + 'Pk_rbm': Pk_rbm, + 'Pk_cam': Pk_cam, + 'Pk_aero': Pk_aero, + 'Pk_cs': Pk_cs, + 'Pk_f': Pk_f, + 'Pk_gust': Pk_gust, + 'Pk_unsteady': Pk_unsteady, + 'Pk_idrag': Pk_idrag, + 'Pmac_vdrag': Pmac_vdrag, + 'q_dyn': np.array([q_dyn]), + 'Pb': Pb, + 'Pmac': Pmac, + 'Pf': Pf, + 'alpha': np.array([alpha]), + 'beta': np.array([beta]), + 'Ux2': Ux2, + 'dUcg_dt': dUcg_dt, + 'd2Ucg_dt2': d2Ucg_dt2, + 'Uf': Uf, + 'dUf_dt': dUf_dt, + 'd2Uf_dt2': d2Uf_dt2, + 'Nxyz': Nxyz, + 'g_cg': g_cg, + 'Pextra': Pextra, + } + return response + elif modus in ['sim_full_output']: + # For time domain simulations, typically not all results are required. To reduce the amount of data while + # maintaining compatibility with the trim, empty arrays are used. + response = {'X': X, + 'Y': Y, + 't': np.array([t]), + 'Pk_aero': Pk_aero, + 'q_dyn': np.array([q_dyn]), + 'Pb': Pb, + 'Pmac': Pmac, + 'alpha': np.array([alpha]), + 'beta': np.array([beta]), + 'Ux2': Ux2, + 'dUcg_dt': dUcg_dt, + 'd2Ucg_dt2': d2Ucg_dt2, + 'Uf': Uf, + 'dUf_dt': dUf_dt, + 'd2Uf_dt2': d2Uf_dt2, + 'Nxyz': Nxyz, + 'g_cg': g_cg, + 'Pextra': Pextra, + } + return response + + +class Landing(Common): + + def equations(self, X, t, modus): + self.counter += 1 + # recover states + Tgeo2body, Tbody2geo = self.geo2body(X) + dUcg_dt, Uf, dUf_dt = self.recover_states(X) + Vtas, q_dyn = self.recover_Vtas(X) + onflow = self.recover_onflow(X) + alpha, beta, gamma = self.windsensor(X, Vtas, Uf, dUf_dt) + Ux2 = self.get_Ux2(X) + + # aerodynamics + Pk_rbm, wj_rbm = self.rbm(onflow, q_dyn, Vtas) + Pk_cam, wj_cam = self.camber_twist(q_dyn) + Pk_cs, wj_cs = self.cs(X, Ux2, q_dyn) + Pk_f, wj_f = self.flexible(Uf, dUf_dt, onflow, q_dyn, Vtas) + Pk_gust, wj_gust = self.gust(X, q_dyn) + + wj = wj_rbm + wj_cam + wj_cs + wj_f + wj_gust + Pk_idrag = self.idrag(wj, q_dyn) + Pk_unsteady = Pk_rbm * 0.0 + + # correction coefficients + Pb_corr = self.correctioon_coefficients(alpha, beta, q_dyn) + + # landing gear + Pextra, _, dp2, ddp2, F1, F2 = self.landinggear(X, Tbody2geo) + + # summation of forces + Pk_aero = Pk_rbm + Pk_cam + Pk_cs + Pk_f + Pk_gust + Pk_idrag + Pk_unsteady + Pmac = np.dot(self.Dkx1.T, Pk_aero) + Pb = np.dot(self.PHImac_cg.T, Pmac) + Pb_corr + np.dot(self.PHIextra_cg.T, Pextra) + + g_cg = gravitation_on_earth(self.PHInorm_cg, Tgeo2body) + + # EoM + d2Ucg_dt2, Nxyz = self.rigid_EoM(dUcg_dt, Pb, g_cg, modus) + Pf = np.dot(self.PHIkf.T, Pk_aero) + self.Mfcg.dot( + np.hstack((d2Ucg_dt2[0:3] - g_cg, d2Ucg_dt2[3:6]))) + np.dot(self.PHIf_extra, Pextra) # viel schneller! + d2Uf_dt2 = self.flexible_EoM(dUf_dt, Uf, Pf) + + # CS derivatives + dcommand = self.get_command_derivatives(t, X, Vtas, gamma, alpha, beta, Nxyz, np.dot(Tbody2geo, X[6:12])[0:3]) + + # output + if modus in ['trim', 'trim_full_output']: + Y = np.hstack((np.dot(Tbody2geo, X[6:12]), + np.dot(self.PHIcg_norm, d2Ucg_dt2), + dUf_dt, + d2Uf_dt2, + dcommand, + Nxyz[2], + Vtas, + beta, + )) + elif modus in ['sim', 'sim_full_output']: + Y = np.hstack((np.dot(Tbody2geo, X[6:12]), + np.dot(self.PHIcg_norm, d2Ucg_dt2), + dUf_dt, + d2Uf_dt2, + dcommand, + np.hstack((dp2, ddp2)), + Nxyz[2], + Vtas, + beta, + )) + if modus in ['trim', 'sim']: + return Y + elif modus in ['trim_full_output', 'sim_full_output']: + # calculate translations, velocities and accelerations of some additional points + # (might also be used for sensors in a closed-loop system + # position LG attachment point over ground + p1 = -self.cggrid['offset'][:, 2] + self.extragrid['offset'][:, 2] \ + + self.PHIextra_cg.dot(np.dot(self.PHInorm_cg, X[0:6]))[self.extragrid['set'][:, 2]] \ + + self.PHIf_extra.T.dot(X[12:12 + self.n_modes])[self.extragrid['set'][:, 2]] + # velocity LG attachment point + dp1 = self.PHIextra_cg.dot(np.dot(self.PHInorm_cg, np.dot(Tbody2geo, X[6:12])))[self.extragrid['set'][:, 2]] \ + + self.PHIf_extra.T.dot(X[12 + self.n_modes:12 + self.n_modes * 2])[self.extragrid['set'][:, 2]] + # acceleration LG attachment point + ddp1 = self.PHIextra_cg.dot(np.dot(self.PHInorm_cg, np.dot(Tbody2geo, Y[6:12])))[self.extragrid['set'][:, 2]] \ + + self.PHIf_extra.T.dot(Y[12 + self.n_modes:12 + self.n_modes * 2])[self.extragrid['set'][:, 2]] + + response = {'X': X, + 'Y': Y, + 't': np.array([t]), + 'Pk_rbm': Pk_rbm, + 'Pk_cam': Pk_cam, + 'Pk_aero': Pk_aero, + 'Pk_cs': Pk_cs, + 'Pk_f': Pk_f, + 'Pk_gust': Pk_gust, + 'Pk_unsteady': Pk_unsteady, + 'Pk_idrag': Pk_idrag, + 'q_dyn': np.array([q_dyn]), + 'Pb': Pb, + 'Pmac': Pmac, + 'Pf': Pf, + 'alpha': np.array([alpha]), + 'beta': np.array([beta]), + 'Ux2': Ux2, + 'dUcg_dt': dUcg_dt, + 'd2Ucg_dt2': d2Ucg_dt2, + 'Uf': Uf, + 'dUf_dt': dUf_dt, + 'd2Uf_dt2': d2Uf_dt2, + 'Nxyz': Nxyz, + 'g_cg': g_cg, + 'Pextra': Pextra, + 'p1': p1, + 'dp1': dp1, + 'ddp1': ddp1, + 'F1': F1, + 'F2': F2, + } + return response + + def eval_equations(self, X_free, time, modus='trim_full_output'): + + if modus in ['trim', 'trim_full_output']: + # get inputs from trimcond and apply inputs from fsolve + X = np.array(self.trimcond_X[:, 2], dtype='float') + X[np.where((self.trimcond_X[:, 1] == 'free'))[0]] = X_free + elif modus in ['sim', 'sim_full_output']: + X = X_free + + # evaluate model equations + if modus == 'trim': + Y = self.equations(X, time, 'trim') + # get the current values from Y and substract tamlab.figure() + # fsolve only finds the roots; Y = 0 + Y_target_ist = Y[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] + Y_target_soll = np.array(self.trimcond_Y[:, 2], dtype='float')[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] + out = Y_target_ist - Y_target_soll + return out + + elif modus == 'sim': + Y = self.equations(X, time, 'sim') + return Y[self.solution.idx_state_derivatives + self.solution.idx_input_derivatives + + self.solution.idx_lg_derivatives] # Nz ist eine Rechengroesse und keine Simulationsgroesse! + + elif modus == 'sim_full_output': + response = self.equations(X, time, 'sim_full_output') + return response + + elif modus == 'trim_full_output': + response = self.equations(X, time, 'trim_full_output') + return response diff --git a/loadskernel/equations/nonlin_steady.py b/loadskernel/equations/nonlin_steady.py deleted file mode 100644 index 707ec5f5..00000000 --- a/loadskernel/equations/nonlin_steady.py +++ /dev/null @@ -1,117 +0,0 @@ -import numpy as np - -from loadskernel.equations.steady import Steady -from loadskernel.solution_tools import gravitation_on_earth - - -class NonlinSteady(Steady): - - def equations(self, X, t, modus): - self.counter += 1 - # recover states - Tgeo2body, Tbody2geo = self.geo2body(X) - dUcg_dt, Uf, dUf_dt = self.recover_states(X) - Vtas, q_dyn = self.recover_Vtas(X) - onflow = self.recover_onflow(X) - alpha, beta, gamma = self.windsensor(X, Vtas, Uf, dUf_dt) - Ux2 = self.get_Ux2(X) - - # aerodynamics - Pk_rbm, wj_rbm = self.rbm_nonlin(onflow) - Pk_cam, wj_cam = self.camber_twist_nonlin(onflow, Vtas) - Pk_cs, wj_cs = self.cs_nonlin(onflow, X, Ux2, Vtas) - Pk_f, wj_f = self.flexible_nonlin(onflow, Uf, dUf_dt, Vtas) - Pk_gust = Pk_rbm * 0.0 - - wj = (wj_rbm + wj_cam + wj_cs + wj_f) / Vtas - Pk_idrag = self.idrag(wj, q_dyn) - - Pk_unsteady = Pk_rbm * 0.0 - - Pextra, Pb_ext, Pf_ext = self.engine(X, Vtas, q_dyn, Uf, dUf_dt, t) - - # correction coefficients - Pb_corr = self.correctioon_coefficients(alpha, beta, q_dyn) - Pmac_vdrag = self.vdrag(alpha, q_dyn) - - # summation of forces - Pk_aero = Pk_rbm + Pk_cam + Pk_cs + Pk_f + Pk_gust + Pk_idrag + Pk_unsteady - Pmac = np.dot(self.Dkx1.T, Pk_aero) + Pmac_vdrag - Pb = np.dot(self.PHImac_cg.T, Pmac) + Pb_corr + Pb_ext - - g_cg = gravitation_on_earth(self.PHInorm_cg, Tgeo2body) - - # EoM - d2Ucg_dt2, Nxyz = self.rigid_EoM(dUcg_dt, Pb, g_cg, modus) - Pf = np.dot(self.PHIkf.T, Pk_aero) + self.Mfcg.dot(np.hstack((d2Ucg_dt2[0:3] - g_cg, d2Ucg_dt2[3:6]))) + Pf_ext - d2Uf_dt2 = self.flexible_EoM(dUf_dt, Uf, Pf) - - # CS derivatives - dcommand = self.get_command_derivatives(t, X, Vtas, gamma, alpha, beta, Nxyz, np.dot(Tbody2geo, X[6:12])[0:3]) - - # output - Y = np.hstack((np.dot(Tbody2geo, X[6:12]), - np.dot(self.PHIcg_norm, d2Ucg_dt2), - dUf_dt, - d2Uf_dt2, - dcommand, - Nxyz[2], - Vtas, - beta, - )) - - if modus in ['trim', 'sim']: - return Y - elif modus in ['trim_full_output']: - response = {'X': X, - 'Y': Y, - 't': np.array([t]), - 'Pk_rbm': Pk_rbm, - 'Pk_cam': Pk_cam, - 'Pk_aero': Pk_aero, - 'Pk_cs': Pk_cs, - 'Pk_f': Pk_f, - 'Pk_gust': Pk_gust, - 'Pk_unsteady': Pk_unsteady, - 'Pk_idrag': Pk_idrag, - 'Pmac_vdrag': Pmac_vdrag, - 'q_dyn': np.array([q_dyn]), - 'Pb': Pb, - 'Pmac': Pmac, - 'Pf': Pf, - 'alpha': np.array([alpha]), - 'beta': np.array([beta]), - 'Ux2': Ux2, - 'dUcg_dt': dUcg_dt, - 'd2Ucg_dt2': d2Ucg_dt2, - 'Uf': Uf, - 'dUf_dt': dUf_dt, - 'd2Uf_dt2': d2Uf_dt2, - 'Nxyz': Nxyz, - 'g_cg': g_cg, - 'Pextra': Pextra, - } - return response - elif modus in ['sim_full_output']: - # For time domain simulations, typically not all results are required. To reduce the amount of data while - # maintaining compatibility with the trim, empty arrays are used. - response = {'X': X, - 'Y': Y, - 't': np.array([t]), - 'Pk_aero': Pk_aero, - 'q_dyn': np.array([q_dyn]), - 'Pb': Pb, - 'Pmac': Pmac, - 'alpha': np.array([alpha]), - 'beta': np.array([beta]), - 'Ux2': Ux2, - 'dUcg_dt': dUcg_dt, - 'd2Ucg_dt2': d2Ucg_dt2, - 'Uf': Uf, - 'dUf_dt': dUf_dt, - 'd2Uf_dt2': d2Uf_dt2, - 'Nxyz': Nxyz, - 'g_cg': g_cg, - 'Pextra': Pextra, - } - return response diff --git a/loadskernel/equations/steady.py b/loadskernel/equations/steady.py deleted file mode 100644 index f04ef860..00000000 --- a/loadskernel/equations/steady.py +++ /dev/null @@ -1,214 +0,0 @@ -import logging -import numpy as np -from scipy import linalg - -from loadskernel.equations.common import Common, ConvergenceError -from loadskernel.solution_tools import gravitation_on_earth - - -class Steady(Common): - - def equations(self, X, t, modus): - self.counter += 1 - # recover states - Tgeo2body, Tbody2geo = self.geo2body(X) - dUcg_dt, Uf, dUf_dt = self.recover_states(X) - Vtas, q_dyn = self.recover_Vtas(X) - onflow = self.recover_onflow(X) - alpha, beta, gamma = self.windsensor(X, Vtas, Uf, dUf_dt) - Ux2 = self.get_Ux2(X) - - # aerodynamics - Pk_rbm, wj_rbm = self.rbm(onflow, q_dyn, Vtas) - Pk_cam, wj_cam = self.camber_twist(q_dyn) - Pk_cs, wj_cs = self.cs(X, Ux2, q_dyn) - Pk_f, wj_f = self.flexible(Uf, dUf_dt, onflow, q_dyn, Vtas) - Pk_gust, wj_gust = self.gust(X, q_dyn) - - wj = wj_rbm + wj_cam + wj_cs + wj_f + wj_gust - Pk_idrag = self.idrag(wj, q_dyn) - Pk_unsteady = Pk_rbm * 0.0 - - Pextra, Pb_ext, Pf_ext = self.engine(X, Vtas, q_dyn, Uf, dUf_dt, t) - - # correction coefficients - Pb_corr = self.correctioon_coefficients(alpha, beta, q_dyn) - - # summation of forces - Pk_aero = Pk_rbm + Pk_cam + Pk_cs + Pk_f + Pk_gust + Pk_idrag + Pk_unsteady - Pmac = np.dot(self.Dkx1.T, Pk_aero) - Pb = np.dot(self.PHImac_cg.T, Pmac) + Pb_corr + Pb_ext - - g_cg = gravitation_on_earth(self.PHInorm_cg, Tgeo2body) - - # EoM - d2Ucg_dt2, Nxyz = self.rigid_EoM(dUcg_dt, Pb, g_cg, modus) - Pf = np.dot(self.PHIkf.T, Pk_aero) + self.Mfcg.dot(np.hstack((d2Ucg_dt2[0:3] - g_cg, d2Ucg_dt2[3:6]))) + Pf_ext - d2Uf_dt2 = self.flexible_EoM(dUf_dt, Uf, Pf) - - # CS derivatives - dcommand = self.get_command_derivatives(t, X, Vtas, gamma, alpha, beta, Nxyz, np.dot(Tbody2geo, X[6:12])[0:3]) - - # output - Y = np.hstack((np.dot(Tbody2geo, X[6:12]), - np.dot(self.PHIcg_norm, d2Ucg_dt2), - dUf_dt, - d2Uf_dt2, - dcommand, - Nxyz[2], - Vtas, - beta, - )) - - if modus in ['trim', 'sim']: - return Y - elif modus in ['trim_full_output']: - response = {'X': X, - 'Y': Y, - 't': np.array([t]), - 'Pk_rbm': Pk_rbm, - 'Pk_cam': Pk_cam, - 'Pk_aero': Pk_aero, - 'Pk_cs': Pk_cs, - 'Pk_f': Pk_f, - 'Pk_gust': Pk_gust, - 'Pk_unsteady': Pk_unsteady, - 'Pk_idrag': Pk_idrag, - 'q_dyn': np.array([q_dyn]), - 'Pb': Pb, - 'Pmac': Pmac, - 'Pf': Pf, - 'alpha': np.array([alpha]), - 'beta': np.array([beta]), - 'Ux2': Ux2, - 'dUcg_dt': dUcg_dt, - 'd2Ucg_dt2': d2Ucg_dt2, - 'Uf': Uf, - 'dUf_dt': dUf_dt, - 'd2Uf_dt2': d2Uf_dt2, - 'Nxyz': Nxyz, - 'g_cg': g_cg, - 'Pextra': Pextra, - } - return response - elif modus in ['sim_full_output']: - # For time domain simulations, typically not all results are required. To reduce the amount of data while - # maintaining compatibility with the trim, empty arrays are used. - response = {'X': X, - 'Y': Y, - 't': np.array([t]), - 'Pk_aero': Pk_aero, - 'q_dyn': np.array([q_dyn]), - 'Pb': Pb, - 'Pmac': Pmac, - 'alpha': np.array([alpha]), - 'beta': np.array([beta]), - 'Ux2': Ux2, - 'dUcg_dt': dUcg_dt, - 'd2Ucg_dt2': d2Ucg_dt2, - 'Uf': Uf, - 'dUf_dt': dUf_dt, - 'd2Uf_dt2': d2Uf_dt2, - 'Nxyz': Nxyz, - 'g_cg': g_cg, - 'Pextra': Pextra, - } - return response - - def eval_equations(self, X_free, time, modus='trim_full_output'): - # this is a wrapper for the model equations 'eqn_basic' - if modus in ['trim', 'trim_full_output']: - # get inputs from trimcond and apply inputs from fsolve - X = np.array(self.trimcond_X[:, 2], dtype='float') - X[np.where((self.trimcond_X[:, 1] == 'free'))[0]] = X_free - elif modus in ['sim', 'sim_full_output']: - X = X_free - - # evaluate model equations - if modus == 'trim': - Y = self.equations(X, time, 'trim') - # get the current values from Y and substract tamlab.figure() - # fsolve only finds the roots; Y = 0 - Y_target_ist = Y[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] - Y_target_soll = np.array(self.trimcond_Y[:, 2], dtype='float')[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] - out = Y_target_ist - Y_target_soll - return out - - elif modus == 'sim': - Y = self.equations(X, time, 'sim') - return Y[self.solution.idx_state_derivatives + self.solution.idx_input_derivatives] - - elif modus == 'sim_full_output': - response = self.equations(X, time, 'sim_full_output') - return response - - elif modus == 'trim_full_output': - response = self.equations(X, time, 'trim_full_output') - return response - - def eval_equations_iteratively(self, X_free, time, modus='trim_full_output'): - # this is a wrapper for the model equations - - # get inputs from trimcond and apply inputs from fsolve - X = np.array(self.trimcond_X[:, 2], dtype='float') - X[np.where((self.trimcond_X[:, 1] == 'free'))[0]] = X_free - logging.debug('X_free: {}'.format(X_free)) - converged = False - inner_loops = 0 - while not converged: - inner_loops += 1 - response = self.equations(X, time, 'trim_full_output') - logging.info('Inner iteration {}, calculate structural deformations.'.format(self.counter)) - Uf_new = linalg.solve(self.Kff, response['Pf']) - - # Add a relaxation factor between each loop to reduce overshoot and/or oscillations. - # In case the solution hasn't converged in a reasonable time (say 8 loops), increase the relaxation. - # A low relaxation factor is slower but more robust, f_relax = 1.0 implies no relaxation. - # After 16 inner loops, decrease the convergence criterion epsilon. - if inner_loops < 8: - f_relax = 0.8 - else: - f_relax = 0.5 - if inner_loops < 16: - epsilon = self.jcl.general['b_ref'] * 1.0e-5 - else: - epsilon = self.jcl.general['b_ref'] * 1.0e-4 - logging.info(' - Relaxation factor: {}'.format(f_relax)) - logging.info(' - Epsilon: {:0.6g}'.format(epsilon)) - - # recover Uf_old from last step and blend with Uf_now - Uf_old = [self.trimcond_X[np.where((self.trimcond_X[:, 0] == 'Uf' + str(i_mode)))[0][0], 2] for i_mode in - range(1, self.n_modes + 1)] - Uf_old = np.array(Uf_old, dtype='float') - Uf_new = Uf_new * f_relax + Uf_old * (1.0 - f_relax) - - # set new values for Uf in trimcond for next loop and store in response - for i_mode in range(self.n_modes): - self.trimcond_X[np.where((self.trimcond_X[:, 0] == 'Uf' + str(i_mode + 1)))[0][0], 2] = '{:g}'.format( - Uf_new[i_mode]) - response['X'][12 + i_mode] = Uf_new[i_mode] - - # convergence parameter for iterative evaluation - Ug_f_body = np.dot(self.PHIf_strc.T, Uf_new.T).T - defo_new = Ug_f_body[self.strcgrid['set'][:, :3]].max() # Groesste Verformung, meistens Fluegelspitze - ddefo = defo_new - self.defo_old - self.defo_old = np.copy(defo_new) - if np.abs(ddefo) < epsilon: - converged = True - logging.info(' - Max. deformation: {:0.6g}, delta: {:0.6g} smaller than epsilon, converged.'.format( - defo_new, ddefo)) - else: - logging.info(' - Max. deformation: {:0.6g}, delta: {:0.6g}'.format(defo_new, ddefo)) - if inner_loops > 20: - raise ConvergenceError('No convergence of structural deformation achieved after {} inner loops. \ - Check convergence of CFD solution and/or convergence criterion "delta".'.format(inner_loops)) - # get the current values from Y and substract tamlab.figure() - # fsolve only finds the roots; Y = 0 - Y_target_ist = response['Y'][np.where((self.trimcond_Y[:, 1] == 'target'))[0]] - Y_target_soll = np.array(self.trimcond_Y[:, 2], dtype='float')[np.where((self.trimcond_Y[:, 1] == 'target'))[0]] - out = Y_target_ist - Y_target_soll - - if modus in ['trim']: - return out - elif modus == 'trim_full_output': - return response diff --git a/loadskernel/equations/unsteady.py b/loadskernel/equations/unsteady.py deleted file mode 100644 index 255e030d..00000000 --- a/loadskernel/equations/unsteady.py +++ /dev/null @@ -1,130 +0,0 @@ -import numpy as np -from loadskernel.equations.common import Common -from loadskernel.solution_tools import gravitation_on_earth - - -class Unsteady(Common): - - def equations(self, X, t, modus): - self.counter += 1 - # recover states - Tgeo2body, Tbody2geo = self.geo2body(X) - dUcg_dt, Uf, dUf_dt = self.recover_states(X) - Vtas, q_dyn = self.recover_Vtas(self.X0) - onflow = self.recover_onflow(X) - alpha, beta, gamma = self.windsensor(X, Vtas, Uf, dUf_dt) - Ux2 = self.get_Ux2(X) - - # aerodynamics - Pk_rbm, wj_rbm = self.rbm(onflow, q_dyn, Vtas) - Pk_cam, wj_cam = self.camber_twist(q_dyn) - Pk_cs, wj_cs = self.cs(X, Ux2, q_dyn) - Pk_f, wj_f = self.flexible(Uf, dUf_dt, onflow, q_dyn, Vtas) - Pk_gust, wj_gust = self.gust(X, q_dyn) - - wj = wj_rbm + wj_cam + wj_cs + wj_f + wj_gust - Pk_idrag = self.idrag(wj, q_dyn) - Pk_unsteady, dlag_states_dt = self.unsteady(X, t, wj, Uf, dUf_dt, onflow, q_dyn, Vtas) - - Pextra, Pb_ext, Pf_ext = self.engine(X, Vtas, q_dyn, Uf, dUf_dt, t) - - # correction coefficients - Pb_corr = self.correctioon_coefficients(alpha, beta, q_dyn) - - # summation of forces - Pk_aero = Pk_rbm + Pk_cam + Pk_cs + Pk_f + Pk_gust + Pk_idrag + Pk_unsteady - Pmac = np.dot(self.Dkx1.T, Pk_aero) - Pb = np.dot(self.PHImac_cg.T, Pmac) + Pb_corr + Pb_ext - - g_cg = gravitation_on_earth(self.PHInorm_cg, Tgeo2body) - - # EoM - d2Ucg_dt2, Nxyz = self.rigid_EoM(dUcg_dt, Pb, g_cg, modus) - Pf = np.dot(self.PHIkf.T, Pk_aero) + self.Mfcg.dot(np.hstack((d2Ucg_dt2[0:3] - g_cg, d2Ucg_dt2[3:6]))) + Pf_ext - d2Uf_dt2 = self.flexible_EoM(dUf_dt, Uf, Pf) - - # CS derivatives - dcommand = self.get_command_derivatives(t, X, Vtas, gamma, alpha, beta, Nxyz, np.dot(Tbody2geo, X[6:12])[0:3]) - - # output - Y = np.hstack((np.dot(Tbody2geo, X[6:12]), - np.dot(self.PHIcg_norm, d2Ucg_dt2), - dUf_dt, - d2Uf_dt2, - dcommand, - dlag_states_dt, - Nxyz[2], - Vtas, - beta, - )) - - if modus in ['trim', 'sim']: - return Y - elif modus in ['trim_full_output']: - response = {'X': X, - 'Y': Y, - 't': np.array([t]), - 'Pk_rbm': Pk_rbm, - 'Pk_cam': Pk_cam, - 'Pk_aero': Pk_aero, - 'Pk_cs': Pk_cs, - 'Pk_f': Pk_f, - 'Pk_gust': Pk_gust, - 'Pk_unsteady': Pk_unsteady, - 'Pk_idrag': Pk_idrag, - 'q_dyn': np.array([q_dyn]), - 'Pb': Pb, - 'Pmac': Pmac, - 'Pf': Pf, - 'alpha': np.array([alpha]), - 'beta': np.array([beta]), - 'Ux2': Ux2, - 'dUcg_dt': dUcg_dt, - 'd2Ucg_dt2': d2Ucg_dt2, - 'Uf': Uf, - 'dUf_dt': dUf_dt, - 'd2Uf_dt2': d2Uf_dt2, - 'Nxyz': Nxyz, - 'g_cg': g_cg, - 'Pextra': Pextra, - } - return response - elif modus in ['sim_full_output']: - # For time domain simulations, typically not all results are required. To reduce the amount of data while - # maintaining compatibility with the trim, empty arrays are used. - response = {'X': X, - 'Y': Y, - 't': np.array([t]), - 'Pk_aero': Pk_aero, - 'Pk_gust': Pk_gust, - 'Pk_unsteady': Pk_unsteady, - 'q_dyn': np.array([q_dyn]), - 'Pb': Pb, - 'Pmac': Pmac, - 'alpha': np.array([alpha]), - 'beta': np.array([beta]), - 'Ux2': Ux2, - 'dUcg_dt': dUcg_dt, - 'd2Ucg_dt2': d2Ucg_dt2, - 'Uf': Uf, - 'dUf_dt': dUf_dt, - 'd2Uf_dt2': d2Uf_dt2, - 'Nxyz': Nxyz, - 'g_cg': g_cg, - 'Pextra': Pextra, - } - return response - - def eval_equations(self, X_free, time, modus='trim_full_output'): - if modus in ['sim', 'sim_full_output']: - X = X_free - - # evaluate model equations - if modus == 'sim': - Y = self.equations(X, time, 'sim') - return Y[self.solution.idx_state_derivatives + self.solution.idx_input_derivatives - + self.solution.idx_lag_derivatives] - - elif modus == 'sim_full_output': - response = self.equations(X, time, 'sim_full_output') - return response diff --git a/loadskernel/gather_loads.py b/loadskernel/gather_loads.py index 6ff87b14..a7c25fe1 100644 --- a/loadskernel/gather_loads.py +++ b/loadskernel/gather_loads.py @@ -35,81 +35,90 @@ def __init__(self, jcl, model): } def gather_monstations(self, trimcase, response): - logging.info('gathering information on monitoring stations from response(s)...') - for i_station in range(self.mongrid['n']): - name = self.mongrid['name'][i_station] - subcase = str(trimcase['subcase']) - # Unterscheidung zwischen Trim und Zeit-Simulation, da die Dimensionen der response anders sind (n_step x n_value) - if len(response['t']) > 1: - loads = response['Pmon_local'][:, self.mongrid['set'][i_station, :]] - # save time data per subcase - self.monstations[name][str(response['i'])] = {'subcase': subcase, - 'loads': loads, - 't': response['t']} - else: - loads = response['Pmon_local'][0, self.mongrid['set'][i_station, :]] - self.monstations[name]['subcases'].append(subcase) - self.monstations[name]['loads'].append(loads) - self.monstations[name]['t'].append(response['t'][0]) + if 'Pmon_local' in response: + logging.info('gathering information on monitoring stations from response(s)...') + for i_station in range(self.mongrid['n']): + name = self.mongrid['name'][i_station] + subcase = str(trimcase['subcase']) + # Unterscheidung zwischen Trim und Zeit-Simulation, da die Dimensionen der response anders sind + # (n_step x n_value) + if len(response['t']) > 1: + loads = response['Pmon_local'][:, self.mongrid['set'][i_station, :]] + # save time data per subcase + self.monstations[name][str(response['i'])] = {'subcase': subcase, + 'loads': loads, + 't': response['t']} + else: + loads = response['Pmon_local'][0, self.mongrid['set'][i_station, :]] + self.monstations[name]['subcases'].append(subcase) + self.monstations[name]['loads'].append(loads) + self.monstations[name]['t'].append(response['t'][0]) - if 'Pmon_turb' in response: - # Check if there are any limit turbulence loads available in the response. - # If yes, map them into the monstations. - loads = response['Pmon_turb'][0, self.mongrid['set'][i_station, :]] - correlations = response['correlations'][self.mongrid['set'][i_station, - :], :][:, self.mongrid['set'][i_station, :]] + if 'Pmon_turb' in response: + # Check if there are any limit turbulence loads available in the response. + # If yes, map them into the monstations. + loads = response['Pmon_turb'][0, self.mongrid['set'][i_station, :]] + correlations = response['correlations'][self.mongrid['set'] + [i_station, :], :][:, self.mongrid['set'][i_station, :]] - self.monstations[name]['turbulence_loads'].append(loads) - self.monstations[name]['correlations'].append(correlations) + self.monstations[name]['turbulence_loads'].append(loads) + self.monstations[name]['correlations'].append(correlations) def gather_dyn2stat(self, response): """ - Schnittlasten an den Monitoring Stationen raus schreiben (z.B. zum Plotten) - Knotenlasten raus schreiben (weiterverarbeitung z.B. als FORCE und MOMENT Karten fuer Nastran) + Write out section loads as well and nodal loads, for example for plotting and for FORCE and MOMENT + cards in Nastran. The procedure involves searching for min and max loads over time at each monitoring + station and gathering the corresponding time slices / snapshots. + The dyn2stat loads are also useful if a dynamic simulations are performed and the response file becomes + too large to handle. """ - logging.info('searching min/max in time data at {} monitoring stations ' - 'and gathering loads (dyn2stat)...'.format(len(self.monstations.keys()))) - if len(response['t']) > 1: - i_case = str(response['i']) - timeslices_dyn2stat = np.array([], dtype=int) - for key, monstation in self.monstations.items(): - pos_max_loads_over_time = np.argmax(monstation[i_case]['loads'], 0) - pos_min_loads_over_time = np.argmin(monstation[i_case]['loads'], 0) - """ - Although the time-based approach considers all DoFs, it might lead to fewer time slices / snapshots - compared to Fz,min/max, Mx,min/max, ..., because identical time slices are avoided. - """ - # Remember identified time slices - timeslices_dyn2stat = np.concatenate((timeslices_dyn2stat, pos_max_loads_over_time, pos_min_loads_over_time)) - logging.info('reducing dyn2stat data...') - timeslices_dyn2stat = np.unique(timeslices_dyn2stat) - nastran_subcase_running_number = 1 - for pos in timeslices_dyn2stat: - # save nodal loads Pg for this time slice - self.dyn2stat['Pg'].append(response['Pg'][pos, :]) - subcases_dyn2stat_string = str(self.monstations[key][i_case]['subcase']) + '_t={:.3f}'.format( - self.monstations[key][i_case]['t'][pos, 0]) - self.dyn2stat['subcases'].append(subcases_dyn2stat_string) - """ - Generate unique IDs for subcases: - Take first digits from original subcase, then add a running number. This is a really stupid approach, - but we are limited to 8 digits and need to make the best out of that... Using 5 digits for the subcase - and 3 digits for the running number appears to work for most cases. This is only important for Nastran - Force and Moment cards as Nastran does not like subcases in the style 'xy_t=0.123' and requires numbers - in ascending order. - """ - self.dyn2stat['subcases_ID'].append(int(self.monstations[key][i_case]['subcase']) - * 1000 + nastran_subcase_running_number) - # save section loads to monstations + if 'Pg' in response: + logging.info('searching min/max in time data at %s monitoring stations ' + 'and gathering loads (dyn2stat)...', len(self.monstations.keys())) + if len(response['t']) > 1: + i_case = str(response['i']) + timeslices_dyn2stat = np.array([], dtype=int) for key, monstation in self.monstations.items(): - monstation['subcases'].append(subcases_dyn2stat_string) - monstation['loads'].append(monstation[i_case]['loads'][pos, :]) - monstation['t'].append(monstation[i_case]['t'][pos, :]) - # increase counter of running number by 1 - nastran_subcase_running_number += 1 + pos_max_loads_over_time = np.argmax(monstation[i_case]['loads'], 0) + pos_min_loads_over_time = np.argmin(monstation[i_case]['loads'], 0) + """ + Although the time-based approach considers all DoFs, it might lead to fewer time slices / snapshots + compared to Fz,min/max, Mx,min/max, ..., because identical time slices are avoided. + """ + # Remember identified time slices + timeslices_dyn2stat = np.concatenate((timeslices_dyn2stat, pos_max_loads_over_time, + pos_min_loads_over_time)) + logging.info('reducing dyn2stat data...') + timeslices_dyn2stat = np.unique(timeslices_dyn2stat) + nastran_subcase_running_number = 1 + for pos in timeslices_dyn2stat: + # Save nodal loads Pg for this time slice + self.dyn2stat['Pg'].append(response['Pg'][pos, :]) + # For the next operation, just take the first monstation to get the subcase name and time + key = list(self.monstations.keys())[0] + subcases_dyn2stat_string = f'{self.monstations[key][i_case]['subcase']}_t={self.monstations[key] + [i_case]['t'][pos, 0]:.3f}' + self.dyn2stat['subcases'].append(subcases_dyn2stat_string) + """ + Generate unique IDs for subcases: + Take first digits from original subcase, then add a running number. This is a really stupid approach, + but we are limited to 8 digits and need to make the best out of that... Using 5 digits for the subcase + and 3 digits for the running number appears to work for most cases. This is only important for Nastran + Force and Moment cards as Nastran does not like subcases in the style 'xy_t=0.123' and requires numbers + in ascending order. + """ + self.dyn2stat['subcases_ID'].append(int(self.monstations[key][i_case]['subcase']) + * 1000 + nastran_subcase_running_number) + # Save section loads to monstations + for key, monstation in self.monstations.items(): + monstation['subcases'].append(subcases_dyn2stat_string) + monstation['loads'].append(monstation[i_case]['loads'][pos, :]) + monstation['t'].append(monstation[i_case]['t'][pos, :]) + # Increase counter of running number by 1 + nastran_subcase_running_number += 1 - else: - i_case = response['i'] - self.dyn2stat['Pg'].append(response['Pg'][0, :]) - self.dyn2stat['subcases'].append(str(self.jcl.trimcase[i_case]['subcase'])) - self.dyn2stat['subcases_ID'].append(int(self.jcl.trimcase[i_case]['subcase'])) + else: + i_case = response['i'] + self.dyn2stat['Pg'].append(response['Pg'][0, :]) + self.dyn2stat['subcases'].append(str(self.jcl.trimcase[i_case]['subcase'])) + self.dyn2stat['subcases_ID'].append(int(self.jcl.trimcase[i_case]['subcase'])) diff --git a/loadskernel/interpolate.py b/loadskernel/interpolate.py index 5ac8c38e..1cfecd15 100644 --- a/loadskernel/interpolate.py +++ b/loadskernel/interpolate.py @@ -31,9 +31,9 @@ def __call__(self, i): def calc_gradients(self): # calculate the gradients between all samples for m in range(self.n_samples - 1): - self.gradients[m, :, :] = (self.data[m + 1, :, :] - self.data[m, :, :]) / (self.x[m + 1] - self.x[m]) + self.gradients[[m], :] = (self.data[[m + 1], :] - self.data[[m], :]) / (self.x[m + 1] - self.x[m]) # repeat the last set of gradients to make sure there are right-sided gradients available in case of extrapolation - self.gradients[-1, :, :] = self.gradients[-2, :, :] + self.gradients[[-1], :] = self.gradients[[-2], :] def interpolate(self, i): # find the nearest neighbor @@ -46,4 +46,5 @@ def interpolate(self, i): else: pos_gradients = pos # perform linear interpolation and return result - return self.data[pos, :, :] + self.gradients[pos_gradients, :, :] * delta + result = self.data[[pos], :] + self.gradients[[pos_gradients], :] * delta + return result.squeeze() diff --git a/loadskernel/io_functions/data_handling.py b/loadskernel/io_functions/data_handling.py index 7be1a665..b0ceeca7 100644 --- a/loadskernel/io_functions/data_handling.py +++ b/loadskernel/io_functions/data_handling.py @@ -40,6 +40,10 @@ def open_hdf5(filename): return h5py.File(filename, 'w') +def append_hdf5(filename): + return h5py.File(filename, 'a') + + def write_hdf5(fid, dic, path=''): recursively_save_dict_to_hdf5(fid, dic, path) @@ -85,7 +89,7 @@ def recursively_save_dict_to_hdf5(fid, dic, path=''): def load_hdf5_responses(job_name, path_output): - logging.info('--> Opening response(s).') + logging.info('Opening response(s).') filename = path_output + 'response_' + job_name + '.hdf5' fid = load_hdf5(filename) response = [fid[key] for key in sorted(fid.keys(), key=int) if fid[key]['successful'][()]] @@ -130,7 +134,7 @@ def load_hdf5_sparse_matrix(hdf5_group): def load_jcl(job_name, path_input, jcl): if jcl is None: - logging.info('--> Reading parameters from JCL.') + logging.info('Reading parameters from JCL.') # import jcl dynamically by filename # this is the newer way used in Python 3 spec = importlib.util.spec_from_file_location('jcl', os.path.join(path_input, job_name + '.py')) @@ -146,7 +150,7 @@ def load_jcl(job_name, path_input, jcl): def gather_responses(job_name, path): - logging.info('--> Gathering response(s).') + logging.info('Gathering response(s).') filenames = os.listdir(path) filenames.sort() response = [] diff --git a/loadskernel/model.py b/loadskernel/model.py index ffe29373..972e6abb 100644 --- a/loadskernel/model.py +++ b/loadskernel/model.py @@ -14,6 +14,7 @@ from loadskernel import spline_rules from loadskernel import spline_functions from loadskernel import build_splinegrid +from loadskernel.io_functions import data_handling from loadskernel.io_functions import read_mona, read_op4, read_bdf from loadskernel.io_functions import read_cfdgrids from loadskernel import grid_trafo @@ -30,6 +31,8 @@ def __init__(self, jcl, path_output): self.path_output = path_output # init the bdf reader self.bdf_reader = read_bdf.Reader() + # initialize empty dict for GAFs + self.GAFs = {} def build_model(self): # run all build function/stages @@ -200,8 +203,9 @@ def build_atmo(self): def build_aero(self): logging.info('Building aero model...') - if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'hybrid', 'nonlin_steady', - 'cfd_steady', 'cfd_unsteady', 'freq_dom']: + if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'nonlin_steady', + 'freq_dom', 'mona_freq_dom', + 'cfd_steady', 'cfd_unsteady', 'cfd_freq_dom']: self.build_aerogrid() self.build_aero_matrices() self.build_W2GJ() @@ -209,12 +213,38 @@ def build_aero(self): self.build_cs() self.build_AICs_steady() self.build_AICs_unsteady() + self.read_precomputed_GAFs() else: logging.error('Unknown aero method: ' + str(self.jcl.aero['method'])) logging.info('The aerodynamic model consists of {} panels and {} control surfaces.'.format( self.aerogrid['n'], len(self.x2grid['ID_surf']))) + def read_precomputed_GAFs(self): + # The CFD-based frequency domain method requires precomputed GAFs from the responses of another job. + # The corresponding job name is given in the jcl. + if self.jcl.aero['method'] in ['cfd_freq_dom'] and 'job_name_gafs' in self.jcl.aero: + responses = data_handling.load_hdf5_responses(self.jcl.aero['job_name_gafs'], self.path_output) + logging.info('Moving GAFs from responses into model:') + for response in responses: + if response['successful'] and 'pulse_signal' in response: + # Write info about which GAFs we found in the response + key = '.'.join(response['desc'].asstr()[()].split('.')[:-1]) + logging.info(' - %s', key) + # Pick GAF matrices from response and copy into model + self.GAFs[key] = {} + gaf_items = ['k_red', 'Qhk', 'Qhh', 'Qk_gust', 'q_dyn'] + for item in gaf_items: + self.GAFs[key][item] = response[item][()] + # Copy the linearization point + self.GAFs[key]['X0'] = response['X'][()].squeeze() + # Pick relevant data from the linearization point and copy into model + self.GAFs[key]['response'] = {} + resp_items = ['X', 't', 'Pk_aero', 'Pk_gust', 'Pk_unsteady', + 'dUcg_dt', 'd2Ucg_dt2', 'Uf', 'dUf_dt', 'd2Uf_dt2', 'g_cg'] + for item in resp_items: + self.GAFs[key]['response'][item] = response[item][()].squeeze() + def build_aerogrid(self): # To avoid interference with other CQUAD4 cards parsed earlier, clear those dataframes first if 'method_caero' in self.jcl.aero and self.jcl.aero['method_caero'] == 'CQUAD4': @@ -351,7 +381,7 @@ def build_AICs_unsteady(self): self.build_rfa() else: logging.error('Unknown AIC method: ' + str(self.jcl.aero['method_AIC'])) - elif self.jcl.aero['method'] in ['freq_dom']: + elif self.jcl.aero['method'] in ['freq_dom', 'mona_freq_dom']: if self.jcl.aero['method_AIC'] == 'dlm': self.build_AICs_DLM() elif self.jcl.aero['method_AIC'] == 'nastran': diff --git a/loadskernel/plotting_standard.py b/loadskernel/plotting_standard.py index 187dcb0c..1b991b16 100755 --- a/loadskernel/plotting_standard.py +++ b/loadskernel/plotting_standard.py @@ -4,6 +4,7 @@ import os import numpy as np from scipy.spatial import ConvexHull +from scipy.fftpack import fft, fftfreq from matplotlib import pyplot as plt from matplotlib.backends.backend_pdf import PdfPages @@ -80,7 +81,7 @@ def plot_monstations(self, filename_pdf): if self.cuttingforces_fuselage: self.cuttingforces_along_axis_plots(monstations=self.cuttingforces_fuselage, axis=0) self.pp.close() - logging.info('Plots saved as ' + filename_pdf) + logging.info('Plots saved as %s', filename_pdf) def potato_plot(self, station, desc, color, dof_xaxis, dof_yaxis, show_hull=True, show_labels=False, show_minmax=False): loads = np.array(self.monstations[station]['loads']) @@ -102,7 +103,7 @@ def potato_plot(self, station, desc, color, dof_xaxis, dof_yaxis, show_hull=True self.subplot.plot(points[simplex, 0], points[simplex, 1], color=color, linewidth=2.0, linestyle='--') crit_trimcases = [subcases[i] for i in hull.vertices] if show_labels: - for i_case in range(crit_trimcases.__len__()): + for i_case in range(len(crit_trimcases)): self.subplot.text(points[hull.vertices[i_case], 0], points[hull.vertices[i_case], 1], str(subcases[hull.vertices[i_case]]), fontsize=8) except Exception: @@ -197,7 +198,7 @@ def potato_plots(self): def cuttingforces_along_axis_plots(self, monstations, axis): assert axis in [0, 1, 2], 'Plotting along an axis only supported for axis 0, 1 or 2!' - logging.info('Start plotting cutting forces along axis {}...'.format(axis)) + logging.info('Start plotting cutting forces along axis %d...', axis) # Read the data required for plotting. loads = [] offsets = [] @@ -345,7 +346,7 @@ def plot_monstations_time(self, filename_pdf): pp.savefig() plt.close() pp.close() - logging.info('plots saved as ' + filename_pdf) + logging.info('plots saved as %s', filename_pdf) class FlutterPlots(LoadPlots): @@ -442,7 +443,7 @@ def plot_eigenvalues(self): # this kind of plot is only feasible for methods which iterate over Vtas, e.g. not the K- or KE-methods if 'flutter' in simcase and simcase['flutter_para']['method'] not in ['pk_rodden', 'pk_schwochow', 'pk', 'statespace']: - logging.info('skip plotting of eigenvalues and -vectors for {}'.format(trimcase['desc'])) + logging.info('skip plotting of eigenvalues and -vectors for %s', trimcase['desc']) continue # Plot boundaries @@ -455,13 +456,14 @@ def plot_eigenvalues(self): colors = itertools.cycle((plt.cm.tab20c(np.linspace(0, 1, 20)))) markers = itertools.cycle(('+', 'o', 'v', '^', '<', '>', '8', 's', 'p', '*', 'x', 'D',)) desc = [str(mode) for mode in range(response['eigenvalues'].shape[1])] - + # clear all axes for next plot ax[0].cla() ax[1].cla() ax[2].cla() - ax_cbar.cla(), ax_freq.cla() # clear all axes for next plot + ax_cbar.cla() + ax_freq.cla() # plot eigenvector - im_eig = ax[2].imshow(response['eigenvectors'][i].__abs__(), cmap='hot_r', aspect='auto', + im_eig = ax[2].imshow(np.abs(response['eigenvectors'][i]), cmap='hot_r', aspect='auto', origin='upper', vmin=0.0, vmax=1.0) # add colorbar to plot fig.colorbar(im_eig, cax=ax_cbar, orientation="horizontal") @@ -477,12 +479,12 @@ def plot_eigenvalues(self): color=color, linestyle='--') ax[1].plot(response['eigenvalues'][i, j].real, response['eigenvalues'][i, j].imag, marker=marker, markersize=8.0, color=color, label=desc[j]) - ax[2].plot(j, response['states'].__len__(), + ax[2].plot(j, len(response['states']), marker=marker, markersize=8.0, c=color) # make plots nice - fig.suptitle(t='{}, Veas={:.2f} m/s, Vtas={:.2f} m/s'.format( - trimcase['desc'], tas2eas(response['Vtas'][i, 0], h), response['Vtas'][i, 0]), fontsize=16) + fig.suptitle(t=f'{trimcase['desc']}, Veas={tas2eas(response['Vtas'][i, 0], h):.2f} m/s, \ + Vtas={response['Vtas'][i, 0]:.2f} m/s', fontsize=16) ax[0].set_position([0.12, 0.1, 0.25, 0.8]) ax[0].set_xlabel('Real') ax[0].set_ylabel('Imag') @@ -514,7 +516,7 @@ def plot_eigenvalues(self): ax[1].legend(bbox_to_anchor=(1.10, 1), loc='upper left', borderaxespad=0.0, fontsize=10) ax[2].set_position([0.60, 0.1, 0.35, 0.8]) - ax[2].yaxis.set_ticks(np.arange(0, response['states'].__len__(), 1)) + ax[2].yaxis.set_ticks(np.arange(0, len(response['states']), 1)) ax[2].yaxis.set_ticklabels(response['states'].asstr(), fontsize=10) ax[2].yaxis.set_tick_params(rotation=0) ax[2].xaxis.set_ticks(np.arange(0, response['eigenvalues'].shape[1], 1)) @@ -530,13 +532,13 @@ def plot_fluttercurves_to_pdf(self, filename_pdf): self.pp = PdfPages(filename_pdf) self.plot_fluttercurves() self.pp.close() - logging.info('plots saved as ' + filename_pdf) + logging.info('plots saved as %s', filename_pdf) def plot_eigenvalues_to_pdf(self, filename_pdf): self.pp = PdfPages(filename_pdf) self.plot_eigenvalues() self.pp.close() - logging.info('plots saved as ' + filename_pdf) + logging.info('plots saved as %s', filename_pdf) class TurbulencePlots(LoadPlots): @@ -547,7 +549,7 @@ def plot_monstations(self, filename_pdf): self.pp = PdfPages(filename_pdf) self.potato_plots() self.pp.close() - logging.info('plots saved as ' + filename_pdf) + logging.info('plots saved as %s', filename_pdf) def potato_plot(self, station, desc, color, dof_xaxis, dof_yaxis, show_hull=True, show_labels=False, show_minmax=False): loads = np.array(self.monstations[station]['loads']) @@ -608,7 +610,7 @@ def fit_ellipse(self, X0, Y0, X, Y, color): # Print the equation of the ellipse in standard form logging.debug( - 'The ellipse is given by {0:.3}x^2 + {1:.3}*2xy+{2:.3}y^2 = 1'.format(x[0], x[1], x[2])) + 'The ellipse is given by %.3fx^2 + %.3f*2xy+%.3fy^2 = 1', x[0], x[1], x[2]) # Calculate the parameters of the ellipse alpha = -0.5 * np.arctan(2 * x[1] / (x[2] - x[0])) @@ -617,9 +619,9 @@ def fit_ellipse(self, X0, Y0, X, Y, color): major = (2.0 / (eta - zeta)) ** 0.5 minor = (2.0 / (eta + zeta)) ** 0.5 - logging.debug('Major axis = {:.3f}'.format(major)) - logging.debug('Minor axis = {:.3f}'.format(minor)) - logging.debug('Rotation = {:.3f} deg'.format(alpha / np.pi * 180.0)) + logging.debug('Major axis = %.3f', major) + logging.debug('Minor axis = %.3f', minor) + logging.debug('Rotation = %.3f deg', alpha / np.pi * 180.0) # Plot the given samples # self.subplot.scatter(X, Y, label='Data Points') @@ -631,3 +633,140 @@ def ellipse_polar(self, major, minor, alpha, phi=np.linspace(0, 2.0 * np.pi, 360 X = 0.0 + major * np.cos(phi) * np.cos(alpha) - minor * np.sin(phi) * np.sin(alpha) Y = 0.0 + major * np.cos(phi) * np.sin(alpha) + minor * np.sin(phi) * np.cos(alpha) return X, Y + + +class GAFPlots(LoadPlots): + + def make_as_nice(self, ax): + for a in ax: + a.get_yaxis().set_label_coords(x=-0.10, y=0.5) + a.grid(visible=True, which='major', axis='both') + a.minorticks_on() + a.legend(loc='upper right') + + def plot_pluse(self, filename_pdf): + logging.info('start plotting pulse signals...') + pp = PdfPages(filename_pdf) + for response in self.responses: + # Get/load data from the response + trimcase = self.jcl.trimcase[response['i'][()]] + simcase = self.jcl.simcase[response['i'][()]] + pulse_signal = response['pulse_signal'][()] + gust_signal = response['gust_signal'][()] + t = response['t_pulse'][()] + Pb_pulse = response['Pb_pulse'][()] + Pb_gust = response['Pb_gust'][()] + # Qhh = response['Qhh'][()] + k_red = response['k_red'][()] + Vtas = sum(response['X'][0, 6:9] ** 2) ** 0.5 + q_dyn = response['q_dyn'][0] + c_ref = self.jcl.general['c_ref'] + + # Step 1: plot pulse in frequency domain + pulse_f = fft(pulse_signal) + gust_f = fft(gust_signal) + n_freqs = int(simcase['gaf_para']['fmax'] / simcase['gaf_para']['df']) + if n_freqs % 2 != 0: # n_freq is odd + n_freqs += 1 # make even + # Calculate all parameters from the number of freqs + fmax = n_freqs * simcase['gaf_para']['df'] + dt = 1.0 / fmax + # Whole frequency space including negative frequencies + fftfreqs = fftfreq(n_freqs, dt) + # Positive only frequencies where we need to calculate the TFs and excitations + positiv_fftfreqs = np.abs(fftfreqs[:n_freqs // 2 + 1]) + # Do the actual plotting + fig, ax = plt.subplots(2, sharex=True, figsize=(8, 10)) + fig.suptitle(f'{trimcase['desc']}', fontsize=16) + ax[0].set_position([0.15, 0.55, 0.75, 0.35]) + ax[0].plot(positiv_fftfreqs, np.abs(pulse_f[0, :len(positiv_fftfreqs)]), '-', label='FFT Pulse') + ax[1].set_position([0.15, 0.15, 0.75, 0.35]) + ax[1].plot(positiv_fftfreqs, np.abs(gust_f[:len(positiv_fftfreqs)]), '-', label='FFT Gust') + # the shared x-axis + ax[1].set_xlabel('Frequency [Hz]') + f_max = 50.0 # I don't like hard-coded limits here, but 50 Hz is ok for most applications + k_max = 2.0 * np.pi * f_max * c_ref / 2.0 / Vtas + ax[1].set_xlim((0.0, f_max)) + self.make_as_nice(ax) + # set second x-axis + ax_k = ax[-1].twiny() + ax_k.set_position([0.15, 0.15, 0.75, 0.35]) # Move the axis on top of the lower plot + ax_k.xaxis.set_ticks_position('bottom') + ax_k.xaxis.set_label_position('bottom') + ax_k.spines['bottom'].set_position(('outward', 60)) + ax_k.set_xlim((0.0, k_max)) + ax_k.minorticks_on() + ax_k.set_xlabel('$k_{red}$') + pp.savefig() + plt.close() + + # Step 2: plot lift coefficient derivative from heave and pitch motions as suggesrted by Marc-Johan + # Fourier transformation + Pb_f = fft(Pb_pulse) + Qhb = Pb_f / pulse_f + # Calculate dCl from the lift force Fz at CG + omega = k_red * Vtas * 2.0 / c_ref + dCl_heave = Qhb[2, 1, :len(k_red)] / (1j * omega) * Vtas / (q_dyn * self.jcl.general['A_ref']) + dCl_pitch = Qhb[2, 3, :len(k_red)] / (q_dyn * self.jcl.general['A_ref']) + # ToDo: calculate dCl from Qhh as cross-check + + # Do the actual plotting + fig = plt.figure() + ax = fig.add_axes([0.15, 0.15, 0.75, 0.75]) + ax.plot(k_red, np.abs(dCl_heave), '--x', label='dCl from heave') + ax.plot(k_red, np.abs(dCl_pitch), '--x', label='dCl from pitch') + fig.suptitle(f'{trimcase['desc']}', fontsize=16) + self.make_as_nice([ax]) + ax.set_xlabel('$k_{red}$') + ax.set_ylabel('[-]') + pp.savefig() + plt.close() + + # Step 3: plot each mode separately + fig, ax = plt.subplots(7, sharex=True, figsize=(8, 10)) + for i in range(Pb_pulse.shape[1]): + fig.suptitle(f'{trimcase['desc']}, Mode {i + 2}', fontsize=16) + for a in ax: + a.cla() + ax[0].plot(t, pulse_signal[i, :], '-', label='Pulse') + ax[1].plot(t, Pb_pulse[0, i, :], '-', label='Fx') + ax[2].plot(t, Pb_pulse[1, i, :], '-', label='Fy') + ax[3].plot(t, Pb_pulse[2, i, :], '-', label='Fz') + ax[4].plot(t, Pb_pulse[3, i, :], '-', label='Mx') + ax[5].plot(t, Pb_pulse[4, i, :], '-', label='My') + ax[6].plot(t, Pb_pulse[5, i, :], '-', label='Mz') + self.make_as_nice(ax) + ax[-1].set_xlabel('Time [s]') + ax[0].set_ylabel('[-]') + ax[1].set_ylabel('[N]') + ax[2].set_ylabel('[N]') + ax[3].set_ylabel('[N]') + ax[4].set_ylabel('[Nm]') + ax[5].set_ylabel('[Nm]') + ax[6].set_ylabel('[Nm]') + fig.tight_layout(pad=1.0) + pp.savefig() + + # Step 4: plot gust + fig, ax = plt.subplots(7, sharex=True, figsize=(8, 10)) + fig.suptitle(f'{trimcase['desc']}, Gust', fontsize=16) + ax[0].plot(t, gust_signal, '-', label='Gust') + ax[1].plot(t, Pb_gust[0, :], '-', label='Fx') + ax[2].plot(t, Pb_gust[1, :], '-', label='Fy') + ax[3].plot(t, Pb_gust[2, :], '-', label='Fz') + ax[4].plot(t, Pb_gust[3, :], '-', label='Mx') + ax[5].plot(t, Pb_gust[4, :], '-', label='My') + ax[6].plot(t, Pb_gust[5, :], '-', label='Mz') + self.make_as_nice(ax) + ax[-1].set_xlabel('Time [s]') + ax[0].set_ylabel('[-]') + ax[1].set_ylabel('[N]') + ax[2].set_ylabel('[N]') + ax[3].set_ylabel('[N]') + ax[4].set_ylabel('[Nm]') + ax[5].set_ylabel('[Nm]') + ax[6].set_ylabel('[Nm]') + fig.tight_layout(pad=1.0) + pp.savefig() + pp.close() + plt.close() diff --git a/loadskernel/post_processing.py b/loadskernel/post_processing.py index 61e57639..e8fc54d9 100755 --- a/loadskernel/post_processing.py +++ b/loadskernel/post_processing.py @@ -42,7 +42,7 @@ def __init__(self, jcl, model, trimcase, response): self.PHIcfd_strc = load_hdf5_sparse_matrix(self.model['PHIcfd_strc']) def force_summation_method(self): - logging.info('calculating forces & moments on structural set (force summation method)...') + logging.info('Calculating forces & moments on structural set (force summation method)...') response = self.response response['Pg_iner'] = np.zeros((len(response['t']), 6 * self.strcgrid['n'])) @@ -86,10 +86,8 @@ def force_summation_method(self): response['d2Ug_dt2'][i_step, :] = d2Ug_dt2_r + d2Ug_dt2_f def modal_displacement_method(self): - logging.info( - 'calculating forces & moments on structural set (modal displacement method)...') - logging.warning( - 'using the modal displacement method is not recommended, use force summation method instead.') + logging.info('Calculating forces & moments on structural set (modal displacement method)...') + logging.warning('Using the modal displacement method is not recommended, use force summation method instead.') response = self.response Kgg = load_hdf5_sparse_matrix(self.model['KGG']) @@ -101,7 +99,7 @@ def modal_displacement_method(self): response['Pg'][i_step, :] = Kgg.dot(Ug_f_body) def euler_transformation(self): - logging.info('apply euler angles...') + logging.info('Apply euler angles...') response = self.response response['Pg_iner_global'] = np.zeros((len(response['t']), 6 * self.strcgrid['n'])) @@ -156,7 +154,7 @@ def euler_transformation(self): response['Ug'][i_step, :] = response['Ug_r'][i_step, :] + response['Ug_f'][i_step, :] def cuttingforces(self): - logging.info('calculating cutting forces & moments...') + logging.info('Calculating cutting forces & moments...') response = self.response response['Pmon_local'] = np.zeros((len(response['t']), 6 * self.mongrid['n'])) for i_step in range(len(response['t'])): diff --git a/loadskernel/program_flow.py b/loadskernel/program_flow.py index c7dcd865..75bcf268 100755 --- a/loadskernel/program_flow.py +++ b/loadskernel/program_flow.py @@ -36,7 +36,6 @@ def __init__(self, self.post = post # True/False # debug options self.debug = False # True/False - self.restart = False # True/False # advanced options self.test = test # True/False # job control options @@ -128,25 +127,25 @@ def create_logfile_and_console_output(self, filename): # add the handler(s) to the root logger logger.addHandler(console) - logger.info('This is the log for process {}.'.format(self.myid)) + logger.info('This is the log for process %s.', self.myid) class Kernel(ProgramFlowHelper): def run(self): self.setup_logger() - logging.info('Starting Loads Kernel with job: ' + self.job_name) - logging.info('User ' + getpass.getuser() + ' on ' + platform.node() + ' (' + platform.platform() + ')') + logging.info('Starting Loads Kernel with job: %s', self.job_name) + logging.info('User %s on %s (%s)', getpass.getuser(), platform.node(), platform.platform()) logging.info('Software versions:') - logging.info(' - Loads Kernel ' + metadata.version('loadskernel') + ' (' + repr(loadskernel) + ')') - logging.info(' - Panel Aero ' + metadata.version('panelaero')) - logging.info(' - Python ' + platform.python_version()) - logging.info(' - Numpy ' + metadata.version('numpy')) - logging.info(' - Scipy ' + metadata.version('scipy')) - logging.info('pre: ' + str(self.pre)) - logging.info('main: ' + str(self.main)) - logging.info('post: ' + str(self.post)) - logging.info('test: ' + str(self.test)) + logging.info(' - Loads Kernel %s (%s)', metadata.version('loadskernel'), repr(loadskernel)) + logging.info(' - Panel Aero %s', metadata.version('panelaero')) + logging.info(' - Python %s', platform.python_version()) + logging.info(' - Numpy %s', metadata.version('numpy')) + logging.info(' - Scipy %s', metadata.version('scipy')) + logging.info('pre: %s', self.pre) + logging.info('main: %s', self.main) + logging.info('post: %s', self.post) + logging.info('test: %s', self.test) self.jcl = data_handling.load_jcl(self.job_name, self.path_input, self.jcl) # add machinefile to jcl self.jcl.machinefile = self.machinefile @@ -166,79 +165,86 @@ def run(self): self.print_logo() def run_pre(self): - logging.info('--> Starting preprocessing.') + logging.info('Starting preprocessing.') t_start = time.time() model = model_modul.Model(self.jcl, self.path_output) model.build_model() - logging.info('--> Saving model data.') + logging.info('Saving model data.') data_handling.dump_hdf5(self.path_output + 'model_' + self.job_name + '.hdf5', model.__dict__) - logging.info('--> Done in {}.'.format(seconds2string(time.time() - t_start))) + logging.info('Done in %s.', seconds2string(time.time() - t_start)) - def main_common(self, model, jcl, i): + def main_core(self, model, jcl, i): + """ + This functions performs the main processing for one trimcase / subcase. + It can be called either from various run_main_... functions. + """ logging.info('') logging.info('========================================') - logging.info('trimcase: ' + jcl.trimcase[i]['desc']) - logging.info('subcase: ' + str(jcl.trimcase[i]['subcase'])) - logging.info('(case ' + str(i + 1) + ' of ' + str(len(jcl.trimcase)) + ')') + logging.info('trimcase: %s', jcl.trimcase[i]['desc']) + logging.info('subcase: %s', jcl.trimcase[i]['subcase']) + logging.info('(case %d of %d)', i + 1, len(jcl.trimcase)) logging.info('========================================') solution_i = solution_sequences.SolutionSequences(model, jcl, jcl.trimcase[i], jcl.simcase[i]) + # Set trim conditions and execute trim. solution_i.set_trimcond() solution_i.exec_trim() - # solution_i.iterative_trim() + # In case of successful trim, do time simulation if requested in the simcase. if solution_i.successful and 't_final' and 'dt' in jcl.simcase[i].keys(): solution_i.exec_sim() - elif solution_i.successful and 'flutter' in jcl.simcase[i] and jcl.simcase[i]['flutter']: + # The post processing applies only to trim and time/frequency simulations. Thus, an alternative place for + # post processing could be in the solution_sequence.py. Not sure which place is better... + # Also, the name 'post_processing' might be misleading here, as it is not post processing of the entire + # job (post=True), but only of the trim / sim solution sequence. + if solution_i.successful: + post_processing_i = post_processing.PostProcessing(jcl, model, jcl.trimcase[i], solution_i.response) + post_processing_i.force_summation_method() + post_processing_i.euler_transformation() + post_processing_i.cuttingforces() + del post_processing_i + # Look if any other special analyses are requested (such as flutter, derivatives, GAFs) in the simcase. + if 'flutter' in jcl.simcase[i] and jcl.simcase[i]['flutter']: solution_i.exec_flutter() elif solution_i.successful and 'derivatives' in jcl.simcase[i] and jcl.simcase[i]['derivatives']: solution_i.calc_jacobian() solution_i.calc_derivatives() + elif solution_i.successful and 'gaf' in jcl.simcase[i] and jcl.simcase[i]['gaf']: + solution_i.calc_gafs() + # Collect response from solution sequence, then destroy it to free memory. response = solution_i.response response['i'] = i response['successful'] = solution_i.successful del solution_i - if response['successful']: - post_processing_i = post_processing.PostProcessing(jcl, model, jcl.trimcase[i], response) - post_processing_i.force_summation_method() - post_processing_i.euler_transformation() - post_processing_i.cuttingforces() - del post_processing_i return response def run_main_sequential(self): - logging.info('--> Starting Main in sequential mode for {} trimcase(s).'.format(len(self.jcl.trimcase))) + logging.info('Starting Main in sequential mode for %d trimcase(s).', len(self.jcl.trimcase)) t_start = time.time() model = data_handling.load_hdf5(self.path_output + 'model_' + self.job_name + '.hdf5') - if self.myid == 0: - mon = gather_loads.GatherLoads(self.jcl, model) - if self.restart: - logging.info('Restart option: loading existing responses.') - # open response - responses = data_handling.load_hdf5_responses(self.job_name, self.path_output) - fid = data_handling.open_hdf5(self.path_output + 'response_' + self.job_name + '.hdf5') # open response - - for i in range(len(self.jcl.trimcase)): - if self.restart and i in [response['i'][()] for response in responses]: - logging.info('Restart option: found existing response.') - response = responses[[response['i'][()] for response in responses].index(i)] - else: - jcl = copy.deepcopy(self.jcl) - response = self.main_common(model, jcl, i) - if self.myid == 0 and response['successful']: - mon.gather_monstations(self.jcl.trimcase[i], response) - mon.gather_dyn2stat(response) - logging.info('--> Saving response(s).') + # Initialize module to gather loads + loads = gather_loads.GatherLoads(self.jcl, model) + # Open response + fid = data_handling.open_hdf5(self.path_output + 'response_' + self.job_name + '.hdf5') + # Loop over all trimcases + for i, trimcase in enumerate(self.jcl.trimcase): + jcl = copy.deepcopy(self.jcl) + # Perfome the actual calculations for this trimcase + response = self.main_core(model, jcl, i) + if response['successful']: + # Gathers loads into monstations and dyn2stat data structures + loads.gather_monstations(trimcase, response) + loads.gather_dyn2stat(response) + # Save response to file + logging.info('Saving response(s).') data_handling.write_hdf5(fid, response, path='/' + str(response['i'])) - if self.myid == 0: - # close response - data_handling.close_hdf5(fid) - - logging.info('--> Saving monstation(s).') - data_handling.dump_hdf5(self.path_output + 'monstations_' + self.job_name + '.hdf5', mon.monstations) - - logging.info('--> Saving dyn2stat.') - data_handling.dump_hdf5(self.path_output + 'dyn2stat_' + self.job_name + '.hdf5', mon.dyn2stat) - logging.info('--> Done in {}.'.format(seconds2string(time.time() - t_start))) + # Close response + data_handling.close_hdf5(fid) + # Save gathered loads, even if the data structures might be empty (e.g. if flutter analysis was requested). + logging.info('Saving monstation(s).') + data_handling.dump_hdf5(self.path_output + 'monstations_' + self.job_name + '.hdf5', loads.monstations) + logging.info('Saving dyn2stat.') + data_handling.dump_hdf5(self.path_output + 'dyn2stat_' + self.job_name + '.hdf5', loads.dyn2stat) + logging.info('Done in %s.', seconds2string(time.time() - t_start)) def run_main_multiprocessing(self): """ @@ -248,11 +254,9 @@ def run_main_multiprocessing(self): This concept is adapted from Jörg Bornschein (see https://github.com/jbornschein/mpi4py-examples/blob/master/ 09-task-pull.py) """ - logging.info( - '--> Starting Main in multiprocessing mode for %d trimcase(s).', len(self.jcl.trimcase)) + logging.info('Starting Main in multiprocessing mode for %d trimcase(s).', len(self.jcl.trimcase)) t_start = time.time() - model = data_handling.load_hdf5( - self.path_output + 'model_' + self.job_name + '.hdf5') + model = data_handling.load_hdf5(self.path_output + 'model_' + self.job_name + '.hdf5') # MPI tags can be any integer values tags = {'ready': 0, 'start': 1, @@ -262,9 +266,9 @@ def run_main_multiprocessing(self): # The master process runs on the first processor if self.myid == 0: n_workers = self.comm.Get_size() - 1 - logging.info('--> I am the master with %d worker(s).', n_workers) + logging.info('I am the master with %d worker(s).', n_workers) - mon = gather_loads.GatherLoads(self.jcl, model) + loads = gather_loads.GatherLoads(self.jcl, model) # open response fid = data_handling.open_hdf5(self.path_output + 'response_' + self.job_name + '.hdf5') @@ -281,7 +285,7 @@ def run_main_multiprocessing(self): # Worker is ready, send out a new subcase. if i_subcase < len(self.jcl.trimcase): self.comm.send(i_subcase, dest=source, tag=tags['start']) - logging.info('--> Sending case %d of %d to worker %d', i_subcase + 1, len(self.jcl.trimcase), source) + logging.info('Sending case %d of %d to worker %d', i_subcase + 1, len(self.jcl.trimcase), source) i_subcase += 1 else: # No more task to do, send the exit signal. @@ -291,13 +295,13 @@ def run_main_multiprocessing(self): # The worker has returned a response. response = data if response['successful']: - logging.info("--> Received response ('successful') from worker %d.", source) - mon.gather_monstations(self.jcl.trimcase[response['i']], response) - mon.gather_dyn2stat(response) + logging.info("Received response ('successful') from worker %d.", source) + loads.gather_monstations(self.jcl.trimcase[response['i']], response) + loads.gather_dyn2stat(response) else: # Trim failed, no post processing, save the empty response - logging.info("--> Received response ('failed') from worker %d.", source) - logging.info('--> Saving response(s).') + logging.info("Received response ('failed') from worker %d.", source) + logging.info('Saving response(s).') data_handling.write_hdf5(fid, response, path='/' + str(response['i'])) elif tag == tags['exit']: @@ -306,13 +310,11 @@ def run_main_multiprocessing(self): closed_workers += 1 # close response data_handling.close_hdf5(fid) - logging.info('--> Saving monstation(s).') - data_handling.dump_hdf5(self.path_output + 'monstations_' + self.job_name + '.hdf5', - mon.monstations) + logging.info('Saving monstation(s).') + data_handling.dump_hdf5(self.path_output + 'monstations_' + self.job_name + '.hdf5', loads.monstations) - logging.info('--> Saving dyn2stat.') - data_handling.dump_hdf5(self.path_output + 'dyn2stat_' + self.job_name + '.hdf5', - mon.dyn2stat) + logging.info('Saving dyn2stat.') + data_handling.dump_hdf5(self.path_output + 'dyn2stat_' + self.job_name + '.hdf5', loads.dyn2stat) # The worker process runs on all other processors else: logging.info('I am worker on process %d.', self.myid) @@ -324,7 +326,7 @@ def run_main_multiprocessing(self): if tag == tags['start']: # Start a new job - response = self.main_common(model, self.jcl, i_subcase) + response = self.main_core(model, self.jcl, i_subcase) self.comm.send(response, dest=0, tag=tags['done']) elif tag == tags['exit']: # Received an exit signal. @@ -332,17 +334,17 @@ def run_main_multiprocessing(self): # Confirm the exit signal. self.comm.send(None, dest=0, tag=tags['exit']) - logging.info('--> Done in {}.'.format(seconds2string(time.time() - t_start))) + logging.info('Done in %s.', seconds2string(time.time() - t_start)) def run_post(self): model = data_handling.load_hdf5(self.path_output + 'model_' + self.job_name + '.hdf5') responses = data_handling.load_hdf5_responses(self.job_name, self.path_output) - logging.info('--> Loading monstations(s).') + logging.info('Loading monstations(s).') monstations = data_handling.load_hdf5(self.path_output + 'monstations_' + self.job_name + '.hdf5') - logging.info('--> Loading dyn2stat.') + logging.info('Loading dyn2stat.') dyn2stat_data = data_handling.load_hdf5(self.path_output + 'dyn2stat_' + self.job_name + '.hdf5') - logging.info('--> Drawing some standard plots.') + logging.info('Drawing some standard plots.') if 'flutter' in self.jcl.simcase[0] and self.jcl.simcase[0]['flutter']: plt = plotting_standard.FlutterPlots(self.jcl, model) plt.add_responses(responses) @@ -355,8 +357,11 @@ def run_post(self): elif 'limit_turbulence' in self.jcl.simcase[0] and self.jcl.simcase[0]['limit_turbulence']: plt = plotting_standard.TurbulencePlots(self.jcl, model) plt.add_monstations(monstations) - plt.plot_monstations( - self.path_output + 'monstations_turbulence_' + self.job_name + '.pdf') + plt.plot_monstations(self.path_output + 'monstations_turbulence_' + self.job_name + '.pdf') + elif 'gaf' in self.jcl.simcase[0] and self.jcl.simcase[0]['gaf']: + plt = plotting_standard.GAFPlots(self.jcl, model) + plt.add_responses(responses) + plt.plot_pluse(self.path_output + 'pulse_' + self.job_name + '.pdf') else: # Here come the loads plots plt = plotting_standard.LoadPlots(self.jcl, model) @@ -367,7 +372,7 @@ def run_post(self): 'cs_signal', 'controller']]): plt.plot_monstations_time(self.path_output + 'monstations_time_' + self.job_name + '.pdf') # nur sim - logging.info('--> Saving auxiliary output data.') + logging.info('Saving auxiliary output data.') aux_out = auxiliary_output.AuxiliaryOutput(self.jcl, model, self.jcl.trimcase) aux_out.crit_trimcases = plt.crit_trimcases aux_out.dyn2stat_data = dyn2stat_data @@ -395,13 +400,13 @@ def run_test(self): """ # Import plotting_extra not before here, as the import of graphical libraries such as mayavi takes a long time and # fails of systems without graphical display (such as HPS clusters). - from loadskernel import plotting_extra + from loadskernel import plotting_extra # pylint: disable=import-outside-toplevel # Load the model and the response as usual model = data_handling.load_hdf5(self.path_output + 'model_' + self.job_name + '.hdf5') responses = data_handling.load_hdf5_responses(self.job_name, self.path_output) - logging.info('--> Drawing some more detailed plots.') + logging.info('Drawing some more detailed plots.') plt = plotting_extra.DetailedPlots(self.jcl, model) plt.add_responses(responses) if 't_final' and 'dt' in self.jcl.simcase[0].keys(): @@ -449,7 +454,7 @@ class ClusterMode(Kernel): k = program_flow.ClusterMode('jcl_name', ...) k.run_cluster(sys.argv[2]) or - k.gather_cluster() + k.gather_responses() """ def run_cluster(self, i): @@ -458,8 +463,8 @@ def run_cluster(self, i): # add machinefile to jcl self.jcl.machinefile = self.machinefile self.setup_logger_cluster(i=i) - logging.info('Starting Loads Kernel with job: ' + self.job_name) - logging.info('User ' + getpass.getuser() + ' on ' + platform.node() + ' (' + platform.platform() + ')') + logging.info('Starting Loads Kernel with job: %s', self.job_name) + logging.info('User %s on %s (%s)', getpass.getuser(), platform.node(), platform.platform()) logging.info('Cluster array mode') self.run_main_single(i) @@ -471,7 +476,7 @@ def run_main_single(self, i): """ This function calculates one single load case, e.g. using CFD with mpi hosts on a cluster. """ - logging.info('--> Starting main in single mode for {} trimcase(s).'.format(len(self.jcl.trimcase))) + logging.info('Starting main in single mode for %d trimcase(s).', len(self.jcl.trimcase)) t_start = time.time() model = data_handling.load_hdf5(self.path_output + 'model_' + self.job_name + '.hdf5') jcl = copy.deepcopy(self.jcl) @@ -490,47 +495,44 @@ def run_main_single(self, i): + str(self.jcl.trimcase[i]['subcase']) + '.pickle', 'wb') as f: data_handling.dump_pickle(empty_response, f) # Start the simulation - response = self.main_common(model, jcl, i) + response = self.main_core(model, jcl, i) # Overwrite the empty response from above if self.myid == 0: - logging.info('--> Saving response(s).') + logging.info('Saving response(s).') path_responses = data_handling.check_path(self.path_output + 'responses/') with open(path_responses + 'response_' + self.job_name + '_subcase_' + str(self.jcl.trimcase[i]['subcase']) + '.pickle', 'wb') as f: data_handling.dump_pickle(response, f) - logging.info('--> Done in {}.'.format(seconds2string(time.time() - t_start))) + logging.info('Done in %s.', seconds2string(time.time() - t_start)) - def gather_cluster(self): + def gather_responses(self): self.setup_logger() - t_start = time.time() - logging.info('Starting Loads Kernel with job: ' + self.job_name) - logging.info('user ' + getpass.getuser() + ' on ' + platform.node() + ' (' + platform.platform() + ')') - logging.info('cluster gather mode') + logging.info('Starting Loads Kernel with job: %s', self.job_name) + logging.info('User %s on %s (%s)', getpass.getuser(), platform.node(), platform.platform()) + logging.info('response gather mode') self.jcl = data_handling.load_jcl(self.job_name, self.path_input, self.jcl) model = data_handling.load_hdf5(self.path_output + 'model_' + self.job_name + '.hdf5') responses = data_handling.gather_responses(self.job_name, data_handling.check_path(self.path_output + 'responses')) mon = gather_loads.GatherLoads(self.jcl, model) fid = data_handling.open_hdf5(self.path_output + 'response_' + self.job_name + '.hdf5') # open response - for i in range(len(self.jcl.trimcase)): + for i, (trimcase) in enumerate(self.jcl.trimcase): response = responses[[response['i'] for response in responses].index(i)] if response['successful']: - mon.gather_monstations(self.jcl.trimcase[i], response) + mon.gather_monstations(trimcase, response) mon.gather_dyn2stat(response) - logging.info('--> Saving response(s).') + logging.info('Saving response(s).') data_handling.write_hdf5(fid, response, path='/' + str(response['i'])) # close response data_handling.close_hdf5(fid) - logging.info('--> Saving monstation(s).') + logging.info('Saving monstation(s).') data_handling.dump_hdf5(self.path_output + 'monstations_' + self.job_name + '.hdf5', mon.monstations) - logging.info('--> Saving dyn2stat.') + logging.info('Saving dyn2stat.') data_handling.dump_hdf5(self.path_output + 'dyn2stat_' + self.job_name + '.hdf5', mon.dyn2stat) - logging.info( - '--> Done in {}.'.format(seconds2string(time.time() - t_start))) logging.info('Loads Kernel finished.') self.print_logo() @@ -551,7 +553,7 @@ def str2bool(v): def seconds2string(seconds): m, s = divmod(seconds, 60) h, m = divmod(m, 60) - return '{:n}:{:02n}:{:02n} [h:mm:ss]'.format(h, m, round(s)) + return f'{h:n}:{m:02n}:{round(s):02n} [h:mm:ss]' def command_line_interface(): diff --git a/loadskernel/solution_sequences.py b/loadskernel/solution_sequences.py index 4683c699..08840e4b 100755 --- a/loadskernel/solution_sequences.py +++ b/loadskernel/solution_sequences.py @@ -3,25 +3,17 @@ import numpy as np import scipy.optimize as so from scipy.integrate import ode +from scipy.fftpack import fft, fftfreq from loadskernel.integrate import RungeKutta4, ExplicitEuler, AdamsBashforth -from loadskernel.equations.steady import Steady -from loadskernel.equations.cfd import CfdSteady, CfdUnsteady -from loadskernel.equations.nonlin_steady import NonlinSteady -from loadskernel.equations.unsteady import Unsteady -from loadskernel.equations.landing import Landing +from loadskernel.equations.mona_time_domain import Steady, Unsteady, NonlinSteady, Landing +from loadskernel.equations.cfd_time_domain import CfdSteady, CfdUnsteady from loadskernel.equations.common import ConvergenceError -from loadskernel.equations.frequency_domain import GustExcitation -from loadskernel.equations.frequency_domain import TurbulenceExcitation, LimitTurbulence -from loadskernel.equations.frequency_domain import KMethod -from loadskernel.equations.frequency_domain import KEMethod -from loadskernel.equations.frequency_domain import PKMethodRodden -from loadskernel.equations.frequency_domain import PKMethodSchwochow -from loadskernel.equations.state_space import StateSpaceAnalysis -from loadskernel.equations.state_space import JacobiAnalysis +from loadskernel.equations import mona_frequency_domain, cfd_frequency_domain, mona_state_space from loadskernel.trim_conditions import TrimConditions from loadskernel.cfd_interfaces.tau_interface import TauError -from loadskernel.io_functions.data_handling import load_hdf5_dict +from loadskernel.io_functions.data_handling import load_hdf5_sparse_matrix, load_hdf5_dict +from loadskernel.solution_tools import polynomial_pulse, one_m_cosine_pulse class SolutionSequences(TrimConditions): @@ -54,12 +46,11 @@ def calc_jacobian(self): if self.jcl.aero['method'] in ['mona_steady']: equations = Steady(self) else: - logging.error('Unknown aero method: ' + str(self.jcl.aero['method'])) + logging.error('Unknown aero method: %s', self.jcl.aero['method']) + equations = None - # flight condition X0 = self.response['X'][0, :] - # X0 = np.array(self.trimcond_X[:,2], dtype='float') - logging.info('Calculating jacobian for ' + str(len(X0)) + ' variables...') + logging.info('Calculating jacobian for %d variables...', len(X0)) # epsilon sollte klein sein, dt sollte 1.0s sein jac = self.approx_jacobian(X0=X0, func=equations.equations, epsilon=0.01, dt=1.0) self.response['X0'] = X0 # Linearisierungspunkt @@ -89,12 +80,10 @@ def calc_jacobian(self): self.response['desc'] = self.trimcase['desc'] # perform analysis on jacobian matrix - equations = JacobiAnalysis(self.response) + equations = mona_state_space.JacobiAnalysis(self.response) equations.eval_equations() def calc_derivatives(self): - self.macgrid = load_hdf5_dict(self.model['macgrid']) - self.calc_flexible_derivatives() self.calc_rigid_derivatives() self.calc_additional_derivatives('rigid') @@ -106,12 +95,13 @@ def calc_derivatives(self): logging.info('--------------------------------------------------------------------------------------') def calc_rigid_derivatives(self): - if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'hybrid']: + if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady']: equations = Steady(self) elif self.jcl.aero['method'] in ['nonlin_steady']: equations = NonlinSteady(self) else: - logging.error('Unknown aero method: ' + str(self.jcl.aero['method'])) + logging.error('Unknown aero method: %s', self.jcl.aero['method']) + equations = None A = self.jcl.general['A_ref'] delta = 0.01 @@ -125,8 +115,8 @@ def calc_rigid_derivatives(self): xi[i] += delta response = equations.equations(xi, 0.0, 'trim_full_output') Pmac_c = (response['Pmac'] - response0['Pmac']) / response['q_dyn'] / A / delta - derivatives.append([Pmac_c[0], Pmac_c[1], Pmac_c[2], Pmac_c[3] / self.macgrid['b_ref'], - Pmac_c[4] / self.macgrid['c_ref'], Pmac_c[5] / self.macgrid['b_ref']]) + derivatives.append([Pmac_c[0], Pmac_c[1], Pmac_c[2], Pmac_c[3] / self.model['macgrid']['b_ref'], + Pmac_c[4] / self.model['macgrid']['c_ref'], Pmac_c[5] / self.model['macgrid']['b_ref']]) # write back original response and store results self.response['rigid_parameters'] = self.trimcond_X[:, 0].tolist() self.response['rigid_derivatives'] = derivatives @@ -165,8 +155,8 @@ def calc_flexible_derivatives(self): # re-calculate new trim self.exec_trim() Pmac_c = (self.response['Pmac'] - response0['Pmac']) / response0['q_dyn'] / A / delta - derivatives.append([Pmac_c[0, 0], Pmac_c[0, 1], Pmac_c[0, 2], Pmac_c[0, 3] / self.macgrid['b_ref'], - Pmac_c[0, 4] / self.macgrid['c_ref'], Pmac_c[0, 5] / self.macgrid['b_ref']]) + derivatives.append([Pmac_c[0, 0], Pmac_c[0, 1], Pmac_c[0, 2], Pmac_c[0, 3] / self.model['macgrid']['b_ref'], + Pmac_c[0, 4] / self.model['macgrid']['c_ref'], Pmac_c[0, 5] / self.model['macgrid']['b_ref']]) # restore trim condition for next loop self.trimcond_X = copy.deepcopy(trimcond_X0) # write back original response and store results @@ -177,14 +167,13 @@ def calc_flexible_derivatives(self): def calc_NP(self): pos = self.response['flexible_parameters'].index('theta') self.response['NP_flex'] = np.zeros(3) - self.response['NP_flex'][0] = self.macgrid['offset'][0, 0] - self.jcl.general['c_ref'] \ + self.response['NP_flex'][0] = self.model['macgrid']['offset'][0, 0] - self.jcl.general['c_ref'] \ * self.response['flexible_derivatives'][pos][4] / self.response['flexible_derivatives'][pos][2] - self.response['NP_flex'][1] = self.macgrid['offset'][0, 1] + self.jcl.general['b_ref'] \ + self.response['NP_flex'][1] = self.model['macgrid']['offset'][0, 1] + self.jcl.general['b_ref'] \ * self.response['flexible_derivatives'][pos][3] / self.response['flexible_derivatives'][pos][2] logging.info('--------------------------------------------------------------------------------------') logging.info('Aeroelastic neutral point / aerodynamic center:') - logging.info('NP_flex (x,y) = {:0.4g},{:0.4g}'.format( - self.response['NP_flex'][0], self.response['NP_flex'][1])) + logging.info('NP_flex (x,y) = %0.4g,%0.4g', self.response['NP_flex'][0], self.response['NP_flex'][1]) def calc_cs_effectiveness(self): logging.info('--------------------------------------------------------------------------------------') @@ -195,8 +184,7 @@ def calc_cs_effectiveness(self): pos_flex = self.response['flexible_parameters'].index(p) d = np.array(self.response['flexible_derivatives'][pos_flex]) \ / np.array(self.response['rigid_derivatives'][pos_rigid]) - tmp = '{:>20} {:< 10.4g} {:< 10.4g} {:< 10.4g} {:< 10.4g} {:< 10.4g} {:< 10.4g}'.format( - p, d[0], d[1], d[2], d[3], d[4], d[5]) + tmp = f'{p:>20} {d[0]:< 10.4g} {d[1]:< 10.4g} {d[2]:< 10.4g} {d[3]:< 10.4g} {d[4]:< 10.4g} {d[5]:< 10.4g}' logging.info(tmp) def calc_additional_derivatives(self, key): @@ -227,26 +215,36 @@ def print_derivatives(self, key): parameters = self.response[key + '_parameters'] derivatives = self.response[key + '_derivatives'] logging.info('--------------------------------------------------------------------------------------') - logging.info('Calculated ' + key + ' derivatives for ' + str(len(parameters)) + ' variables.') - logging.info('MAC_ref = {}'.format(self.jcl.general['MAC_ref'])) - logging.info('A_ref = {}'.format(self.jcl.general['A_ref'])) - logging.info('b_ref = {}'.format(self.jcl.general['b_ref'])) - logging.info('c_ref = {}'.format(self.jcl.general['c_ref'])) - logging.info('q_dyn = {}'.format(self.response['q_dyn'][0])) + logging.info('Calculated %s derivatives for %d variables.', key, len(parameters)) + logging.info('MAC_ref = %s', self.jcl.general['MAC_ref']) + logging.info('A_ref = %s', self.jcl.general['A_ref']) + logging.info('b_ref = %s', self.jcl.general['b_ref']) + logging.info('c_ref = %s', self.jcl.general['c_ref']) + logging.info('q_dyn = %s', self.response['q_dyn'][0]) logging.info('Derivatives given in body axis (aft-right-up):') logging.info(' Cx Cy Cz Cmx Cmy Cmz') for p, d in zip(parameters, derivatives): - tmp = '{:>20} {:< 10.4g} {:< 10.4g} {:< 10.4g} {:< 10.4g} {:< 10.4g} {:< 10.4g}'.format( - p, d[0], d[1], d[2], d[3], d[4], d[5]) + tmp = f'{p:>20} {d[0]:< 10.4g} {d[1]:< 10.4g} {d[2]:< 10.4g} {d[3]:< 10.4g} {d[4]:< 10.4g} {d[5]:< 10.4g}' logging.info(tmp) def exec_trim(self): - if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'hybrid', 'nonlin_steady', 'freq_dom']: + if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'nonlin_steady', 'freq_dom', 'mona_freq_dom']: self.direct_trim() elif self.jcl.aero['method'] in ['cfd_steady', 'cfd_unsteady']: self.iterative_trim() + elif self.jcl.aero['method'] in ['cfd_freq_dom']: + logging.info('Using response / trim data form CFD-based GAF computation.') + # Fetch data from the GAF computation and fill the response. + # The assumption is that the gust is superposed with the linearization point used in the GAF computation. + key = '.'.join(self.trimcase['desc'].split('.')[:-1]) + if key in self.model['GAFs']: + self.response = load_hdf5_dict(self.model['GAFs'][key]['response']) + self.successful = True + else: + logging.error('No response / trim data found for "%s" in model!', key) + self.successful = False else: - logging.error('Unknown aero method: ' + str(self.jcl.aero['method'])) + logging.error('Unknown aero method: %s', str(self.jcl.aero['method'])) if self.successful: # To align the trim results with the time/frequency simulations, we expand the response by one dimension. # Notation: (n_timesteps, n_dof) --> the trim results can be considered as the solution at time step zero. @@ -262,27 +260,28 @@ def direct_trim(self): # ward-difference approximation. # http://www.math.utah.edu/software/minpack/minpack/hybrd.html - if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'hybrid', - 'freq_dom'] and not hasattr(self.jcl, 'landinggear'): + if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', + 'freq_dom', 'mona_freq_dom'] and not hasattr(self.jcl, 'landinggear'): equations = Steady(self) elif self.jcl.aero['method'] in ['nonlin_steady']: equations = NonlinSteady(self) elif self.simcase['landinggear'] and self.jcl.landinggear['method'] in ['generic', 'skid']: equations = Landing(self) else: - logging.error('Unknown aero method: ' + str(self.jcl.aero['method'])) + logging.error('Unknown aero method: %s', self.jcl.aero['method']) + equations = None xfree_0 = np.array(self.trimcond_X[:, 2], dtype='float')[np.where((self.trimcond_X[:, 1] == 'free'))[0]] if self.trimcase['maneuver'] == 'bypass': - logging.info('Running bypass...') + logging.info('Bypassing trim.') self.response = equations.eval_equations(xfree_0, time=0.0, modus='trim_full_output') self.successful = True else: - logging.info('Running trim for ' + str(len(xfree_0)) + ' variables...') + logging.info('Running trim for %d variables...', len(xfree_0)) xfree, info, status, msg = so.fsolve(equations.eval_equations, xfree_0, args=(0.0, 'trim'), full_output=True) - logging.info(msg) - logging.debug('Function evaluations: ' + str(info['nfev'])) + logging.info('%s', msg) + logging.debug('Function evaluations: %d', info['nfev']) # no errors, check trim status for success if status == 1: @@ -292,29 +291,30 @@ def direct_trim(self): else: self.response = {} self.successful = False - logging.warning('SolutionSequences failed for subcase {}. The SolutionSequences solver reports: {}'.format( - self.trimcase['subcase'], msg)) + logging.warning('SolutionSequences failed for subcase %s. The SolutionSequences solver reports: %s', + self.trimcase['subcase'], msg) equations.finalize() return def iterative_trim(self): - if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'hybrid']: + if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady']: equations = Steady(self) - elif self.jcl.aero['method'] in ['cfd_steady', 'cfd_unsteady']: + elif self.jcl.aero['method'] in ['cfd_steady', 'cfd_unsteady', 'cfd_freq_dom']: equations = CfdSteady(self) else: - logging.error('Unknown aero method: ' + str(self.jcl.aero['method'])) + logging.error('Unknown aero method: %s', self.jcl.aero['method']) + equations = None self.set_modal_states_fix() # start trim from scratch xfree_0 = np.array(self.trimcond_X[:, 2], dtype='float')[np.where((self.trimcond_X[:, 1] == 'free'))[0]] if self.trimcase['maneuver'] == 'bypass': - logging.info('running bypass...') + logging.info('Bypassing trim.') self.response = equations.eval_equations(xfree_0, time=0.0, modus='trim_full_output') self.successful = True else: - logging.info('Running trim for ' + str(len(xfree_0)) + ' variables...') + logging.info('Running trim for %d variables...', len(xfree_0)) """ Because the iterative trim is typically used in combination with CFD, some solver settings need to be modified. - The jacobian matrix is constructed using finite differences. With CFD, a sufficiently large step size should @@ -325,43 +325,42 @@ def iterative_trim(self): parameter 'factor=0.1'. """ try: - xfree, info, status, msg = so.fsolve(equations.eval_equations_iteratively, xfree_0, args=(0.0, 'trim'), - full_output=True, epsfcn=1.0e-3, xtol=1.0e-3, factor=0.1) + xfree, info, status, msg = so.fsolve( + equations.eval_equations_iteratively, xfree_0, args=(0.0, 'trim'), + full_output=True, epsfcn=1.0e-3, xtol=1.0e-3, factor=0.1) except TauError as e: self.response = {} self.successful = False - logging.warning('SolutionSequences failed for subcase {} due to CFDError: {}'.format( - self.trimcase['subcase'], e)) + logging.warning('SolutionSequences failed for subcase %s due to CFDError: %s', + self.trimcase['subcase'], e) except ConvergenceError as e: self.response = {} self.successful = False - logging.warning('SolutionSequences failed for subcase {} due to ConvergenceError: {}'.format( - self.trimcase['subcase'], e)) + logging.warning('SolutionSequences failed for subcase %s due to ConvergenceError: %s', + self.trimcase['subcase'], e) else: - logging.info(msg) - logging.info('function evaluations: ' + str(info['nfev'])) - # no errors, check trim status for success + logging.info('%s', msg) + logging.info('function evaluations: %d', info['nfev']) if status == 1: - # if trim was successful, then do one last evaluation with the final parameters. self.response = equations.eval_equations_iteratively( xfree, time=0.0, modus='trim_full_output') self.successful = True else: self.response = {} self.successful = False - logging.warning('SolutionSequences failed for subcase {}. The SolutionSequences solver reports: {}'.format( - self.trimcase['subcase'], msg)) + logging.warning('SolutionSequences failed for subcase %s. The SolutionSequences solver reports: %s', + self.trimcase['subcase'], msg) equations.finalize() return def exec_sim(self): # select solution sequence - if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'hybrid', 'nonlin_steady', 'cfd_unsteady']: + if self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady', 'nonlin_steady', 'cfd_unsteady']: self.exec_sim_time_dom() - elif self.jcl.aero['method'] in ['freq_dom']: + elif self.jcl.aero['method'] in ['freq_dom', 'mona_freq_dom', 'cfd_freq_dom']: self.exec_sim_freq_dom() else: - logging.error('Unknown aero method: ' + str(self.jcl.aero['method'])) + logging.error('Unknown aero method: %s', str(self.jcl.aero['method'])) def exec_sim_time_dom(self): """ @@ -371,7 +370,7 @@ def exec_sim_time_dom(self): # get initial solution from trim X0 = self.response['X'][0, :] # select solution sequence - if self.jcl.aero['method'] in ['mona_steady', 'hybrid'] and not hasattr(self.jcl, 'landinggear'): + if self.jcl.aero['method'] in ['mona_steady'] and not hasattr(self.jcl, 'landinggear'): equations = Steady(self, X0) elif self.jcl.aero['method'] in ['nonlin_steady']: equations = NonlinSteady(self, X0) @@ -383,7 +382,7 @@ def exec_sim_time_dom(self): equations = Landing(self, X0) elif self.jcl.aero['method'] in ['mona_unsteady']: if 'disturbance' in self.simcase.keys(): - logging.info('Adding disturbance of {} to state(s) '.format(self.simcase['disturbance'])) + logging.info('Adding disturbance of %s to state(s) ', self.simcase['disturbance']) self.response['X'][0, 11 + self.simcase['disturbance_mode']] += self.simcase['disturbance'] # add lag states to system self.add_lagstates() @@ -393,7 +392,8 @@ def exec_sim_time_dom(self): elif self.jcl.aero['method'] in ['cfd_unsteady']: equations = CfdUnsteady(self, X0) else: - logging.error('Unknown aero method: ' + str(self.jcl.aero['method'])) + logging.error('Unknown aero method: %s', self.jcl.aero['method']) + equations = None """ There are two ways of time intergartion. @@ -424,7 +424,7 @@ def exec_sim_time_dom(self): xt = [] t = [] - logging.info('Running time simulation for ' + str(t_final) + ' sec...') + logging.info('Running time simulation for %g sec...', t_final) if self.jcl.aero['method'] in ['cfd_unsteady']: integrator = self.select_integrator(equations, 'AdamsBashforth_FixedTimeStep', dt_integration) integrator.set_initial_value(X0, 0.0) @@ -485,42 +485,245 @@ def select_integrator(self, equations, integration_scheme='AdamsBashforth', step elif integration_scheme == 'RK45': integrator = ode(equations.ode_arg_sorter).set_integrator('dopri5', nsteps=2000, rtol=1e-2, atol=1e-8, max_step=1e-4) + else: + logging.error('Unknown integration scheme: %s.', integration_scheme) + integrator = None return integrator def exec_sim_freq_dom(self): # get initial solution from trim X0 = self.response['X'][0, :] - # select solution sequence - if self.simcase['gust']: - equations = GustExcitation(self, X0) - elif self.simcase['turbulence']: - equations = TurbulenceExcitation(self, X0) - elif self.simcase['limit_turbulence']: - equations = LimitTurbulence(self, X0) - self.response['Pmon_turb'] = 0.0 - self.response['correlations'] = 0.0 + if self.jcl.aero['method'] in ['freq_dom', 'mona_freq_dom']: + # select solution sequence + if self.simcase['gust']: + equations = mona_frequency_domain.GustExcitation(self, X0) + elif self.simcase['turbulence']: + equations = mona_frequency_domain.TurbulenceExcitation(self, X0) + elif self.simcase['limit_turbulence']: + equations = mona_frequency_domain.LimitTurbulence(self, X0) + self.response['Pmon_turb'] = 0.0 + self.response['correlations'] = 0.0 + else: + logging.error('Unknown frequency domain simulation type.') + equations = None + elif self.jcl.aero['method'] in ['cfd_freq_dom']: + # select solution sequence + if self.simcase['gust']: + equations = cfd_frequency_domain.GustExcitation(self, X0) + else: + logging.error('Unknown CFD-based simulation type.') + equations = None + else: + logging.error('Unknown aero method: %s', self.jcl.aero['method']) + equations = None response_sim = equations.eval_equations() - for key in response_sim.keys(): - self.response[key] = response_sim[key] + self.response[key] + for key, item in response_sim.items(): + self.response[key] = item + self.response[key] logging.info('Frequency domain simulation finished.') self.successful = True def exec_flutter(self): - # get initial solution from trim - X0 = self.response['X'][0, :] - # select solution sequence - if self.simcase['flutter_para']['method'] == 'k': - equations = KMethod(self, X0) - elif self.simcase['flutter_para']['method'] == 'ke': - equations = KEMethod(self, X0) - elif self.simcase['flutter_para']['method'] in ['pk', 'pk_schwochow']: - equations = PKMethodSchwochow(self, X0) - elif self.simcase['flutter_para']['method'] in ['pk_rodden']: - equations = PKMethodRodden(self, X0) - elif self.simcase['flutter_para']['method'] == 'statespace': - equations = StateSpaceAnalysis(self, X0) + if self.jcl.aero['method'] in ['freq_dom', 'mona_freq_dom']: + # Get initial solution from trim + X0 = self.response['X'][0, :] + # Select mona-based solution sequence + if self.simcase['flutter_para']['method'] == 'k': + equations = mona_frequency_domain.KMethod(self, X0) + elif self.simcase['flutter_para']['method'] == 'ke': + equations = mona_frequency_domain.KEMethod(self, X0) + elif self.simcase['flutter_para']['method'] in ['pk', 'pk_schwochow']: + equations = mona_frequency_domain.PKMethodSchwochow(self, X0) + elif self.simcase['flutter_para']['method'] in ['pk_rodden']: + equations = mona_frequency_domain.PKMethodRodden(self, X0) + else: + logging.error('Unknown mona-based flutter method: %s', self.simcase['flutter_para']['method']) + equations = None + elif (self.jcl.aero['method'] in ['mona_steady', 'mona_unsteady'] + and self.simcase['flutter_para']['method'] == 'statespace'): + # Get initial solution from trim + X0 = self.response['X'][0, :] + equations = mona_state_space.StateSpaceAnalysis(self, X0) + elif self.jcl.aero['method'] in ['cfd_freq_dom']: + # Select cfd-based solution sequence + if self.simcase['flutter_para']['method'] == 'k': + equations = cfd_frequency_domain.KMethod(self) + elif self.simcase['flutter_para']['method'] == 'ke': + equations = cfd_frequency_domain.KEMethod(self) + elif self.simcase['flutter_para']['method'] in ['pk_rodden']: + equations = cfd_frequency_domain.PKMethodRodden(self) + else: + logging.error('Unknown CFD-based flutter method: %s', self.simcase['flutter_para']['method']) + equations = None + else: + logging.error('Unknown aero method: %s', self.jcl.aero['method']) + equations = None response_flutter = equations.eval_equations() logging.info('Flutter analysis finished.') - for key in response_flutter.keys(): - self.response[key] = response_flutter[key] + for key, item in response_flutter.items(): + self.response[key] = item + self.successful = True + + def calc_gafs(self): + # Get initial solution from trim + X0 = self.response['X'][0, :] + Vtas = sum(X0[6:9] ** 2) ** 0.5 + # In case I decide to scale the GAFs with the dynamic pressure, I can use q_dyn from here + # q_dyn = self.response['q_dyn'][0] + + # Inline function to calculate reduced frequencies, Nastran definition + def f2k(f): + return 2.0 * np.pi * f * self.jcl.general['c_ref'] / 2.0 / Vtas + + # Get number of modes + n_modes_rbm = 5 + n_modes_flex = self.model['mass'][self.trimcase['mass']]['n_modes'][()] + PHIkh = self.model['mass'][self.trimcase['mass']]['PHIkh'][()] + n_modes = n_modes_rbm + n_modes_flex + # This is the index of each mode in the state vector X + idx_modes = list(range(1, n_modes_rbm + 1)) + list(range(1 + n_modes_rbm + 6, 1 + n_modes_rbm + 6 + n_modes_flex)) + logging.info('Calculating GAFs for %d rigid body modes and %d flexible modes...', + n_modes_rbm, n_modes_flex) + # Load matrices + PHIk_cfd = load_hdf5_sparse_matrix(self.model['PHIk_cfd']) + PHIcfd_cg = self.model['mass'][self.trimcase['mass']]['PHIcfd_cg'][()] + + # Step 1: set-up frequency parameters, generate pulse signal, and init storage + n_freqs = int(self.simcase['gaf_para']['fmax'] / self.simcase['gaf_para']['df']) + if n_freqs % 2 != 0: # n_freq is odd + n_freqs += 1 # make even + # Calculate all parameters from the number of freqs + fmax = n_freqs * self.simcase['gaf_para']['df'] + dt = 1.0 / fmax + t_final = 1.0 / self.simcase['gaf_para']['df'] + # Update simcase for time domain simulation + self.simcase['dt'] = dt + self.simcase['t_final'] = t_final + # Whole frequency space including negative frequencies + fftfreqs = fftfreq(n_freqs, dt) + # Positive only frequencies where we need to calculate the TFs and excitations + positiv_fftfreqs = np.abs(fftfreqs[:n_freqs // 2 + 1]) + # Only reduced frequencies < 3.0 are of interest + k = f2k(positiv_fftfreqs) + idx_k = np.where(k < 3.0)[0] + k_red = k[idx_k] + # Generate small-amplitude pulse signal + t, unit_pulse = polynomial_pulse(dt, t_final, eps=1.0) + # Scale the unit pulse for each mode such that the aplitudes are small. + # Right now the scaling is hard-codes based on test with the DC3, but might need to be adjusted + # for different configurations. On the other hand, I'm no fan of too many user-defined parameters... + pulse_factor = [1e-3, 1e-3, 1e-4, 1e-4, 1e-4] + [1e-3] * n_modes_flex + # The pulse's sign is used to align the aicraft rigid body motion with the nastran coordinate + # system (compatibility with DLM-based solutions). + pulse_sign = [1.0, -1.0, -1.0, 1.0, -1.0] + [1.0] * n_modes_flex + pulse_signal = [unit_pulse * factor for factor in pulse_factor] + pulse_signal = np.array(pulse_signal) + gust_f = fft(pulse_signal) + # Init storage for CFD forces + # To avoid an excessive amount of data, e.g. during unsteady cfd simulations, + # keep only the response data on the first mpi process (id = 0). + if self.myid == 0: + n_cfd = len(self.response['Pcfd'].squeeze()) + Pcfd_ref = np.zeros((n_cfd, len(t))) + Pcfd_pulse = np.zeros((n_cfd, len(t))) + Pcfd_gust = np.zeros((n_cfd, len(t))) + Pb_pulse = np.zeros((6, n_modes, len(t))) + Qhk = np.zeros((self.model['aerogrid']['n'][()] * 6, n_modes, len(k_red)), dtype=complex) + Qhh = np.zeros((n_modes, n_modes, len(k_red)), dtype=complex) + Qh_gust = np.zeros((n_modes, len(k_red)), dtype=complex) + + # Step 2: Run reference simulation without pulse + # Select CFD solution sequence and initialize + equations = CfdUnsteady(self, X0) + logging.info('Running reference time simulation for %g sec...', t_final) + # Loop over time steps + for i_step, t_step in enumerate(t): + X = copy.deepcopy(X0) + output_dict = equations.eval_equations(X, t_step, modus='sim_full_output') + if self.myid == 0: + Pcfd_ref[:, i_step] = output_dict['Pcfd'] + equations.finalize() + + # Step 3a: Run pulse simulations for all modes + for i_mode, idx_mode in zip(range(n_modes), idx_modes): + # Re-initialze CFD solution sequence for each mode + equations = CfdUnsteady(self, X0) + logging.info('Running small-amplitude pulse simulation for mode %d for %g sec...', i_mode + 2, t_final) + # Loop over time steps + for i_step, t_step in enumerate(t): + X = copy.deepcopy(X0) + X[idx_mode] += pulse_signal[i_mode, i_step] * pulse_sign[i_mode] + output_dict = equations.eval_equations(X, t_step, modus='sim_full_output') + if self.myid == 0: + Pcfd_pulse[:, i_step] = output_dict['Pcfd'] + equations.finalize() + + # Step 3b: Calculate TF for current mode + logging.info('Calculating transfer functions...') + if self.myid == 0: + # Compensate for initial condition and drift over time, transfer to aero grid 'k' + Pcfd = Pcfd_pulse - Pcfd_ref + Pk = PHIk_cfd.T.dot(Pcfd) + # Calculate transfer functions + Pk_f = fft(Pk, axis=1) + TF = Pk_f / (gust_f[i_mode, :]) + # Store + Qhk[:, i_mode, :] = TF[:, idx_k] + Pb_pulse[:, i_mode, :] = np.dot(PHIcfd_cg.T, Pcfd) + + # Step 4a: Run pulse simulation for gust mode + # Set-up small-amplitude 1-cosine gust with amplitude of 0.003 * Vtas + WG_TAS = 3e-3 + t, gust_signal, half_length = one_m_cosine_pulse(dt, t_final, Vtas, eps=WG_TAS * Vtas) + gust_f = fft(gust_signal) + self.simcase['gust'] = True + self.simcase['gust_orientation'] = 0 + self.simcase['gust_gradient'] = half_length + self.simcase['WG_TAS'] = WG_TAS + self.simcase['gust_para'] = {} + self.simcase['gust_para']['T1'] = 0.0 + # Select CFD solution sequence and initialize + equations = CfdUnsteady(self, X0) + logging.info('Running small-amplitude gust simulation for %g sec...', t_final) + # Loop over time steps + for i_step, t_step in enumerate(t): + X = copy.deepcopy(X0) + output_dict = equations.eval_equations(X, t_step, modus='sim_full_output') + if self.myid == 0: + Pcfd_gust[:, i_step] = output_dict['Pcfd'] + equations.finalize() + + # Step 4b: Calculate TF for gust mode + logging.info('Calculating transfer functions...') + if self.myid == 0: + # Compensate for initial condition and drift over time, transfer to aero grid 'k' + Pcfd = Pcfd_gust - Pcfd_ref + Pk = PHIk_cfd.T.dot(Pcfd) + # Calculate transfer functions + Pk_f = fft(Pk, axis=1) + TF = Pk_f / gust_f + # Store + Qk_gust = TF[:, idx_k] + Pb_gust = np.dot(PHIcfd_cg.T, Pcfd) + + if self.myid == 0: + # Because the CFD-based GAFs are calculated on the VLM/DLM aerogrid 'k', + # project also the initial trim solution on the k-set. + self.response['Pk_aero'] = PHIk_cfd.T.dot(self.response['Pcfd']) + # Apply modal transformation per frequency k_red to obtain Qhh + for i, _ in enumerate(k_red): + Qhh[:, :, i] = PHIkh.T.dot(Qhk[:, :, i]) + Qh_gust[:, i] = PHIkh.T.dot(Qk_gust[:, i]) + # Store results in response dictionary + self.response['desc'] = self.trimcase['desc'] + self.response['k_red'] = k_red + self.response['Qhk'] = Qhk + self.response['Qhh'] = Qhh + self.response['Qk_gust'] = Qk_gust + # The time signals are only saved for plotting / plausibility checking + self.response['pulse_signal'] = pulse_signal + self.response['gust_signal'] = gust_signal + self.response['t_pulse'] = t + self.response['Pb_pulse'] = Pb_pulse + self.response['Pb_gust'] = Pb_gust + self.successful = True diff --git a/loadskernel/solution_tools.py b/loadskernel/solution_tools.py index ce5ec924..8240f556 100644 --- a/loadskernel/solution_tools.py +++ b/loadskernel/solution_tools.py @@ -59,8 +59,7 @@ def design_gust_cs_25_341(simcase, atmo, V): MTOW = float(simcase['gust_para']['MTOW']) # Maximum Take-Off Weight MZFW = float(simcase['gust_para']['MZFW']) # Maximum Zero Fuel Weight fg = calc_fg(altitude, Z_mo, MLW, MTOW, MZFW) - logging.info( - 'CS25_Uds is set up with flight profile alleviation factor Fg = {}'.format(fg)) + logging.info('CS25_Uds is set up with flight profile alleviation factor Fg = %s', fg) # reference gust velocity (EAS) [m/s] if altitude <= 4572: @@ -134,3 +133,42 @@ def calc_fg(altitude, Z_mo, MLW, MTOW, MZFW): else: fg = fg_sl + (1.0 - fg_sl) * altitude / Z_mo return fg + + +def polynomial_pulse(dt, t_final, eps): + # Create a pulse signal with timestep dt up to t_final with magnitude eps. + # The pulse uses a 5th-order polynomial following eq. 2.26 in [1]. + # [1] Koch, C., “Whirl Flutter Stability Analysis Using Propeller Transfer Matrices”, + # Deutsches Zentrum für Luft- und Raumfahrt e. V. (DLR), 2024, https://doi.org/10.57676/BF00-1962. + t = np.arange(0.0, t_final + dt, dt) + # Pulse width in seconds, dt*40 should excite the low frequencies up to 5% of fmax. + tw = dt * 40 + if tw > t_final: + logging.warning('The pulse signal is longer than the simulation time. Please increase fmax and/or decrease df.') + # The pulse is assembled from two half pulses; the up and down strokes. + stroke_up = -4.0 * (2.0 * t / tw - 1.0)**5 - 15 * (2.0 * t / tw - 1.0)**4 - 20 * (2.0 * t / tw - 1.0)**3 \ + - 10 * (2.0 * t / tw - 1.0)**2 + 1 + stroke_down = +4.0 * (2.0 * t / tw - 1.0)**5 - 15 * (2.0 * t / tw - 1.0)**4 + 20 * (2.0 * t / tw - 1.0)**3 \ + - 10 * (2.0 * t / tw - 1.0)**2 + 1 + n_stroke = int(tw / 2 / dt) + n_lead = 10 + pulse = np.zeros(t.shape) + pulse[n_lead:n_lead + n_stroke] = stroke_up[:n_stroke] + pulse[n_lead + n_stroke:n_lead + n_stroke * 2] = stroke_down[n_stroke:n_stroke * 2] + pulse *= eps + return t, pulse + + +def one_m_cosine_pulse(dt, t_final, Vtas, eps): + # Create a pulse signal with timestep dt up to t_final with magnitude eps. + # The pulse uses the 1-cosine gust shape according to CS-25.341. + # The downside of the 1-cosine pulse is that it has zeros in the frequency domain. + # Because for gust analsysis mainly the low-frequency range is of interest, this + # pulse shape is acceptable and is implemented in most CFD codes. In addition, + # a pulse shorter than the shortest gust precribed in CS-25 (9-107m) should be suffcient. + t = np.arange(0.0, t_final + dt, dt) + half_length = 5.0 # meters + tw = half_length * 2.0 / Vtas + pulse = eps * 0.5 * (1 - np.cos(2.0 * np.pi * t / tw)) + pulse[np.where(t > tw)] = 0.0 + return t, pulse, half_length