Dynamic¶
- class sharpy.linear.src.linuvlm.Dynamic(tsdata, dt=None, dynamic_settings=None, integr_order=2, RemovePredictor=True, ScalingDict=None, UseSparse=True, for_vel=numpy.zeros)[source]¶
Class for dynamic linearised UVLM solution. Linearisation around steady-state are only supported. The class is built upon Static, and inherits all the methods contained there.
- Input:
tsdata: aero timestep data from SHARPy solution
dt: time-step
integr_order=2: integration order for UVLM unsteady aerodynamic force
RemovePredictor=True: if true, the state-space model is modified so as to accept in input perturbations, u, evaluated at time-step n rather than n+1.
ScalingDict=None: disctionary containing fundamental reference units:
{'length': reference_length, 'speed': reference_speed, 'density': reference density}
used to derive scaling quantities for the state-space model variables. The scaling factors are stored in
self.ScalingFact
.Note that while time, circulation, angular speeds) are scaled accordingly, FORCES ARE NOT. These scale by \(q_\infty b^2\), where \(b\) is the reference length and \(q_\infty\) is the dynamic pressure.
UseSparse=False: builds the A and B matrices in sparse form. C and D are dense anyway so the sparse format cannot be applied to them.
- - nondimss
normalises a dimensional state-space model based on the scaling factors in self.ScalingFact.
- - dimss
inverse of nondimss.
- - assemble_ss
builds state-space model. See function for more details.
- - assemble_ss_profiling
generate profiling report of the assembly and saves it into self.prof_out. To read the report:
import pstats p = pstats.Stats(self.prof_out)
- - solve_steady
solves for the steady state. Several methods available.
- - solve_step
solves one time-step
- - freqresp
ad-hoc method for fast frequency response (only implemented) for
remove_predictor=False
- Nx¶
Number of states
- Type
int
- Nu¶
Number of inputs
- Type
int
- Ny¶
Number of outputs
- Type
int
- K¶
Number of paneles \(K = MN\)
- Type
int
- K_star¶
Number of wake panels \(K^*=M^*N\)
- Type
int
- Kzeta¶
Number of panel vertices \(K_\zeta=(M+1)(N+1)\)
- Type
int
- Kzeta_star¶
Number of wake panel vertices \(K_{\zeta,w} = (M^*+1)(N+1)\)
- Type
int
To do: Upgrade to linearise around unsteady snapshot (adjoint)
- property Nu¶
Number of inputs \(m\) to the system.
- property Nx¶
Number of states \(n\) of the system.
- property Ny¶
Number of outputs \(p\) of the system.
- assemble_ss(wake_prop_settings=None)[source]¶
Produces state-space model of the form
\[\begin{split}\mathbf{x}_{n+1} &= \mathbf{A}\,\mathbf{x}_n + \mathbf{B} \mathbf{u}_{n+1} \\ \mathbf{y}_n &= \mathbf{C}\,\mathbf{x}_n + \mathbf{D} \mathbf{u}_n\end{split}\]where the state, inputs and outputs are:
\[\mathbf{x}_n = \{ \delta \mathbf{\Gamma}_n,\, \delta \mathbf{\Gamma_{w_n}},\, \Delta t\,\delta\mathbf{\Gamma}'_n,\, \delta\mathbf{\Gamma}_{n-1} \}\]\[\mathbf{u}_n = \{ \delta\mathbf{\zeta}_n,\, \delta\mathbf{\zeta}'_n,\, \delta\mathbf{u}_{ext,n} \}\]\[\mathbf{y} = \{\delta\mathbf{f}\}\]with \(\mathbf{\Gamma}\in\mathbb{R}^{MN}\) being the vector of vortex circulations, \(\mathbf{\zeta}\in\mathbb{R}^{3(M+1)(N+1)}\) the vector of vortex lattice coordinates and \(\mathbf{f}\in\mathbb{R}^{3(M+1)(N+1)}\) the vector of aerodynamic forces and moments. Note that \((\bullet)'\) denotes a derivative with respect to time.
Note that the input is atypically defined at time
n+1
, therefore by defaultself.remove_predictor = True
and the predictor termu_{n+1}
is eliminated through the change of state[1]:\[\begin{split}\mathbf{h}_n &= \mathbf{x}_n - \mathbf{B}\,\mathbf{u}_n \\\end{split}\]such that:
\[\begin{split}\mathbf{h}_{n+1} &= \mathbf{A}\,\mathbf{h}_n + \mathbf{A\,B}\,\mathbf{u}_n \\ \mathbf{y}_n &= \mathbf{C\,h}_n + (\mathbf{C\,B}+\mathbf{D})\,\mathbf{u}_n\end{split}\]which only modifies the equivalent \(\mathbf{B}\) and \(\mathbf{D}\) matrices.
References
[1] Franklin, GF and Powell, JD. Digital Control of Dynamic Systems, Addison-Wesley Publishing Company, 1980
To do: - remove all calls to scipy.linalg.block_diag
- assemble_ss_profiling()[source]¶
Generate profiling report for assembly and save it in self.prof_out.
- To read the report:
import pstats p=pstats.Stats(self.prof_out)
- balfreq(DictBalFreq, wake_prop_settings=None)[source]¶
Low-rank method for frequency limited balancing. The Observability ad controllability Gramians over the frequencies kv are solved in factorised form. Balancd modes are then obtained with a square-root method.
Details:
Observability and controllability Gramians are solved in factorised form through explicit integration. The number of integration points determines both the accuracy and the maximum size of the balanced model.
Stability over all (Nb) balanced states is achieved if:
one of the Gramian is integrated through the full Nyquist range
the integration points are enough.
Note, however, that even when stability is not achieved over the full balanced states, stability of the balanced truncated model with Ns<=Nb states is normally observed even when a low number of integration points is used. Two integration methods (trapezoidal rule on uniform grid and Gauss-Legendre quadrature) are provided.
Input:
DictBalFreq: dictionary specifying integration method with keys:
frequency
: defines limit frequencies for balancing. The balanced model will be accurate in the range [0,F], where F is the value of this key. Note that F units must be consistent with the units specified in the self.ScalingFacts dictionary.method_low
: [‘gauss’,’trapz’] specifies whether to use gauss quadrature or trapezoidal rule in the low-frequency range [0,F]options_low
: options to use for integration in the low-frequencies. These depend on the integration scheme (See below).method_high
: method to use for integration in the range [F,F_N], where F_N is the Nyquist frequency. See ‘method_low’.options_high
: options to use for integration in the high-frequencies.check_stability
: if True, the balanced model is truncated to eliminate unstable modes - if any is found. Note that very accurate balanced model can still be obtained, even if high order modes are unstable. Note that this option is overridden if “”get_frequency_response
: if True, the function also returns the frequency response evaluated at the low-frequency range integration points. If True, this option also allows to automatically tune the balanced model.
Future options:
truncation_tolerance
: ifget_frequency_response
is True, allows to truncate the balanced model so as to achieved a prescribed tolerance in the low-frequwncy range.Ncpu
: for parallel run
The following integration schemes are available:
trapz
: performs integration over equally spaced points using trapezoidal rule. It accepts options dictionaries with keys:points
: number of integration points to use (including domain boundary)
gauss
performs gauss-lobotto quadrature. The domain can be partitioned in Npart sub-domain in which the gauss-lobotto quadrature of order Ord can be applied. A total number of Npart*Ord points is required. It accepts options dictionaries of the form:partitions
: number of partitionsorder
: quadrature order.
Example:
The following dictionary
DictBalFreq={ 'frequency': 1.2, 'method_low': 'trapz', 'options_low': {'points': 12}, 'method_high': 'gauss', 'options_high': {'partitions': 2, 'order': 8}, 'check_stability': True }
balances the state-space model self.SS in the frequency range [0, 1.2] using
12 equally-spaced points integration of the Gramians in the low-frequency range [0,1.2] and
a 2 Gauss-Lobotto 8-th order quadratures of the controllability Gramian in the high-frequency range.
A total number of 28 integration points will be required, which will result into a balanced model with number of states
min{2*28* number_inputs, 2*28* number_outputs}
The model is finally truncated so as to retain only the first Ns stable modes.
- balfreq_profiling(wake_prop_settings=None)[source]¶
Generate profiling report for balfreq function and saves it into
self.prof_out.
The function also returns apstats.Stats
object.- To read the report:
import pstats p=pstats.Stats(self.prof_out).sort_stats('cumtime') p.print_stats(20)
- freqresp(kv, wake_prop_settings=None)[source]¶
Ad-hoc method for fast UVLM frequency response over the frequencies kv. The method, only requires inversion of a K x K matrix at each frequency as the equation for propagation of wake circulation are solved exactly. The algorithm implemented here can be used also upon projection of the state-space model.
Note: This method is very similar to the “minsize” solution option is the steady_solve.
- get_Cw_cpx(zval, settings=None)[source]¶
Produces a sparse matrix
\[\bar{\mathbf{C}}(z)\]where
\[z = e^{k \Delta t}\]such that the wake circulation frequency response at \(z\) is
\[\bar{\boldsymbol{\Gamma}}_w = \bar{\mathbf{C}}(z) \bar{\mathbf{\Gamma}}\]
- solve_steady(usta, method='direct')[source]¶
Steady state solution from state-space model.
Warning: these methods are less efficient than the solver in Static class, Static.solve, and should be used only for verification purposes. The “minsize” method, however, guarantees the inversion of a K x K matrix only, similarly to what is done in Static.solve.
- solve_step(x_n, u_n, u_n1=None, transform_state=False)[source]¶
Solve step.
If the predictor term has not been removed (
remove_predictor = False
) then the system is solved as:\[\begin{split}\mathbf{x}^{n+1} &= \mathbf{A\,x}^n + \mathbf{B\,u}^n \\ \mathbf{y}^{n+1} &= \mathbf{C\,x}^{n+1} + \mathbf{D\,u}^n\end{split}\]Else, if
remove_predictor = True
, the state is modified as\[\mathbf{h}^n = \mathbf{x}^n - \mathbf{B\,u}^n\]And the system solved by:
\[\begin{split}\mathbf{h}^{n+1} &= \mathbf{A\,h}^n + \mathbf{B_{mod}\,u}^{n} \\ \mathbf{y}^{n+1} &= \mathbf{C\,h}^{n+1} + \mathbf{D_{mod}\,u}^{n+1}\end{split}\]Finally, the original state is recovered using the reverse transformation:
\[\mathbf{x}^{n+1} = \mathbf{h}^{n+1} + \mathbf{B\,u}^{n+1}\]where the modifications to the \(\mathbf{B}_{mod}\) and \(\mathbf{D}_{mod}\) are detailed in
Dynamic.assemble_ss()
.Notes
Although the original equations include the term \(\mathbf{u}_{n+1}\), it is a reasonable approximation to take \(\mathbf{u}_{n+1}\approx\mathbf{u}_n\) given a sufficiently small time step, hence if the input at time
n+1
is not parsed, it is estimated from \(u^n\).- Parameters
x_n (np.array) – State vector at the current time step \(\mathbf{x}^n\)
u_n (np.array) – Input vector at time step \(\mathbf{u}^n\)
u_n1 (np.array) – Input vector at time step \(\mathbf{u}^{n+1}\)
transform_state (bool) – When the predictor term is removed, if true it will transform the state vector. If false it will be assumed that the state vector that is parsed is already transformed i.e. it is \(\mathbf{h}\).
- Returns
Updated state and output vector packed in a tuple \((\mathbf{x}^{n+1},\,\mathbf{y}^{n+1})\)
- Return type
Tuple
Notes
- To speed-up the solution and use minimal memory:
solve for bound vorticity (and)
propagate the wake
compute the output separately.
- unpack_state(xvec)[source]¶
Unpacks the state vector into physical constituents for full order models.
The state vector \(\mathbf{x}\) of the form
\[\mathbf{x}_n = \{ \delta \mathbf{\Gamma}_n,\, \delta \mathbf{\Gamma_{w_n}},\, \Delta t\,\delta\mathbf{\Gamma}'_n,\, \delta\mathbf{\Gamma}_{n-1} \}\]Is unpacked into:
\[{\delta \mathbf{\Gamma}_n,\, \delta \mathbf{\Gamma_{w_n}},\, \,\delta\mathbf{\Gamma}'_n}\]- Parameters
xvec (np.ndarray) – State vector
- Returns
Column vectors for bound circulation, wake circulation and circulation derivative packed in a tuple.
- Return type
tuple