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 default self.remove_predictor = True and the predictor term u_{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:

  1. one of the Gramian is integrated through the full Nyquist range

  2. 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: if get_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 partitions

    • order: 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

  1. 12 equally-spaced points integration of the Gramians in the low-frequency range [0,1.2] and

  2. 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 a pstats.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}}\]
nondimss()[source]

Scale state-space model based of self.ScalingFacts

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