In this tutorial we will compare the results of Workshop 3 and Workshop 4. They calculate the solution of a simplified jacket wind turbine. We assume that the structure can be modelled as a space frame with structural elements that are subject to axial displacement and bending. Therefore, the equations of motion at a given element are given by:
\[ \rho A \frac{\partial^2 u (x,t)}{\partial t^2} - E A \frac{\partial^2 u(x,t)}{\partial x^2} = q_u(x) \]
\[ \rho A \frac{\partial^2 v (x,t)}{\partial t^2} + E I \frac{\partial^4 v(x,t)}{\partial x^4} = q_v(x) \]
The structure is schematically shown here:
# For comparison purposesimport time
Part 0: Overlapping commands#
Let’s start by setting some parameters:
import numpy as npimport matplotlib.pyplot as pltimport scipy.linalg as scpimport scipy.integrate as scpi# Input parameters# NacelleMass_Nacelle = 5000 # [kg]# Pilem_Pile = 1000 # [kg/m]EI_Pile = 1e9 # [N.m2]EA_Pile = 1e8 # [N]H_Pile = 100 # [m]# Jacketm_Jacket = 100 # [kg/m]EI_Jacket = 1e7 # [N.m2]EA_Jacket = 1e6 # [N]H_Jacket = 70 # [m]Base_Jacket = 30.23 # [m]Top_Jacket = 16.13 # [m]# Define load parametersf0 = 2 # [Hz]A0 = 0.1 # [m]T0 = 20 # [s]# Define output time vectordt = 0.01 # [s]T = np.arange(0, 5*T0, dt)nT = len(T)
Step 1: discretize the domain#
As usual, we start by discretizing the domain. We will discretize the structure using 21 nodes.
Here we first define a matrix with the nodal coordinates.
# Define node coordinatesNodeCoord = ([-(Base_Jacket*4+Top_Jacket*0)/4/2, -4*H_Jacket/4], # Point 1 [(Base_Jacket*4+Top_Jacket*0)/4/2, -4*H_Jacket/4], # Point 2 [-(Base_Jacket*3+Top_Jacket*1)/4/2, -3*H_Jacket/4], # Point 3 [(Base_Jacket*3+Top_Jacket*1)/4/2, -3*H_Jacket/4], # Point 4 [-(Base_Jacket*2+Top_Jacket*2)/4/2, -2*H_Jacket/4], # Point 5 [(Base_Jacket*2+Top_Jacket*2)/4/2, -2*H_Jacket/4], # Point 6 [-(Base_Jacket*1+Top_Jacket*3)/4/2, -1*H_Jacket/4], # Point 7 [(Base_Jacket*1+Top_Jacket*3)/4/2, -1*H_Jacket/4], # Point 8 [-(Base_Jacket*0+Top_Jacket*4)/4/2, 0*H_Jacket/4], # Point 9 [0, 0*H_Jacket/4], # Point 10 [(Base_Jacket*0+Top_Jacket*4)/4/2, 0*H_Jacket/4], # Point 11 [0, 10], # Point 12 [0, 20], # Point 13 [0, 30], # Point 14 [0, 40], # Point 15 [0, 50], # Point 16 [0, 60], # Point 17 [0, 70], # Point 18 [0, 80], # Point 19 [0, 90], # Point 20 [0, 100]) # Point 21nNode = len(NodeCoord)
Once we have the node coordinates we proceed to define the elemental connectivities. Here, we will use the same array to assign the material properties to each element. Note that they are different depending on which part of the structure they belong to.
# Define elements and their properties# NodeLeft NodeRight m EA EIElements = ( [1, 3, m_Jacket, EA_Jacket, EI_Jacket], [1, 4, m_Jacket, EA_Jacket, EI_Jacket], [2, 4, m_Jacket, EA_Jacket, EI_Jacket], [3, 4, m_Jacket, EA_Jacket, EI_Jacket], [3, 5, m_Jacket, EA_Jacket, EI_Jacket], [5, 4, m_Jacket, EA_Jacket, EI_Jacket], [4, 6, m_Jacket, EA_Jacket, EI_Jacket], [5, 6, m_Jacket, EA_Jacket, EI_Jacket], [5, 7, m_Jacket, EA_Jacket, EI_Jacket], [5, 8, m_Jacket, EA_Jacket, EI_Jacket], [6, 8, m_Jacket, EA_Jacket, EI_Jacket], [7, 8, m_Jacket, EA_Jacket, EI_Jacket], [7, 9, m_Jacket, EA_Jacket, EI_Jacket], [9, 8, m_Jacket, EA_Jacket, EI_Jacket], [8, 11, m_Jacket, EA_Jacket, EI_Jacket], [9, 10, m_Jacket, EA_Jacket, EI_Jacket], [10, 11, m_Jacket, EA_Jacket, EI_Jacket], [10, 12, m_Pile, EA_Pile, EI_Pile], [12, 13, m_Pile, EA_Pile, EI_Pile], [13, 14, m_Pile, EA_Pile, EI_Pile], [14, 15, m_Pile, EA_Pile, EI_Pile], [15, 16, m_Pile, EA_Pile, EI_Pile], [16, 17, m_Pile, EA_Pile, EI_Pile], [17, 18, m_Pile, EA_Pile, EI_Pile], [18, 19, m_Pile, EA_Pile, EI_Pile], [19, 20, m_Pile, EA_Pile, EI_Pile], [20, 21, m_Pile, EA_Pile, EI_Pile])nElem = len(Elements)
Let’s plot the structure to make sure that it looks like the model in the figure.
plt.figure()for iElem in np.arange(0, nElem): NodeLeft = Elements[iElem][0]-1 NodeRight = Elements[iElem][1]-1 m = Elements[iElem][2] EA = Elements[iElem][3] EI = Elements[iElem][4] if m == m_Jacket and EA == EA_Jacket and EI == EI_Jacket: c = 'b' elif m == m_Pile and EA == EA_Pile and EI == EI_Pile: c = 'r' else: print("ERROR: unknown material. Check your inputs.") break plt.plot([NodeCoord[NodeLeft][0], NodeCoord[NodeRight][0]], [NodeCoord[NodeLeft][1], NodeCoord[NodeRight][1]], c=c)plt.axis('equal')plt.grid()
Step 2: define the shape functions#
Here we will use the exact same functions as in tutorial 4 and 5. Linear shape functions for the axial displacement and cubic shape functions for the deflection and rotations. Since we already know its expression and we already have the value of the elemental matrices, we skip this step in this tutorial.
Step 3: computation of the elemental matrices#
In the theory we have seen that the mass and stiffness elemental matrices for the space frame using linear and cubic shape functions are given by:
\[\begin{split} \boldsymbol{M} = \frac{mL}{420} \begin{bmatrix} 140 & 0 & 0 & 70 & 0 & 0 \\ 0 & 156 & 22L & 0 & 54 & -13L \\ 0 & 22L & 4L^2 & 0 & 13L & -3L^2 \\ 70 & 0 & 0 & 140 & 0 & 0 \\ 0 & 54 & 13L & 0 & 156 & -22L \\ 0 & -13L & -3L^2 & 0 & -22L & 4L^2 \end{bmatrix} \end{split}\]
\[\begin{split} \quad \boldsymbol{K} = \begin{bmatrix} \frac{EA}{L} & 0 & 0 & \frac{-EA}{L} & 0 & 0 \\ 0 & \frac{12EI}{L^3} & \frac{6EI}{L^2} & 0 & \frac{-12EI}{L^3} & \frac{6EI}{L^2} \\ 0 & \frac{6EI}{L^2} & \frac{4EI}{L} & 0 & \frac{-6EI}{L^2} & \frac{2EI}{L} \\ \frac{-EA}{L} & 0 & 0 & \frac{EA}{L} & 0 & 0 \\ 0 & \frac{-12EI}{L^3} & \frac{-6EI}{L^2} & 0 & \frac{12EI}{L^3} & \frac{-6EI}{L^2} \\ 0 & \frac{6EI}{L^2} & \frac{2EI}{L} & 0 & \frac{-6EI}{L^2} & \frac{4EI}{L} \end{bmatrix}\end{split}\]
These matrices are used directly when calling the BeamMatrices function within the assembly process.
Step 4: global assembly#
The last step is to compute the global matrices and the global forcing vector. We start by initializing the global matrices as 1-dimensional arrays.
nDof = 3*nNode # 3 Degrees of freedom per nodeK = np.zeros((nDof*nDof))M = np.zeros((nDof*nDof))
Then we loop over elements and perform all the elemental operations.
from BeamMatrices import BeamMatricesJacketfor iElem in np.arange(0, nElem): # Get the nodes of the elements NodeLeft = Elements[iElem][0]-1 NodeRight = Elements[iElem][1]-1 # Get the degrees of freedom that correspond to each node Dofs_Left = 3*(NodeLeft) + np.arange(0, 3) Dofs_Right = 3*(NodeRight) + np.arange(0, 3) # Get the properties of the element m = Elements[iElem][2] EA = Elements[iElem][3] EI = Elements[iElem][4] # Calculate the matrices of the element Me, Ke = BeamMatricesJacket(m, EA, EI, ([NodeCoord[NodeLeft][0], NodeCoord[NodeLeft][1]], [NodeCoord[NodeRight][0], NodeCoord[NodeRight][1]])) # Assemble the matrices at the correct place nodes = np.append(Dofs_Left, Dofs_Right) for i in np.arange(0, 6): for j in np.arange(0, 6): ij = nodes[j] + nodes[i]*nDof M[ij] = M[ij] + Me[i, j] K[ij] = K[ij] + Ke[i, j] # Reshape the global matrix from a 1-dimensional array to a 2-dimensional arrayM = M.reshape((nDof, nDof))K = K.reshape((nDof, nDof))
Now we have the global mass and stiffness matrices. However, in this example we have an additional point mass at the top corresponding to the nacelle. Then, we need to account for this mass adding its value at the corresponding DOFs, in this case the corresponding horizontal and vertical displacements associated to the top node.
nacelle_node = nNodenacelle_dof_h = 3*(nacelle_node-1)nacelle_dof_v = 3*(nacelle_node-1)+1M[nacelle_dof_h, nacelle_dof_h] += Mass_NacelleM[nacelle_dof_v, nacelle_dof_v] += Mass_Nacelle
That completes the filling of the matrices. Let’s have a look at the matrices’ structure.
# Look at the matrix structureplt.figure()plt.spy(M)plt.title("Mass matrix")plt.figure()plt.spy(K)plt.title("Stiffness matrix");
To apply the boundary conditions, we will remove the rows associated to the fixed DOFs and add the contribution to the right-hand-side. First, we obtain the free and fixed DOFs.
DofsP = np.arange(0, 6) # prescribed DOFsDofsF = np.arange(0, nDof) # free DOFsDofsF = np.delete(DofsF, DofsP) # remove the fixed DOFs from the free DOFs array# free & fixed array indicesfx = DofsF[:, np.newaxis]fy = DofsF[np.newaxis, :]bx = DofsP[:, np.newaxis]by = DofsP[np.newaxis, :]
We can re-order the matrices and vectors in blocks, such that it’s easy to operate with the blocks corresponding with the fixed DOFs.
# MassM_FF = M[fx, fy]M_FP = M[fx, by]M_PF = M[bx, fy]M_PP = M[bx, by]# StiffnessK_FF = K[fx, fy]K_FP = K[fx, by]K_PF = K[bx, fy]K_PP = K[bx, by]
Part 1: Full FEM#
We can also solve the final ODE. In this example we assume that we don’t have an external force and the system is excited by the effect of an earthquak, modelled with the following horizontal motion at the boundary nodes:
\[ u_P = A_0 \sin (2\pi f_0t) \]
\[ u'_P = 2 \pi f_0 A_0 \cos(2\pi f_0t) \]
\[ u''_P = - (2 \pi f_0)^2 A_0 \sin(2\pi f_0t) \]
FEM_StartTime = time.time()# Set Dimensions and initial conditions of state vectornDofsF = len(DofsF)udofs = np.arange(0, nDofsF)vdofs = np.arange(nDofsF, 2*nDofsF)q0 = np.zeros((2*nDofsF))# Define the boundary conditionsdef ub(t, T0): if t <= T0: return A0*np.sin(2*np.pi*f0*t)*np.array([1, 0, 0, 1, 0, 0]) else: return np.array([0, 0, 0, 0, 0, 0])def dub_dt2(t, T0): return -(2*np.pi*f0)**2*ub(t, T0)# Time span (output purposes)tf = 5tspan = np.arange(0, tf, 0.01)# Solvedef odefun(t, q): return np.append(q[vdofs], np.linalg.solve(M_FF, -np.dot(K_FP, ub(t, T0)) - np.dot(M_FP, dub_dt2(t, T0)) - np.dot(K_FF, q[udofs]))).tolist()sol = scpi.solve_ivp(fun=odefun, t_span=[tspan[0], tspan[-1]], y0=q0, t_eval=tspan)
# Plot resultsu_total = np.zeros((len(tspan), nDof))u_total[:, DofsF[0]:DofsF[-1]+1] = np.transpose(sol.y[udofs[0]:udofs[-1]+1, :])for it in np.arange(0, len(tspan)): u_total[it, DofsP] = ub(tspan[it], T0)nacelle_result = u_total[:, 60:63]fig, axs = plt.subplots(1, 3, figsize = (10, 3))axs[0].plot(sol.t, nacelle_result[:, 0])axs[0].set_title("u$_h$")axs[0].set_xlabel("Time [s]")axs[0].set_ylabel("u$_h$ [m]")axs[1].plot(sol.t, nacelle_result[:, 1])axs[1].set_title("u$_v$")axs[1].set_xlabel("Time [s]")axs[1].set_ylabel("u$_v$ [m]")axs[2].plot(sol.t, nacelle_result[:, 2])axs[2].set_title("$\phi_{top}$")axs[2].set_xlabel("Time[s]")axs[2].set_ylabel("$\phi_{top}$ [rad]")fig.tight_layout();FEM_EndTime = time.time()# Separated for honest time evaluationsFEM_time = sol.tFEM_uh = nacelle_result[:,0]
Part 2: Model superposition#
Step 5: Solving the ODE#
Here we divert from tutorial 4.3.
Using the matrices associated to the free DOFs, we can perform a modal analysis to get more information on how the structure will deform and determine the natural frequencies.
\[ ( \boldsymbol{K}_{FF} - \omega^2 \boldsymbol{M}_{FF} ) \boldsymbol{\phi} = \boldsymbol{0} \]
To compute the natural frequencies and mode shapes we use the eig command, which is part of the NumPy package. For more information see: https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html
MOD_StartTime = time.time()mat = np.dot(np.linalg.inv(M_FF), K_FF)w2, vr = np.linalg.eig(mat)w = np.sqrt(w2.real)f = w/2/np.pi
We sort the frequencies and mode shapes in descending order:
idx = f.argsort()#[::-1] f = f[idx]vr_sorted = vr[:,idx]ModalShape = np.zeros((nDof, len(f)))ModalShape[6:, :] = vr_sorted
Let’s see what these modes look like. Here, we plot the first 9 modes. Note that the system will have 19 x 3 = 57 modes (as many as the discrete system DOFs).
nMode = 25plt.figure()nCol = int(np.round(np.sqrt(nMode)))nRow = int(np.ceil(nMode/nCol))for iMode in np.arange(1, nMode + 1): plt.subplot(nRow, nCol, iMode) Shape = ModalShape[:, iMode-1] # Scale the mode such that maximum deformation is 10 MaxTranslationx = np.max(np.abs(Shape[0::3])) MaxTranslationy = np.max(np.abs(Shape[1::3])) Shape[0::3] = Shape[0::3]/MaxTranslationx*10 Shape[1::3] = Shape[1::3]/MaxTranslationy*10 # Get the deformed shape DisplacedNode = ([i[0] for i in NodeCoord] + Shape[0::3], [i[1] for i in NodeCoord] + Shape[1::3]) for iElem in np.arange(0, nElem): NodeLeft = Elements[iElem][0]-1 NodeRight = Elements[iElem][1]-1 m = Elements[iElem][2] EA = Elements[iElem][3] EI = Elements[iElem][4] if m == m_Jacket and EA == EA_Jacket and EI == EI_Jacket: c = 'b' elif m == m_Pile and EA == EA_Pile and EI == EI_Pile: c = 'r' else: print("ERROR: unknown material. Check your inputs.") break plt.plot([DisplacedNode[0][NodeLeft], DisplacedNode[0][NodeRight]], [DisplacedNode[1][NodeLeft], DisplacedNode[1][NodeRight]], c=c) plt.title("Mode "+str(iMode)+": f = "+str(round(f[iMode]))+" Hz") plt.axis('equal')# automatically fix subplot spacingplt.tight_layout()
Step 6: calculate modal quantities#
Now we want to compute the mass and stiffness coefficients that go on the diagonal of the respective mode, \(m_{ii}\) and \(k_{ii}\). Here, we also compute the contribution from the external force, inluding the effect of the boundary conditions, \(f_i^{eq}\).
\[\begin{split} \begin{bmatrix} m_{11} \ddot{\Gamma}_1 \\ m_{22} \ddot{\Gamma}_2 \\ \vdots \\ m_{NN} \ddot{\Gamma}_N \end{bmatrix} + \begin{bmatrix} k_{11} \Gamma_1 \\ k_{22} \Gamma_2 \\ \vdots \\ k_{NN} \Gamma_N \end{bmatrix} = \begin{bmatrix} f_1^{eq} \\ f_2^{eq} \\ \vdots \\ f_N^{eq} \end{bmatrix} \rightarrow \begin{cases} m_{11}\ddot{\Gamma}_1 + k_{11}\Gamma_1 = f_1^{eq} \\ m_{22}\ddot{\Gamma}_2 + k_{22}\Gamma_2 = f_2^{eq} \\ \qquad \qquad \vdots \\ m_{NN}\ddot{\Gamma}_N + k_{NN}\Gamma_N = f_N^{eq} \end{cases} \end{split}\]
Additionally, we also add damping into the system. In that case we use the modal damping ratio approach, with a damping ratio of 1% (\(\xi_j=0.01\)).
\[ c_{jj} = \xi_j c_{cr,j} = 2 \xi_j \sqrt{m_{jj}k_{jj}} \]
This will lead to a final system:
\[\begin{split} \begin{cases} m_{11}\ddot{\Gamma}_1 + c_{11}\dot{\Gamma}_1 + k_{11}\Gamma_1 = f_1^{eq} \\ m_{22}\ddot{\Gamma}_2 + c_{22}\dot{\Gamma}_2 + k_{22}\Gamma_2 = f_2^{eq} \\ \qquad \qquad \vdots \\ m_{NN}\ddot{\Gamma}_N + c_{NN}\dot{\Gamma}_N + k_{NN}\Gamma_N = f_N^{eq} \end{cases} \end{split}\]
Attention: we are calculating all the modes, but in general this shouldn’t be the case. You should only calculate the response for these modes that you want to consider. Here, we are calculating all modes, because we will compare the “reduced” and the “full” response. To consider only a few, you should replace nMode by the number of modes you wish to consider.
PHI = vr_sorted[:,0:nMode]Mm = np.zeros(nMode)Km = np.zeros(nMode)Cm = np.zeros(nMode)ModalDampRatio = 0.00# Compute your "nMode" entries of the modal mass, stiffness and dampingfor iMode in np.arange(0,nMode): print('Computing Mode: ',iMode+1) # Starts at 0 off course Mm[iMode] = PHI[:,iMode].T @ M_FF @ PHI[:,iMode] Km[iMode] = PHI[:,iMode].T @ K_FF @ PHI[:,iMode] Cm[iMode] = 2*ModalDampRatio*np.sqrt(Mm[iMode]*Km[iMode]) print('Mm = ',Mm[iMode],', Km = ', Km[iMode],', Cm = ', Cm[iMode])
Computing Mode: 1Mm = 9770.858374615917 , Km = 35.155971711622826 , Cm = 0.0Computing Mode: 2Mm = 5878.202857147078 , Km = 870.3127554379188 , Cm = 0.0Computing Mode: 3Mm = 8174.205181930304 , Km = 1943.0086801171649 , Cm = 0.0Computing Mode: 4Mm = 4775.304243275167 , Km = 9370.730498323164 , Cm = 0.0Computing Mode: 5Mm = 5303.990003284622 , Km = 24397.85071457471 , Cm = 0.0Computing Mode: 6Mm = 4412.846973208809 , Km = 31345.132214564914 , Cm = 0.0Computing Mode: 7Mm = 4040.488162793453 , Km = 47380.65684652576 , Cm = 0.0Computing Mode: 8Mm = 3147.261174569589 , Km = 73548.70340845389 , Cm = 0.0Computing Mode: 9Mm = 3778.9547524130057 , Km = 98975.85673442972 , Cm = 0.0Computing Mode: 10Mm = 3669.3872240330656 , Km = 106568.72904968809 , Cm = 0.0Computing Mode: 11Mm = 3563.95580730377 , Km = 126372.67620113312 , Cm = 0.0Computing Mode: 12Mm = 2999.515033090926 , Km = 138412.5005604382 , Cm = 0.0Computing Mode: 13Mm = 2785.67340049797 , Km = 145600.24421639566 , Cm = 0.0Computing Mode: 14Mm = 2164.8498020365623 , Km = 117410.75951904256 , Cm = 0.0Computing Mode: 15Mm = 2950.421152039472 , Km = 194835.9266121073 , Cm = 0.0Computing Mode: 16Mm = 2271.0953866874515 , Km = 159111.38789952552 , Cm = 0.0Computing Mode: 17Mm = 3076.906852717982 , Km = 255183.07893519773 , Cm = 0.0Computing Mode: 18Mm = 5459.4254369148775 , Km = 488524.059174493 , Cm = 0.0Computing Mode: 19Mm = 2655.684990242143 , Km = 253428.1051428055 , Cm = 0.0Computing Mode: 20Mm = 2093.9658414421606 , Km = 257422.39952900942 , Cm = 0.0Computing Mode: 21Mm = 2395.233999978137 , Km = 348012.8576186892 , Cm = 0.0Computing Mode: 22Mm = 2944.0216626312504 , Km = 505716.88895826286 , Cm = 0.0Computing Mode: 23Mm = 1430.9733626700895 , Km = 258880.45076293507 , Cm = 0.0Computing Mode: 24Mm = 1067.7174249898526 , Km = 227509.665650035 , Cm = 0.0Computing Mode: 25Mm = 1926.0621426582354 , Km = 421993.1210081142 , Cm = 0.0
# Define the boundary conditions (same as before)def ub(t, T0): if t <= T0: return A0*np.sin(2*np.pi*f0*t)*np.array([1, 0, 0, 1, 0, 0]) else: return np.array([0, 0, 0, 0, 0, 0])def dub_dt2(t, T0): return -(2*np.pi*f0)**2*ub(t, T0)def F(t): return -PHI.T @ ( K_FP @ ub(t,T0) + M_FP @ dub_dt2(t,T0) )
# Solve the resulting ODE:def qdot(t,q): Um = q[0:nMode] Vm = q[nMode:2*nMode] Am = ( F(t) - (Km * Um + Cm * Vm) ) / Mm return np.append(Vm,Am)q0 = np.zeros(2*nMode)q = scpi.solve_ivp(fun=qdot,y0=q0,t_span=[tspan[0], tspan[-1]], t_eval=tspan)
# Show the result on the original coordinatesU_F = np.zeros((len(q.t),len(PHI[:,0])))for iMode in np.arange(0,nMode): for it in np.arange(0,len(q.t)): U_F[it,:] += PHI[:,iMode] * q.y[iMode][it] plt.plot(tspan, U_F[:,60-6])plt.xlabel("Time [s]")plt.ylabel("Excursion [m]")plt.title("FEM results from modal analysis");MOD_EndTime = time.time()# Separated for honest time evaluationsMOD_time = q.tMOD_uh = U_F[:,60-6]#MOD_uh = FEM_uh/1.06 # testing
U_F[2]
array([ 4.03467610e-05, 8.19064689e-06, -2.39890237e-07, 6.13652366e-05, 4.71957294e-05, 1.38780849e-05, -9.41204572e-06, 2.22252682e-05, -7.75709145e-06, -2.15863816e-06, 1.29190530e-06, 5.52767593e-07, -6.72413804e-06, -7.42796072e-06, 1.84252939e-06, -1.52470406e-06, 3.94100411e-06, 1.04998073e-06, 7.81377252e-06, 2.20155008e-06, -5.16035042e-07, -6.12361809e-07, -2.15970219e-08, 3.36260327e-08, 7.39020963e-06, -3.37525952e-06, -6.09177384e-07, -7.92303654e-07, 5.73035930e-10, -1.78372774e-08, -2.18187033e-07, 2.19277880e-08, -8.40869558e-08, 5.11433442e-07, 3.78425393e-08, -4.37540619e-08, 4.76176347e-07, 4.47852790e-08, 4.81804372e-08, -1.97904293e-07, 4.12338201e-08, 6.80239540e-08, -5.58955777e-07, 2.81024637e-08, -4.78310089e-09, -1.41710909e-07, 8.56704172e-09, -6.65727901e-08, 4.37440404e-07, -1.26713353e-08, -3.30225674e-08, 3.40137874e-07, -3.05039967e-08, 5.17999288e-08, -4.22561354e-07, -4.06401514e-08, 8.80208095e-08])
Part 3: Performance comparison#
Here we will compare the two methods. We are only comparing step 5 and 6, since step 1 through 4 are identical. While there are different demands for the plots, the general comparison idea can still be carried out.
We can see that:
The modal superposition variant performs 10 times faster
The accuracy xxx …. The ‘real’ answer however can never be known due to the simplifications in both models. Due to its internal assumptions modal superposition will theoretically get less and less accurate for increased damping. In jacket structures however this is mostly not significant.
plt.plot(FEM_time, FEM_uh, label=f"FEM ({FEM_EndTime-FEM_StartTime:.2f} seconds)")plt.plot(MOD_time, MOD_uh, label=f"MOD ({MOD_EndTime-MOD_StartTime:.2f} seconds)")plt.title("FEM to MOD result comparison")plt.legend(loc=3)plt.xlabel("Time [s]")plt.ylabel("Nacelle horizontal displacement [m]");
