Sunday, 20 April 2014

Test post for blog disemination of journal paper



Phasor modelling is not frequency domain modelling. 

There, I said it. I will now describe why this is the case in four steps:
  1. Background on linear WEC modelling: why phasor models are used to provide coefficients for frequency domain models; how radiation memory is encoded in the frequency domain.
  2. Why frequency domain and phasor modelling appear similar: mainly, more subjectivity than desirable required to interpret the maths.
  3. Differences between frequency domain and phasor models: like the difference between cave-aged cheddar and vegan cheese; or Belgium beer vs non-alcoholic lager.
  4. Why this distinction is important: for making sense of the literature, not breaking the laws of physics, and ensuring we use the right tools for research and development of wave energy converters (WECs).




1. Background on linear modelling of WECs


The forces on a WEC due to water waves (wave-body interactions) depend, among other things, on the shape of the surface in contact with the water (its wetted surface). It is certainly possible to run a simulation where the wetted surface is calculated at each time step; this is done for example when modelling parametric roll of ships [1]. However, it is only recently that we have the computational power to run long simulations in this manner. So the most common approach used for hydrodynamic modelling of wave-body interactions is to start with the assumption that the changes in wetted surface during any simulation are negligible. All wave-body forces are calculated on the equilibrium wetted surface.

The next step is to assume that the forces due to wave-body interactions are linear. This assumption allows us to decompose the problem into managable chunks. Superposition allows the wave-body forces to be expressed as the sum of forces due to the incident wave (the Froude-Krylov force), the diffracted wave, and the wave that is caused by the body's motion (the radiated wave). Superposition also allows us to decompose the wave into manageable components. In this case, it is convenient to consider sinusoidal components, because the physics theory of water waves (potential flow theory) has a strong convention of considering sinusoidal waves.


Frequency domain modelling

Another reason why considering sinusoidal components is convenient is that it allows us to use Fourier transforms, a tool that is widely used in many branches of engineering. Most signals describing dynamic systems are time domain: that is, they consist of a parameter that varies with respect to time. The domain is the independent variable for which the signal is defined.

Fourier Transformation allows equivalent representations of signals. The inverse Fourier Transform:

\begin{equation}
f(t) = \int_{-\infty}^{\infty} F(\omega)\mathrm e^{i\omega t}d \omega
\label{eq:IFT}
\end{equation}
tells us that any function can be expressed as a weighted sum (integral to be precise) of complex exponentials. These weighting factors are themselves complex, as they specify the amplitude and phase of each constituent complex exponential. The set of weighting factors is known as the Fourier Transform, which is itself an integral involving the original signal, \(F(\omega)= \mathscr{F}[f(t)] =1/{2 \pi}\int_{-\infty}^{\infty} f(t)\mathrm e^{-i\omega t}dt \), and contains all the information required to reconstruct the original signal.

If you take the Fourier transform of both sides of the equation for a linear time-invariant system: \(\mathscr{F}[y(t)]= \mathscr{F}[h(t) \ast x(t)]\), then this is an equivalent frequency domain model of that system: \(Y(\omega) = H(\omega) X(\omega)\). The frequency domain model can be used to find the steady state response, because this depends on the frequency of excitation. The frequency domain model cannot be used to find the transient response of the (stable) system, because the transient is a decaying signal, and depends on the conditions at one point in time \(t=0\), neither of  which are suited to Fourier Transform representation.

Note that at this point in the modelling process we are considering the waves and response to be composed of many sinusoids; we are simply dealing with them individually. However, in order to to calculate the hydrodynamic forces at each frequency, it is necessary to use potential flow theory. Potential flow theory has a long history of limiting the problem to single frequency waves; Cummins [2] noted that research prior to 1950 was "primarily concerned with sinusoidal responses to sinusoidal waves".  The sinusoidal assumption is certainly useful, as the expressions for fluid potentials and pressures are not initially linear with respect to wave elevation or motion; linearisation is a later step that allows easy numerical solution. Note this starting assumption is a convention rather than a necessity, and Cummins [2] relates the radiation force to the fluid potentials without first assuming sinusoidal motion. Nevertheless, this is a very strong convention, and as we are borrowing from a branch of physics with such a long history, it is convenient to respect the associated assumptions and notation.


Phasor notation

The mathematical notation used in potential flow theory is known as phasor notation. Phasor notation can be used when the problem has been limited to sinusoidal oscillation. Any sinusoidal signal \(f(t)\) of amplitude \(f_0\), frequency \(\omega\), and phase \(\phi\):

\begin{equation}
f(t) = f_0 \cos (\omega t + \phi) = f_0 \Re[\mathrm e^{i\omega t + i \phi}] = \Re[\widehat{F} \mathrm e^{i\omega t }]
\label{eq:phasor}
\end{equation}
can be expressed in terms of a complex amplitude \(\widehat{F} = f_0 \mathrm e^{i \phi}\)  multiplied by a complex exponential. If all the signals in a linear system are expressed in this manner, the complex exponential \(\mathrm e^{i \omega t}\)  term can be cancelled from both sides of the equation. This results in an equation that relates the complex amplitudes only. The sinusoidal nature of the system is ignored; only the vital characteristics of the signals are considered: their amplitudes and phases. These complex amplitudes are known as phasors.


Using hydrodynamic coefficients in linear models

By assuming small amplitude waves and motions, we can linearise the fluid forces given by potential flow theory. This allows the excitation force (Froude-Krylov plus diffracted) to be expressed as a linear function of wave elevation. Likewise the radiation force can be expressed as a linear function of motion (e.g. velocity).

Now that these expressions are linear, we can plug them into a frequency domain model of a WEC \eqref{eq:IFT}. You need to set up and solve many phasor models (of potential flow theory) to provide the hydrodynamic coefficients for one frequency domain model. For example, if you use WAMIT to calculate the hydrodynamic coefficients for a given geometry, the formulation and solution of the equations is repeated for each frequency specified. On the other hand, when the frequency domain problem using these coefficients is solved, this is one equation, which is solved in one step.

There are a couple ways that coefficients calculated for the single frequency problem can be used in frequency domain models. One option is to use a frequency domain model based on the discrete Fourier transform (DFT). Individual phasor models are solved, one for each frequency bin of the DFT. Later the inverse DFT could be used to transform the results back into the time domain. Another option is to solve the potential flow problem at various frequencies of interest, e.g. those where the hydrodynamic coefficients vary most with respect to frequency. Curve-fitting methods (e.g. state-space) are used to approximate these coefficients with a continuous function. The inverse Fourier transform could be used to return to the time domain.


Hydrodynamic added mass and damping

The frequency domain coefficients that relate wave elevation to excitation force, and velocity to radiation force, are complex; i.e. there is a transformation in amplitude and phase at each frequency. There is a long-standing tradition in hydrodynamics to represent the radiation coefficient in a particular manner. The complex radiation coefficient is split into its real and imaginary components. The frequency is factored out from the imaginary term: \(F_r(\omega) = [B(\omega)+ i\omega M(\omega)]U(\omega)\); this combined with velocity now has a similar form to mass-acceleration: \(F_r(\omega) = M(\omega) i\omega U(\omega) + B(\omega)\). It differs of course because \(M(\omega)\)  varies with frequency. However, if we are considering the sinusoidal special case, and using phasor notation, then radiation really does look like the action of mass and damping. Hence the term \(M(\omega)\)  is known as added mass, and \(B(\omega)\)  as radiation (or added) damping.


How added mass and damping represent radiation memory

Recalling that frequency domain system representation is found by taking the Fourier transform of both sides of an equation of (time domain) motion for a linear system, it is clear that \(M(\omega)\) and \(B(\omega)\) must be related to the time domain representation of radiation force via the inverse Fourier transform. The time domain radiation force includes a 'memory term'; a convolution whose intimidating appearance belies its simple implication: the waves created by recent (as well as present) motion place a force on the body. As the radiation force at any instant depends on past motion, we describe this as 'memory'. The most common formulation is widely described as the Cummins equation [2], in which the convolution is between velocity and a 'kernel'. The kernel is a time domain function representing a causal impulse response. Causality is indicated by this kernel having a value of zero for \(t<0\).  

The causal nature of the radiation force is encoded in the values of added mass and damping. Using a clever trick involving odd and even functions, it can be shown that added mass and radiation damping are related to each other via the Kramers-Kronig relations  (See 'related proof from the time domain' on the Wikipedia page). This relationship between added mass and damping confirms (and incodes the information) that the radiation behaviour is causal. The frequency dependency of these coefficients indicates that the system has memory (if they were constants the radiation force would be causal but only related to present motion).




2. Reasons why frequency domain and phasor modelling appear similar


Overlap in practical implementation:

The concept of decomposing a model of a WEC (its equation of motion) into sinusoidal components, performing calculations on parameters for each frequency component and then using superposition to represent the full system once more is the basic premise behind frequency domain modelling. However, it is also possible to use phasor modelling iteratively for this purpose. In fact phasor modelling is required to express the physics theory used to calculate the radiation and excitation coefficients. This is one of the reasons why phasor and Fourier modelling appear so similar.

Same mathematical notation used for equality and assignment:

A major reason why phasor notation and frequency domain modelling are difficult to distinguish between is that we do not use different mathematical notation to distinguish between equality and assignment. This distinction is recognised in computer programming languages. Some branches of mathematics use the \(:=\) symbol to assign values to a variable. However, in engineering the \(=\) sign is used to indicate both assignment and equality. So in equation \eqref{eq:IFT} the \(=\) is an equality. The equality shows that any (integrable) function, \(f(t)\), has an equivalent representation. On the other hand, in equation \eqref{eq:phasor}, the first \(=\) sign is an assignment, and subsequent usage indicates equality. The assignment restricts \(f(t)\) to a special case. The equalities in this case show that the restricted case has equivalent representations.

Same mathematical notation used for a constant and a domain:

Another major reason why phasor notation and frequency domain modelling look so similar is that the notation \(\omega\) is used to indicate a constant in phasor notation, and the domain in a Fourier transform. In equation \eqref{eq:IFT} the integration over \(\omega\)  shows us that this is not a constant. The familiarity of the Fourier transform in engineering applications leads us to interpret the frequency domain representation of a signal \(\mathscr{F}[f(t)]= F(\omega)\) , as a function with respect to the independent variable \(\omega\)  (its domain). In the assignment in equation \eqref{eq:phasor}, \(\omega\) can be seen to be a constant. The function \(f(t)\) has been restricted to the special case of single frequency oscillation; that frequency being \(\omega\).

Similar mathematical form used for representation of systems:

A linear model of a WEC, neglecting losses, can be expressed in the frequency domain as:

\begin{equation}
F_e(\omega) = \left(i \omega[m + M(\omega)] + B(\omega) -ic / \omega +z_{pto}\right) U(\omega)
\label{eq:WEC fd}
\end{equation}
with \(F_e(\omega)\)  the excitation force, \(m\) the physical mass, \(c\) the restoration coefficient, and \(z_{pto}\) the power take off impedance. As mentioned earlier, the added mass and damping are frequency dependent. The corresponding model for the special case of single frequency oscillation at \(\omega\) could be expressed as follows:

\begin{equation}
\Re[\widehat{F_e}\mathrm e^{i\omega t }]= \Re[\left(i \omega[m + M] + B -ic / \omega +z_{pto}\right) \widehat{U}\mathrm e^{i\omega t }]
\label{eq:WEC sin}
\end{equation}
with \(\widehat{F_e}\) the complex amplitude of wave excitation at \(\omega\) and \(\widehat{U}\) the complex amplitude of velocity at \(\omega\) . The added mass and damping are constants, as we are only considering one frequency. The sinusoidal oscillation is represented by \(\Re[\mathrm e^{i\omega t }]\). As this appears on both sides of the equation, they can be cancelled. The relationship between excitation and response may then be expressed in terms of the complex exponentials alone:

\begin{equation}
\widehat{F_e}= \left(i \omega[m + M] + B -ic / \omega +z_{pto}\right) \widehat{U}
\label{eq:WEC phasors}
\end{equation}
While equation \eqref{eq:WEC sin} explicitly states the sinusoidal assumption, the more commonly used phasor model \eqref{eq:WEC phasors} does not. Its form has similarities to the corresponding frequency domain model \eqref{eq:WEC fd}.


Same mathematical notation used for a constant and a function of frequency:

The phasor model \eqref{eq:WEC phasors} can be expressed in a way which makes it look even closer to the frequency domain model \eqref{eq:WEC fd}. It is fairly common practice to express the phasor model as:

\begin{equation}
\widehat{F_e}= \left(i \omega[m + M(\omega)] + B(\omega) -ic / \omega +z_{pto}\right) \widehat{U}
\label{eq:WEC phasor ambiguous notation}
\end{equation}
This notation serves to indicate that the values of added mass and damping depend on the chosen frequency, whereas physical mass and spring do not. Equations \eqref{eq:WEC fd} and \eqref{eq:WEC phasor ambiguous notation} now look very similar indeed. This is because the same notation is used to indicate both constants and functions of frequency. In the frequency domain equation \eqref{eq:WEC fd}, \(M(\omega)\) is a parameter that varies over the independent variable of \(\omega\). In equation \eqref{eq:WEC phasor ambiguous notation} \(M(\omega)\) is used to indicate a constant.

I do not like the notation used in \eqref{eq:WEC phasor ambiguous notation} at all. Pedantry compels me to note its inconsistency; after all, the excitation force and velocity also vary with frequency. However, what I find really problematic is how frighteningly easy it is to read the \(\omega\) as a domain, and interpret \eqref{eq:WEC phasor ambiguous notation} as a 'frequency domain equation'.


3. Differences between frequency domain and phasor modelling


The frequency domain representation of a signal is a transform. The same signal is transformed from one domain to another: \(f(t) \stackrel{\mathscr{F}}{\leftrightarrow} F(\omega)\) . The domain of this equivalent representation is now frequency.

For phasor representation of a signal, we start off with an assignment: \(f(t) = f_0 \cos (\omega t + \phi)\) . The signal \(f(t)\) is now limited to a special case: oscillation at one frequency, \(\omega\) . The domain remains time. In \(f_0 \cos (\omega t + \phi)\)  the signal varies with \(t\), while \(\omega\)  is a constant. When we express it in the following manner: \(f(t) = \Re[\widehat{F} \mathrm e^{i\omega t }]\), then again \(t\) is the domain, and \(\omega\) is a constant.

Equation \eqref{eq:WEC fd} is a frequency domain model of a WEC. The excitation \(F_e(\omega)\)  and response \(U(\omega)\)  are signals whose domains are frequency. They are Fourier transforms of the corresponding time domain signals. The added mass \(M(\omega)\)  and damping \(B(\omega)\)  are also parameters that have domains of frequency. As they are related via the Kramers-Konig relations, and vary with frequency, they convey the information that the time domain radiation behaviour is causal and has 'memory' of past motion.

Equation \eqref{eq:WEC sin} is a time domain model of a WEC limited to a special case. The excitation and response are signals whose domains are time. This is clearer from the sinusoidal form; e.g. wave excitation is a sinusoid with amplitude \(f_{e0}\) and phase \(\phi_f\): \(\Re[\widehat{F_e}\mathrm e^{i\omega t }]= f_{e0} \cos(\omega t + \phi_f)\). The model has been limited to oscillation at a single frequency, and so \(\omega\) is a constant. Added mass and damping are constants related via the Kramer-Kronig relations. Consequently, there is not enough information to indicate whether radiation has memory of past motion, but we do know it does not depend on future motion (it is causal).

Equation \eqref{eq:WEC phasors} is a phasor model of a WEC. The excitation and response are not signals that have domains of time or frequency. They are constants; condensed descriptions of the (limited class of) time domain signals; the time varying parts of which have been canceled from both sides of equation \eqref{eq:WEC sin}. As in \eqref{eq:WEC sin}, the frequency, added mass, and radiation damping are constants. Radiation memory is not represented.


4. Why this distinction is important


Sinusoidal assumption often obscure in literature:

One implication (or cause) of the ambiguity of the \(=\) sign is that assumptions stated mathematically in hydrodynamics articles are often not repeated in writing. For example, when we describe a velocity potential as: \(\varphi (x,y,z,t) = \Re[\Phi(x,y,z)\mathrm e^{-i\omega t}]\), researchers familiar with potential flow theory will understand this as an assignment rather than an equality. As the assumption of a sinusoidal wave is explicitly stated by the equation, often authors will not spell this out in words. If a novice researcher is looking to the text to explain the equations, they will miss this criticial modelling assumption.

Phasor model treated as a frequency domain model:

If the complex amplitudes in equation \eqref{eq:WEC phasor ambiguous notation} are interpreted as Fourier transforms while the added mass and damping are treated as constants that depend on the frequency of a chosen sinusoid, then the corresponding time domain equation of motion would be:

\begin{equation}
\mathscr{F}[F_e(\omega)] = f_e(t) = [m+M(\omega)]\ddot{x} + [B(\omega)+ z_{pto}]\dot{x} + cx
\label{eq:WEC msd}
\end{equation}
When I was new to wave energy, this was my initial interpretation. Diagrams showing systems of masses, springs and dampers can be useful for describing force flow and kinematics, but these also guide thinking towards a mass-spring-damping model.

In Cummins's excellent paper [2] he describes how 'classical' ship-keeping theory had assumed sinusoidal waves and responses using equations of the form \eqref{eq:WEC phasor ambiguous notation}.  In the 1950s, when responses to spectrums were increasingly of interest, a common approach was to use a second-order differential equation taking the form of \eqref{eq:WEC msd}. Cummins's paper was in reponse to this practice, which he discredited from a theoretical perspective, as well as cricticising its replication of experimental results.

Cummins's approach is now widely used in ship-keeping and wave energy, but these ideas took some time to catch on, and there remains literature that uses a mass-spring-damper model with frequency dependent radiation parameters. In some cases, this may a conscious modelling choice, either as a tool for explaining general principles, or as a trade-off between simplicity and accuracy. For instance, the accuracy forfeited by using a mass-spring-damper model could be much less than that forfeited by using a linear model when the behaviour of interest is non-linear. One quick way to build and run a non-linear time domain model is to base it around a mass-spring-damper representation of the linear dynamics. Of course caution is required in such an approach; Greenhow and White [4] demonstrated that in sinusoidal waves, the mass-spring-damper model gave the same results as the model that included radiation memeory when linear PTO was used, but significantly different results when non-linear PTO was used.


PTO optimised at each frequency requires future knowledge of velocity:

However, when the radiation memory is a vital part of the model, then the frequency dependency of the added mass and damping needs to be modelled. In particular, for the optimal control problem, there is a crucial distinction between the single frequency and the general multiple frequency case. In single frequency oscillation, it is possible to choose the PTO impedance \(Z_{pto}\) according to optimal control theory, without 'future knowledge' of the exciting wave. Of course, you know it is a sinusoid of a given frequency, so it could be argued you do indeed have sufficient future knowledge! In the general case of multiple frequencies, the optimal PTO force requires knowledge of the future motion [3]. You would need to know the wave excitation on the system for as far into the future as the radiation force 'remembers' the past motion. This is because the response now is 'remembered' that far into the future; you are not only optimising the present response, but also the future response.

Results from single frequency models require care when applied to multiple frequency cases. For example, you could measure the proportion of incident wave power captured in experimental tests in regular waves, and then use superposition to calculate the linear response to a spectrum. However, if the PTO coefficients were optimised for each frequency, the power capture would be vastly overestimated. In addition to the usual problem that spectra of interest involve non-linear everything, there is the problem that you would need to know future velocity in order to make use of those PTO coefficients.


Type of model often obscure in literature:

Equation \eqref{eq:WEC phasor ambiguous notation} is a representation of the phasor model that looks similar to a frequency domain model. In fact, the term 'frequency domain' modelling is often used to describe phasor models, particularly when they are in this form. Occasionally the description of 'frequency domain' is extended to time domain models limited to single frequency oscillation.This can make some literature ambiguous.


Single frequency models referred to as 'frequency domain':

An example of this type of ambiguity is found in a paper by Greenhow and White [4]. The authors had noted that it was not uncommon for the phasor model to be interpreted and implemented as a mass-spring-damper model \eqref{eq:WEC msd}. They sought to warn readers against using this model for anything other than sinusoidal oscillation. The paper referred to this model, as was common at that time, as the 'frequency domain' model. The equations do not specify domain, but the process of stepping through the equations, with the aid of the guiding text, shows they are indeed referring to a linear time domain equation where the excitation force has been restricted to a single frequency oscillation (of course, I checked with Prof Greenhow too, just to make sure!).


Frequency domain model treated as a single frequency model:

The impact of attaching the label 'frequency domain' to models defined for a limited subset of behavior (single frequency) is that there is a wide-spread mistrust of frequency domain models for use in the general (multiple frequency) case. This mistrust ranges from uncertainty over whether a (slower) time domain simulation is required if the linear model is excited by a spectrum, to the belief that a frequency domain model 'does not represent memory', and so will not give the same steady-state solution as the corresponding linear time domain model.

By casting aspersions on a useful tool that every other branch of engineering values, there is a risk that research on wave energy could be hindered. For instance, when developing a non-linear time domain model using the Cummins' equation as a core, the frequency domain expression of the linear core is an extremely useful tool for verifying model configuration, and for validating the linear time domain model. 


Acknowledgements

Mathjax allowed me to use Latex equations in Blogger. The instructions on Sam Buss's blog worked for me: http://mathjaxtest.blogspot.co.uk/. I also enabled equation numbering and referencing as described in: http://docs.mathjax.org/en/latest/tex.html. Writing this content from my home office has been greatly helped by using Olga Degtyareva's productivity tools: http://www.olgadegtyareva.com/


References

[1] Ribeiro e Silva and Guedes Soares (2013), 'Prediction of parametric rolling in waves with a time domain non-linear strip theory model' , Ocean Engineering 72  453–469


[2] Cummins (1962), 'The impulse response function and ship motions.' Schiffstechnik 9:101–109

[3] Naito and Nakamura (1986), 'Wave energy absorption in irregular waves by feed-forward control system'. In Hydrodynamics of ocean wave-energy utilization, ed. by D.V. Evans & A.F. de O. Falcão, 269-280. IUTAM Symposium, Lisbon.

[4] Greenhow and White (1997), 'Optimal heave motion of some axisymmetric wave energy devices in sinusoidal waves' Applied Ocean Research Volume 19, Issues 3–4, June–August 1997, Pages 141–159


Image credits

Kramers Kronig relations from Wikipedia: http://en.wikipedia.org/wiki/File:KramersKronig.svg
Ying Yang duckies: bloggers own photo (Inverness duck race)

No comments:

Post a Comment