Fluid-memory model (Retardation function): comparison between convolution calculation and linear models

One approach to develop time-domain model of marine structure is using Cummins equations:

(1)    \begin{equation*} (M_{RB}+A_{\infty})\ddot{\zeta}+\int_0^t K(t-\tau)\dot{\zeta}(\tau)d\tau+G\zeta=T_{exc}, \end{equation*}

where M_{RB} is the generalized mass matrix, A_{\infty} is the constant infinite-frequency added mass, G is restoring coefficient, K(t) is the impulse response function usually referred to the retardation or memory functions, T_{exc} is the wave exciting force and \zeta \in R^3 \times S^3 is the position and orientation of the marine structure.

Usually, frequency-dependent added mass A_{\omega} and potential damping B_{\omega} can be calculated by potential theory software, such as, by WAMIT. Then, K in frequency-domain and time-domain can be obtained from

(2)    \begin{align*} K(j\omega)=B(\omega)+j\omega(A(\omega)-A_{\infty}),\\ K(t)=\frac{2}{\pi}\int_0^{\infty}B(\omega)\cos(\omega t) d\omega, \end{align*}

To implement the simulation model, non-parametric fluid-memory model and linear time-invariant parametric model can be applied. The non-parametric fluid-memory model requires a discrete-time approximation of the convolution integral and enough past data. The convolution term becomes

(3)    \begin{equation*} \int_0^tK(t-\tau)\dot{\zeta}(\tau)d\tau\approx\sum_{i=0}^k K(k-i)\dot{\zeta}(i)\delta T  \end{equation*}

where \delta T is sampling time. This approach may be time-consuming and require a significant amount of computer memory.

Besides, the linear time-invariant model can be expressed by

(4)    \begin{align*} \dot{x} = \hat{A}x+\hat{B}\dot{\zeta}, \\ \int_0^t K(t-\tau)\dot{\zeta}(\tau)d\tau\approx \hat{C} x, \end{align*}

The parameters (\hat{A},\,\hat{B},\,\hat{C}) can be obtained by frequency-domain identification or time-domain identification. In frequency-domain, the identification problem is to find (\hat{A},\,\hat{B},\,\hat{C}) that

(5)    \begin{equation*} K(j\omega) \approx \hat{K}(j\omega)= \hat{C}(j\omega I-\hat{A})^{-1}\hat{B}. \end{equation*}

where \hat{K}(s) is matrix of transfer functions with entries

(6)    \begin{equation*} \hat{K}_{ik}=\frac{b_m s^m+b_{m-1}s^{m-1}+\cdots+b_0}{s^n+a_{n-1}s^{n-1}+\cdots+a_0}.  \end{equation*}

Note that \hat{K}_{ik} should be rational that n>m. Least-squares fitting method can be applied to solve the identification problem as follows:

(7)    \begin{align*} \theta^{\star}=\text{arg\ min}\sum_l \epsilon_l^*\epsilon_l, \\ \epsilon_l=K_{ik}(j\omega_l)-\hat{K}_{ik}(j\omega_l,\theta), \end{align*}

where \theta=[a_{n-1},\cdots,a_0,\,b_m,\cdots,b_0]^T. Details can be referred to the research paper [1]. Note that the method may not always obtain a stable estimation of (6). In this case, authors in [1] suggest to reflect the unstable poles about the imaginary axis and recompute the denominator polynomial.

On the other hand, time-domain identification methods are mostly based on the impulse response of the retardation function K_{i,k}. In the following, a method based on realization theory is introduced.

The Hankel matrix of the impulse response K_{ik} is

(8)    \begin{equation*} \math{H}_p=\left[ \begin{array}{rrrr} K_{ik}(1) & K_{ik}(2) & \cdots & K_{ik}(p) \\ K_{ik}(2) & K_{ik}(3) & \cdots & K_{ik}(p+1) \\ \vdots & \vdots & \vdots & \vdots \\ K_{ik}(p) & K_{ik}(p+1)& \cdots & K_{ik}(2p-1) \end{array} \right]. \end{equation*}

The state-space model (\mathbb{A},\  \mathbb{B},\ \mathbb{C},\ \mathbb{D}) can be obtained based on Singular Value Decomposition of \math{H}_p that

(9)    \begin{align*} \mathbb{A}=\Sigma_1^{-1/2}\left[ \begin{array}{ll} U_{11} &U_{22} \end{array} \right]\left[ \begin{array}{l} U_{12} \\ U_{13} \end{array} \right] \Sigma_1^{1/2}, \\ \mathbb{B}=\Sigma_1^{-1/2}V_{11}^*,\\ \mathbb{C}=U_{11}\Sigma_1^{1/2},\\ \mathbb{D}=\math{H}_0  \end{align*}

where \math{H}_0=K_{ik}(0) and \Sigma,\ U,\ V are calculated from Singular Value Decomposition

(10)    \begin{equation*} \math{H}_p=[U_1\ U_2] \left[ \begin{array}{rr} \Sigma_1 & 0 \\ 0 & \Sigma_2 \end{array} \right] [V_1^*\ V_2^*] \end{equation*}

and U_1=[U_{11}^T\ U_{12}^T\ U_{13}^T]^T and V_1=[V_{11}^*\ V_{12}^*\ V_{13}^*]^*.

The identification codes are given as follows:

function [sys]=imp2ss_zhz(t,K,Ord)
% effective only for SISO system
        dt = t(2)-t(1);
        %Hankel matrix and Singular Value Decomposition
        y = dt*K;
        h = hankel(y(2:end));
        [u,svh,v] = svd(h);
        u1 = u(1:length(K)-2,1:Ord);
        v1 = v(1:length(K)-2,1:Ord);
        u2 = u(2:length(K)-1,1:Ord);
        sqs = sqrt(svh(1:Ord));
        ubar = u1.'*u2;
        % discrete-time realization    
        a = ubar.*((1./sqs)*sqs.');
        b = v1(1,:).'.*sqs;
        c = u1(1,:).*sqs.';
        d = y(1);
        % continuous-time realization according to Tustin transform   
        iidd = inv(dt/2*(eye(Ord)+a));  
        ac = (a-eye(Ord))*iidd;         
        bc = dt*(iidd*b);            
        cc = c*iidd;                  
        dc = d-dt/2*((c*iidd)*b);    

In the following, the comparison between the time-domain model (9) and the frequency-domain model (6) using a spar-type floater, which is designed for NREL-5MW wind turbine [2].

The redartation function in heave direction and the results of two identification methods, and the errors of them are compared in the next figures.

From the results, it is observed that the model by time-domain identification method obtained a better result. The source data and MATLAB codes can be found here.

[1] Perez, T. and T. I. Fossen (2009). A Matlab Tool for Parametric Identification of Radiation-Force Models of Ships and Offshore Structures. Modelling, Identification and Control, MIC-30(1):1-15. doi:10.4173/mic.2009.1.1

[2] J. Jonkman (2010). Definition of the Floating System for Phase IV of OC3, Technical Report NREL/TP-500-47535.

Leave a Reply

Your email address will not be published. Required fields are marked *