diff --git a/.calkit/overleaf-sync.json b/.calkit/overleaf-sync.json index ad76b020..96ec707f 100644 --- a/.calkit/overleaf-sync.json +++ b/.calkit/overleaf-sync.json @@ -9,7 +9,7 @@ }, "pubs/applied-ocean-research-model": { "project_id": "69ecde8851d14790c59f03d5", - "last_sync_commit": "07a22aeee8dc9f531a0e62af5ffe54e15e8a6c82" + "last_sync_commit": "32c12bee2947430417b247513669c0738085fcde" }, "pubs/defense": { "project_id": "69ef01941916c106f7f5c565", diff --git a/pubs/applied-ocean-research-model/commands-aor.tex b/pubs/applied-ocean-research-model/commands-aor.tex index d5e63994..7fd3abaf 100644 --- a/pubs/applied-ocean-research-model/commands-aor.tex +++ b/pubs/applied-ocean-research-model/commands-aor.tex @@ -4,7 +4,41 @@ \usepackage{microtype} \usepackage{subcaption} % for subfigures +\usepackage{enumitem} +\usepackage{tikz} +%\usepackage[section]{placeins} % prevent figs from going in wrong section +\newcolumntype{P}[1]{>{\centering\arraybackslash}p{#1}} +\newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}} + +\newcommand{\singleColMacro}[1]{ + \ifdefined\DISSERATION + \else + \onecolumn + \fi + {#1} + \ifdefined\DISSERATION + \else + \twocolumn + \fi +} +\newcommand{\journal}{Applied Ocean Research} + +% Source - https://tex.stackexchange.com/a/515729 +% Posted by user194703, modified by community. See post 'Timeline' for change history +% Retrieved 2026-06-21, License - CC BY-SA 4.0 +\usetikzlibrary{shapes.geometric} +\newcommand{\Stars}[2][fill=black,draw=black]{ + \begin{tikzpicture}[baseline=-0.35em,#1] + \foreach \X in {1,...,3} + {\pgfmathsetmacro{\xfill}{min(1,max(1+#2-\X,0))} + \path (\X*1.1em,0) + node[star,draw,star point height=0.25em,minimum size=1em,inner sep=0pt, + path picture={\fill (path picture bounding box.south west) + rectangle ([xshift=\xfill*1em]path picture bounding box.north west);}]{}; + } + \end{tikzpicture} +} %\usepackage{graphicx} %\usepackage{float} %\usepackage{placeins} diff --git a/pubs/applied-ocean-research-model/main.tex b/pubs/applied-ocean-research-model/main.tex index f7616a0b..acc27552 100644 --- a/pubs/applied-ocean-research-model/main.tex +++ b/pubs/applied-ocean-research-model/main.tex @@ -99,18 +99,13 @@ \end{highlights} \begin{keywords} -wave energy conversion \sep marine renewable energy \sep semi-analytical hydrodynamics \sep -multi-port circuit \sep -constrained optimal control \sep -structural design \sep -survivability analysis \sep +linearized pseudo-spectral optimal control \sep +structural survivability \sep techno-economic modeling \sep -two-body point absorber \sep -power take-off modeling \sep -validation \sep -benchmarking +model validation \sep +computational benchmarking \end{keywords} \maketitle @@ -129,9 +124,9 @@ \section*{Acknowledgements} -The authors thank Kapil Khanal, En Lo, Yinghui Bimali, and John Fernandez for assistance with hydrodynamics; +The authors thank Kapil Khanal, Yinghui Bimali, En Lo, and John Fernandez for assistance with hydrodynamics; Fabien Royer for guidance on structures; -and Ryan Coe, Jacob Mays, Patrick Reed, Nate DeGeode, and Alaa Ahmed for providing valuable manuscript feedback. +Ryan Coe, Jacob Mays, Patrick Reed, Nate DeGeode, and Alaa Ahmed for technical feedback on a draft manuscript; and Nola McCabe for proofreading support. R.M. acknowledges funding from the National Science Foundation Graduate Research Fellowship. M.D. acknowledges funding from the Fund for Undergraduate Research on Solutions to Climate Change and the Bill Nye ’77 Award in Undergraduate Research. @@ -139,7 +134,7 @@ \section*{Acknowledgements} This material is based on work supported by the National Science Foundation Graduate Research Fellowship under Grant No.~DGE--2139899. Any opinions, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. -The graphical abstract uses icons from \url{www.flaticon.com} in accordance with the Flaticon license. +%The graphical abstract uses icons from \url{www.flaticon.com} in accordance with the Flaticon license. \section*{Data availability statement} @@ -151,9 +146,9 @@ \section*{Data availability statement} and accessed at \url{https://calkit.io/symbiotic-engineering/mdocean}. \appendix - +\singleColMacro{ \include{appendices} - +} \printcredits \bibliographystyle{cas-model2-names} diff --git a/pubs/applied-ocean-research-model/sections/benchmarking.tex b/pubs/applied-ocean-research-model/sections/benchmarking.tex index e4e08d5f..cfcc5538 100644 --- a/pubs/applied-ocean-research-model/sections/benchmarking.tex +++ b/pubs/applied-ocean-research-model/sections/benchmarking.tex @@ -21,14 +21,12 @@ \subsection{Dynamic Validation Using WEC-Sim}\label{sec:dynamic-validation} \fi The absolute error in average power compared to the WEC-Sim power is less than \resultsAOR[wecsimAvgPowerErrorBestCase] in the best case and \resultsAOR[wecsimAvgPowerErrorWorstCase] in the worst case, with an error breakdown for all simulation scenarios and sea states shown in \Cref{fig:error-histogram}. - -\begin{figure}[htbp] +\begin{figure*}[htbp] \centering \includegraphics[width=1\linewidth]{figs/from-matlab/wecsim_wcsm_multi__histogram.pdf} - \caption{Error breakdown based on WEC-Sim Validation Runs} + \caption{Error breakdown based on WEC-Sim Verification Runs} \label{fig:error-histogram} -\end{figure} - +\end{figure*} The detailed error breakdown across drag-on/drag-off and MEEM/WAMIT coefficient configurations is provided in \Cref{sec:appendix-dynamic-validation}, revealing that the dominant error sources are interactions between drag, hydrodynamic-coefficient mismatch, and the inter-body phase relationship in the 2-DOF model. \Cref{sec:appendix-dynamic-validation} also validates the describing-function approximation itself, showing total harmonic distortion below 1\% in the worst sea state and excellent agreement between the assumed and actual drag force waveforms at all four corners of the JPD. @@ -107,12 +105,12 @@ \subsection{Static Validation} % \label{tab:validation} % \end{table} -\begin{table}[htbp] -\centering -\input{tables/from-matlab/validation.tex} -\caption{Validation} -\label{tab:validation} -\end{table} + \begin{table*}[htbp] + \centering + \input{tables/from-matlab/validation.tex} + \caption{Verification} + \label{tab:validation} + \end{table*} %\hl{Explain sources of error and rough uncertainty and the implications of what we can trust} @@ -136,13 +134,13 @@ \subsection{Runtime Benchmarking} \begin{figure}[b!] \centering -\includegraphics[width=0.5\linewidth]{figs/from-matlab/sim_runtime.pdf} +\includegraphics[width=\linewidth]{figs/from-matlab/sim_runtime.pdf} \caption{Bar chart showing simulation runtime breakdown between modules}\label{fig:runtime-modules} \end{figure} \begin{figure}[t!] \centering -\includegraphics[width=0.5\linewidth]{figs/from-matlab/hydro_runtime_logscale.pdf} +\includegraphics[width=\linewidth]{figs/from-matlab/hydro_runtime_logscale.pdf} \caption{Bar chart demonstrating the speed improvement of MDOcean's hydro module over baseline solver Capytaine}\label{fig:runtime-hydro} \end{figure} @@ -154,7 +152,7 @@ \subsection{Runtime Benchmarking} \begin{figure}[t!] \centering -\includegraphics[width=0.5\linewidth]{figs/from-matlab/dynamics_runtime.pdf} +\includegraphics[width=\linewidth]{figs/from-matlab/dynamics_runtime.pdf} \caption{Bar chart demonstrating the speed improvement of MDOcean's dynamics module over baseline solver WEC-Sim}\label{fig:runtime-dynamics} \end{figure} diff --git a/pubs/applied-ocean-research-model/sections/discussion.tex b/pubs/applied-ocean-research-model/sections/discussion.tex index dfd48c3c..7ff062c5 100644 --- a/pubs/applied-ocean-research-model/sections/discussion.tex +++ b/pubs/applied-ocean-research-model/sections/discussion.tex @@ -34,7 +34,7 @@ \subsubsection{Effect of Bulk Dimensions on Hydrodynamic Efficiency} \fi \paragraph{Best Design as a function of Wave Environment} -\begin{figure}[htbp] +\begin{figure*}[htbp] \centering \begin{subfigure}[t]{0.48\linewidth} \includegraphics[width=\linewidth]{figs/from-matlab/sweep_geoms_line.pdf} @@ -49,7 +49,7 @@ \subsubsection{Effect of Bulk Dimensions on Hydrodynamic Efficiency} \end{subfigure} \caption{Effect of wave environment and hydrodynamic design variables on (a) radiation efficiency and (b) capture width ratio} \label{fig:meem-sweep-m0h} -\end{figure} +\end{figure*} \ifdefined\DISSERTATION Starting with \Cref{fig:meem-sweep-eff}, we observe that the most efficient design depends on the frequency-depth regime $m_0h$ of the wave environment, shown on the x-axis. @@ -77,7 +77,7 @@ \subsubsection{Effect of Bulk Dimensions on Hydrodynamic Efficiency} \begin{figure}[htbp] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/sweep_geoms_pareto_nondim.pdf} + \includegraphics[width=.9\linewidth]{figs/from-matlab/sweep_geoms_pareto_nondim.pdf} \caption{Effect of hydrodynamic design variables on radiation efficiency and power per unit surface area} \label{fig:meem-sweep-pareto} \end{figure} @@ -88,7 +88,7 @@ \subsubsection{Damping Plate Size} This is because the same force has a lower lever arm to the column and therefore creates less bending moment. \begin{figure}[htbp] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/damping_plate_aspect_ratio.pdf} + \includegraphics[width=.85\linewidth]{figs/from-matlab/damping_plate_aspect_ratio.pdf} \caption{Normalized effect of damping plate aspect ratio on maximum stress and deflection. The dashed lines indicate the nominal design point at $b/a = 0.2$.} \label{fig:damping-plate-maxs} @@ -139,7 +139,7 @@ \subsubsection{Effect of PTO Force and Power Limits} \begin{figure}[htbp] \centering -\includegraphics[width=.8\linewidth]{figs/from-matlab/pto_sweep.pdf} +\includegraphics[width=.9\linewidth]{figs/from-matlab/pto_sweep.pdf} \caption{Effect of Generator Force Limit and Power Limit on Annual Average Power and LCOE} \label{fig:force-power-limit-sweep} \end{figure} @@ -153,7 +153,7 @@ \subsubsection{Effect of PTO Force and Power Limits} The theoretical analysis in \Cref{sec:appendix-constraint-sensitivity} derives mathematical conditions for when this is the case. Convexity suggests that local optimization methods should be effective for finding the optimal PTO design. \else - Power and LCOE both appear to be convex with respect to the force and power limits across the feasible region, a structural property exploited by the optimization strategy in the companion paper \citep{mccabe_leveraging_2026}. + Power and LCOE both appear to be convex with respect to the force and power limits across the feasible region, a structural property that future optimizations can exploit. \fi %%%%%%%%%%%%%%%%%%% @@ -172,7 +172,7 @@ \subsubsection{Design Space Exploration} \begin{figure}[htbp] \centering -\includegraphics[width=\linewidth]{figs/from-matlab/experiments_ratios.pdf} +\includegraphics[width=.9\linewidth]{figs/from-matlab/experiments_ratios.pdf} \caption{Design of experiments}\label{fig:experiments} \end{figure} @@ -211,6 +211,8 @@ \subsection{Multidisciplinary Insights} This multiplicative decomposition allows isolation of the effect of each design variable and parameter on power, intuitively identifying which design variables are most important for improving power. \Cref{tab:power-matrix-dependence} maps the dependence explicitly. +{ +\renewcommand{\arraystretch}{1.35} \begin{table}[h] \centering \caption{Dependence of efficiencies on inputs}\label{tab:power-matrix-dependence} @@ -225,6 +227,7 @@ \subsection{Multidisciplinary Insights} \hline \end{tabular} \end{table} +} \ifdefined\DISSERTATION Importantly, the later matrices in the product depend on the design variables and parameters that affect the earlier matrices. @@ -243,12 +246,12 @@ \subsection{Multidisciplinary Insights} Low-period sea states achieve higher radiation efficiency, while high-period sea states have lower radiation efficiency but contribute significantly to total production due to their higher energy content; despite the most probable period being 6--8~s, 10~s waves contribute the most to annual energy production \citep{zou_practical_2023}. \fi -\begin{figure}[htbp] +\begin{figure*}[htbp] \centering -\includegraphics[width=\linewidth]{figs/from-matlab/nominal_power_matrix.pdf} +\includegraphics[width=.95\linewidth]{figs/from-matlab/nominal_power_matrix.pdf} \caption{Power matrix decomposition} \label{fig:power-matrix-decomposition} -\end{figure} +\end{figure*} The effect of drag is most significant in the 11-12 second range, corresponding to the spar's natural frequency. \ifdefined\DISSERTATION Interestingly, this is also where the force limit has the strongest effect, indicating that the large amplitudes at the natural frequency drive up the force more than the larger stiffnesses required for reactive control at frequencies far from the natural frequency. @@ -278,19 +281,21 @@ \subsection{Multidisciplinary Insights} The result is: \begin{equation}\label{eq:power-double-sum-main} \begin{aligned} - \overline{P}_{elec} = \overline{P}_{elec,0} + \overline{P}_{elec} =~&\overline{P}_{elec,0} - \eta \sum_{\mu\nu\beta} - &\left[ + \left[ \sqrt{\textrm{aff}_{\mu\nu\beta}(b_{\mu})} - + \textrm{quad}_{\mu\nu\beta}(b_{\mu}, b_{\nu}) + ~+ \right. \\ + & \left. \textrm{quad}_{\mu\nu\beta}(b_{\mu}, b_{\nu}) \right.+ \\ &\left. \left(\textrm{aff}_{\mu\nu\beta}(b_{\mu}) + \textrm{aff}_{\mu\nu\beta}(b_{\nu})\right) - \sqrt{\textrm{quad}_{\mu\nu\beta}(b_{\mu}, b_{\nu})} + \sqrt{\textrm{quad}_{\mu\nu\beta}(b_{\mu}, b_{\nu})}~ \right] \end{aligned}\end{equation} where $\overline{P}_{\text{elec},0}$ is the unconstrained power and $\textrm{aff}_{\mu\nu\beta}$, $\textrm{quad}_{\mu\nu\beta}$ are affine and quadratic functions of $b_\mu$ and $b_\nu$, respectively. +$\mu$ and $\nu$ are constraint indices, expanded in \Cref{sec:appendix-constraint-sensitivity}. Under conditions derived in \Cref{sec:appendix-constraint-sensitivity}, this scaling law is convex in the constraint coefficients $b_\mu$. While satisfaction of these conditions is not guaranteed, convexity is observed to hold across the PTO sweep of \Cref{fig:force-power-limit-sweep}. @@ -328,12 +333,12 @@ \subsection{Limitations and Future Work}\label{sec:unmodeled-effects} \Cref{tab:future-work} summarizes potential future improvements to the model, distinguishing between model improvements that would enhance the accuracy or realism of studies that can be conducted with the present model and those that would unlock the ability to answer design questions that the current model cannot. This section describes the relevance and possible implementation paths for each development area. \else - \Cref{tab:future-work} summarizes principal limitations and future work, distinguishing model improvements that would enhance the accuracy of currently-achievable studies from extensions that would unlock new design questions. + \Cref{tab:future-work} summarizes principal limitations and future work, distinguishing model improvements to enhance the accuracy of currently-achievable studies from extensions to unlock new design questions. MDOcean is open-source \citep{mccabe_mdocean_2024}; community contributions are welcome, and the authors are pursuing the multi-objective RM3 optimization in the companion paper \citep{mccabe_leveraging_2026}. \fi \newcommand{\modelTrustBuilders}{ - \begin{enumerate} + \begin{enumerate}[leftmargin=*] \item Surge force and mooring cost \item Nonlinear storm wave forces \item Irregular waves @@ -343,7 +348,7 @@ \subsection{Limitations and Future Work}\label{sec:unmodeled-effects} } \newcommand{\modelStudyEnablers}{ - \begin{enumerate} + \begin{enumerate}[leftmargin=*] \item Spectral power, load, amplitude \item Different WEC archetypes \item Lifetime and sea state contours @@ -357,8 +362,8 @@ \subsection{Limitations and Future Work}\label{sec:unmodeled-effects} \caption{Future model improvements} \label{tab:future-work} \begin{tabular}{ - >{\centering\arraybackslash}p{0.5\linewidth} - >{\centering\arraybackslash}p{0.5\linewidth}} + >{\raggedright\arraybackslash}m{0.35\linewidth} + >{\raggedright\arraybackslash}m{0.5\linewidth}} Enhance Trust in Achievable Studies & Unlock New Studies \\ \hline \modelTrustBuilders & \modelStudyEnablers \\ \end{tabular} @@ -451,7 +456,7 @@ \subsection{Limitations and Future Work}\label{sec:unmodeled-effects} Whether semi-analytical or numerical, extension to other archetypes would require significant development effort but unlock new comparative insights. \else The most consequential current limitations are: -\begin{itemize} +\begin{itemize}[leftmargin=*] \item \textbf{Neglect of surge force in the structures module and absence of a mooring cost model.} Surge forces on the nominal RM3 are \resultsAOR[surgeForceFloatNominal] (float) and \resultsAOR[surgeForceSparNominal] (spar), and mooring/foundation accounts for 12\% of CAPEX (\Cref{tab:CBS}); incorporating these would likely affect optimal designs. \item \textbf{Regular-wave assumption in storm load cases.} Storm waves are nonlinear, and the regular-wave equivalent cannot capture transient peaks; second-order MEEM \citep{cong_novel_2020,mavrakos_second-order_2009} or the slender-body approximation in RAFT \citep{carmo_slender-body_2025} are candidate extensions. \item \textbf{Regular-wave assumption in operational loads.} Stochastic linearization \citep{da_silva_statistical_2020,da_silva_stochastic_2023,kluger_synergistic_2017,folley_spectral-domain_2016,spanos_efficient_2016,neshat_enhancing_2024} could replace the describing function to handle spectral inputs and enable spectral fatigue, grid-integration, and storage-sizing analyses \citep{mccabe_wec_2025}. diff --git a/pubs/applied-ocean-research-model/sections/meem-appendix.tex b/pubs/applied-ocean-research-model/sections/meem-appendix.tex index e77276c0..6b448ac6 100644 --- a/pubs/applied-ocean-research-model/sections/meem-appendix.tex +++ b/pubs/applied-ocean-research-model/sections/meem-appendix.tex @@ -21,7 +21,7 @@ \subsection{Linear Hydrodynamics and Eigenfunctions} Extension to many regions is discussed in section 3.4. \Cref{fig:meem-regions} illustrates the regions and dimensions: two internal regions \textit{i1} and \textit{i2}, and an external region \textit{e} extending to infinity. -\begin{figure}[htbp] +\begin{figure}[htb] \centering \includegraphics{figs/from-matlab/meem_dims.pdf} \caption{Fluid regions and dimensions used in the dual concentric cylinder MEEM} @@ -34,10 +34,14 @@ \subsection{Linear Hydrodynamics and Eigenfunctions} By construction, this potential obeys all boundary conditions except for zero radial velocity on radial body surfaces. The unknown coefficients must be computed to enforce this final condition as well as continuity across regions, which will be the subject of \Cref{sec:appendix-meem-matching}. -\begin{landscape} \begin{table}[htbp] \centering - \begin{tabular}{|>{\centering\arraybackslash}p{0.085\linewidth}|>{\centering\arraybackslash}p{0.26\linewidth}|>{\centering\arraybackslash}p{0.34\linewidth}|>{\centering\arraybackslash}p{0.32\linewidth}|} \hline + \begin{tabular}{| + M{0.069\linewidth}| + M{0.29\linewidth}| + M{0.29\linewidth}| + M{0.25\linewidth}| + } \hline Region& $i1$& $i2$& $e$\\ \hline Homog. potential $\phi_h(r,z)$& $\displaystyle\sum_n C_{1n}^{i1} R_{1n}^{i1}(r) Z_n^{i1}(z)$& $\displaystyle\sum_m \left(C_{1m}^{i2} R_{1m}^{i2}(r) + C_{2m}^{i2} R_{2m}^{i2}(r) \right) Z_m^{i2}(z)$& $\displaystyle\sum_k B_k^e \Lambda_k(r) Z_k^e(z)$\\ \hline Partic. potential $\phi_p(r,z)$& $\displaystyle\frac{1}{2(h-d_1)}\left[ (z+h)^2 - \frac{r^2}{2}\right] $& $\displaystyle\frac{1}{2(h-d_2)}\left[ (z+h)^2 - \frac{r^2}{2}\right]$& $0$\\ \hline @@ -69,9 +73,7 @@ \subsection{Linear Hydrodynamics and Eigenfunctions} \end{tabular} \caption{Equations for potential (homogeneous and particular) and eigenfunctions (radial and vertical) for each region.} \label{tab:MEEM-eigenfunctions} - \fillandplacepagenumber \end{table} -\end{landscape} In \Cref{tab:MEEM-eigenfunctions}, $\textrm{I}_0$, $\textrm{K}_0$, and $\textrm{H}_0^1$ are different Bessel functions of order zero; 1M and 2M mean body 1 and 2 (spar and float) are moving respectively, while 1S and 2S mean each is stationary; and the $N_k$ expression is defined as: @@ -269,7 +271,6 @@ \subsection{Linear Hydrodynamics and Eigenfunctions} \Cref{tab:meem-integrals} lists several integrals of the radial and vertical eigenfunctions, $\boldsymbol{\mathcal{R}}$ and $\boldsymbol{\mathcal{Z}}$ respectively, that will be needed in the calculations to follow. -\begin{landscape} \begin{table}[htbp] \centering \caption{Eigenfunction integrals} @@ -283,9 +284,8 @@ \subsection{Linear Hydrodynamics and Eigenfunctions} \multirow{2}{*}{\ZmkDefn} & $k=0$ & \ZmZeroKZero & \ZmOneKZero \\ \cline{2-4} & $k\geq1$ & \ZkOneMZero & \ZkOneMOne \\ \hline \end{tabular} - \fillandplacepagenumber \end{table} -\end{landscape} + \subsection{Matching Across Fluid Boundaries} \label{sec:appendix-meem-matching} @@ -333,12 +333,11 @@ \subsection{Block Matrix Structure} They are written in compact notation using row vectors of basis functions, so $\vec{R_{1n}^{i1}}=[R_{10}^{i1}, R_{11}^{i1}, ..., R_{1(N-1)}^{i1}]$ and so on. Each basis function is evaluated at the radius described to the left of its row in the table. $0_{ij}$ and $1_{ij}$ are the $i$ x $j$ matrices of zeros and ones respectively; diag($\cdot$) constructs a diagonal matrix from a vector; and $\odot$ is the Hadamard (element-wise) product. -\begin{landscape} \begin{table}[htbp] \centering \caption{MEEM A-matrix} \label{tab:MEEM-A-matrix} - \begin{tabular}{|>{\centering\arraybackslash}p{0.18\linewidth}|c||c|c|c|c|}\hline + \begin{tabular}{|M{0.12\linewidth}|c||c|c|c|c|}\hline & & $\vec{C}_{1n}^{i1}$& $\vec{C}_{1m}^{i2}$& $\vec{C}_{2m}^{i2}$&$\vec{B}_k^e$\\\hline &size& N& M& M& K\\ \hline \hline \shortstack{$\phi^{i1}=\phi^{i2}$ \\ at $r=a_1$}&N& $(h-d_1)~\mathrm{diag}(\vec{R}_{1n}^{i1})$& $-\boldsymbol{\mathcal{Z}}_{nm}\odot 1_{N1}\vec{R}_{1m}^{i2}$& $-\boldsymbol{\mathcal{Z}}_{nm}\odot 1_{N1}\vec{R}_{2m}^{i2}$& $0_{NK}$\\ \hline @@ -346,9 +345,8 @@ \subsection{Block Matrix Structure} \shortstack{$\frac{\partial}{\partial r}\phi^{i1}=\frac{\partial}{\partial r}\phi^{i2}$ \\ at $r=a_1$}&M& $- \boldsymbol{\mathcal{Z}}_{mn} \odot 1_{M1} \frac{\partial}{\partial r}\vec{R}_{1n}^{i1}$& $(h-d_2)~\mathrm{diag}(\frac{\partial}{\partial r}\vec{R}_{1m}^{i2})$& $(h-d_2)~\mathrm{diag}(\frac{\partial}{\partial r}\vec{R}_{2m}^{i2})$& $0_{MK}$\\ \hline \shortstack{$\frac{\partial}{\partial r}\phi^{i2}=\frac{\partial}{\partial r}\phi^{e}$ \\ at $r=a_2$}&K& $0_{KN}$& $-\boldsymbol{\mathcal{Z}}_{km} \odot 1_{K1}\frac{\partial}{\partial r}\vec{R}_{1m}^{i2}$& $-\boldsymbol{\mathcal{Z}}_{km}\odot 1_{K1}\frac{\partial}{\partial r}\vec{R}_{2m}^{i2}$& $h~\mathrm{diag}(\frac{\partial}{\partial r}\vec{\Lambda}_k)$\\ \hline \end{tabular} - \fillandplacepagenumber \end{table} -\end{landscape} + \begin{table}[htbp] \centering \caption{MEEM b-vector} @@ -375,7 +373,7 @@ \subsection{Block Matrix Structure} An even sparser matrix could be obtained with the alternate eigenfunction scaling for the second region described in \citet{chau_inertia_2012}. \begin{figure}[htbp] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/meem_sparsity.pdf} + \includegraphics[width=0.45\linewidth]{figs/from-matlab/meem_sparsity.pdf} \caption{A-matrix sparsity pattern, shown for $N=M=K=4$} \label{fig:sparsity} \end{figure} @@ -441,7 +439,7 @@ \subsection{Calculation of Outputs} The first and second term are independent of and scale linearly with the eigencoefficients $\vec{x}$, respectively. The radiation coefficients are thus computed as -\begin{equation} +\begin{equation}\label{eq:hydro} A_h + \frac{iB_h}{\omega}=2\pi \rho h^3(c_0 + \vec{c}\cdot\vec{x}) =2\pi \rho h^3(c_0 + \vec{c} \cdot A^{-1} \vec{b}) \end{equation} with output constant $c_0$ and output row vector $\vec{c}$ defined in \Cref{eq:c0,tab:MEEM-c-vector} respectively: @@ -482,22 +480,23 @@ \subsection{Calculation of Outputs} \begin{equation}\label{eq:gamma-K} |\gamma| = \sqrt{\frac{ 4 \rho_w g V_g B_h} {m_0}}, \quad % excitation \angle \gamma = -\frac{\pi}{2} + \angle\frac{ B_{k=0}^e}{\textrm{H}_0^{1}(m_0 a_2)},\quad - K_h = \rho_w g A_w % hydrostatic stiffness + K_h = \rho_w g \underbrace{\frac{\pi}{4} (D_f^2 - D_s^2)}_{A_{w,f}},\quad % hydrostatic stiffness + K_s = \rho_w g \underbrace{\frac{\pi}{4} D_s^2}_{A_{w,s}} \end{equation} -where $g$ is the acceleration due to gravity, $\rho_w$ is the density of water, $m_0$ is the wavenumber, $V_g$ is the finite depth group velocity, $\textrm{H}_0^1$ is the zeroth-order Hankel function of the first kind, and $A_w= \frac{\pi}{4} (D_f^2 - D_s^2)$ is the waterplane area \citep{newman}. +where $g$ is the acceleration due to gravity, $\rho_w$ is the density of water, $m_0$ is the wavenumber, $V_g$ is the finite depth group velocity, $\textrm{H}_0^1$ is the zeroth-order Hankel function of the first kind, and $A_w$ is the waterplane area \citep{newman}. This method of calculating excitation from damping is the well-known Haskind relation. Note that while the excitation magnitude $|\gamma|$ depends on the radiation damping $B_h$ which in turn depends on all the inner region eigencoefficients ($\vec{C}_{m}^{i2}$ for float excitation and $\vec{C}_{1n}^{i1}$ for the spar excitation), the excitation phase $\angle\gamma$ depends only on the first exterior eigencoefficient, $B_{k=0}^e$. -\subsection{Validation} +\subsection{Verification and Validation} -Hydro coefficient results are validated by comparing to a benchmark shallow-water concentric-cylinder MEEM solution in \citet{chau_inertia_2012}. +Hydro coefficient results are verified by comparing to a benchmark shallow-water concentric-cylinder MEEM solution in \citet{chau_inertia_2012}. Excellent agreement is observed, shown in \Cref{fig:meem-yeung-validation}. The results are also experimentally validated in the studies \citet{chau_inertia_2012,son_performance_2016}. Previously in \Cref{sec:dynamic-validation}, the N=M=K=11 case was compared to WAMIT results for RM3 in deep water. \begin{figure}[htbp] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/meem_validation.pdf} + \includegraphics[width=0.5\linewidth]{figs/from-matlab/meem_validation.pdf} \caption{Nondimensional added mass and damping coefficient validation against \citep{chau_inertia_2012}} \label{fig:meem-yeung-validation} \end{figure} @@ -510,7 +509,7 @@ \subsection{Convergence} \Cref{fig:meem-matching} shows the matching behavior for $N=M=K=11$, where potential matches well but velocity still has noticeable mismatch. \begin{figure}[htbp] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/meem_matching.pdf} + \includegraphics[width=0.5\linewidth]{figs/from-matlab/meem_matching.pdf} \caption{Matching for $N=M=K=11$ for benchmark geometry} \label{fig:meem-matching} \end{figure} @@ -519,7 +518,7 @@ \subsection{Convergence} There, damping converges well at low frequencies but requires more harmonics at higher frequencies, while added mass has similar convergence across frequencies. \begin{figure}[htbp] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/meem_convergence_vs_omega.pdf} + \includegraphics[width=0.5\linewidth]{figs/from-matlab/meem_convergence_vs_omega.pdf} \caption{Convergence for $N=M=K=(5,10,20,30)$ for RM3} \label{fig:meem-convergence} \end{figure} @@ -548,7 +547,7 @@ \subsection{Numerical Notes}\label{sec:appendix-meem-numerics} \end{equation} \begin{figure}[htbp] \centering - \includegraphics[width=.75\linewidth]{figs/from-matlab/asymptotic_b_vector.pdf} + \includegraphics[width=.5\linewidth]{figs/from-matlab/asymptotic_b_vector.pdf} \caption{Asymptotic b-vector for large $m_0h$} \label{fig:meem-b-limit} \end{figure} diff --git a/pubs/applied-ocean-research-model/sections/model-intro.tex b/pubs/applied-ocean-research-model/sections/model-intro.tex index f09c5265..a60ff345 100644 --- a/pubs/applied-ocean-research-model/sections/model-intro.tex +++ b/pubs/applied-ocean-research-model/sections/model-intro.tex @@ -23,8 +23,8 @@ \subsection{Wave Energy Overview} Automating this process through optimization allows WEC designers to more quickly and systematically explore tradeoffs of major design decisions, and could eventually allow standardized comparisons between different WEC architectures. Thus, performing a system-level WEC optimization that simultaneously considers the coupled techno-economic goals, design decisions, and requirements is highly advantageous. \else - The strong interdisciplinary coupling of WEC design — between powertrain, controller, bulk dimensions, hydrodynamic shape, and structural thicknesses — complicates conventional design techniques. - Systematically exploring this coupled tradeoff requires a fast multidisciplinary simulation that can be embedded in an optimization loop. + The strong interdisciplinary coupling of between powertrain, controller, bulk dimensions, hydrodynamic shape, and structural thicknesses complicates conventional WEC design techniques. + Systematically exploring these coupled tradeoffs requires a fast multidisciplinary simulation that can be embedded in an optimization loop. \fi %\Cref{fig:disciplines} is a non-exhaustive outline of the most common disciplines in a wave energy converter, and the level and direction of coupling between them. @@ -62,21 +62,11 @@ \subsection{WEC System Modeling: State of the Art and Gaps} In industry, multidisciplinary modeling is common but typically sequential rather than concurrent; a survey of 25 WEC designers \citep{trueworthy_wave_2020} found that only 50\% design all subsystems concurrently, and identified concurrent design as an under-utilized technique worthy of further study. \fi -The academic state of the art typically focuses on modeling and optimizing one or a few disciplines at a time, with the most common being hydrodynamics, dynamics, and controls. +The academic state-of-the-art typically focuses on modeling and optimizing one or a few disciplines at a time, with the most common being hydrodynamics, dynamics, and controls. Major modeling choices include the representation (or lack thereof) of the WEC's hydrodynamics, dynamics, controls, mooring, powertrain, structures, and economic viability, as well as the sea states considered. % figure of taxonomy \Cref{fig:model-taxonomy} provides a simplified taxonomy of the various WEC system modeling approaches used in the literature, and \Cref{tab:lit} compares the most relevant studies along this taxonomy. -\begin{landscape} -\begingroup -\begin{figure}[htbp] -\centering -\includegraphics[width=1.1\linewidth]{figs/from-matlab/taxonomy.pdf} -\caption{WEC system modeling approaches, with MDOcean's capabilities outlined in black.}\label{fig:model-taxonomy} -\fillandplacepagenumber -\end{figure} -\endgroup -\end{landscape} % table summary @@ -104,93 +94,6 @@ \subsection{WEC System Modeling: State of the Art and Gaps} The dominant trends across these studies are: point absorbers as the most common architecture; boundary element method (BEM) hydrodynamics with occasional semi-analytical (MEEM) treatments; drag handled via linear damping or describing functions; dynamics in either frequency or time domain (with extensions of either to handle constraints); and limited inclusion of structures, mooring, and cost modeling. \fi -%Need to copy sources over from slide 27 and 35 of F24 reseach slides with original charts (with color) -\ifdefined\DISSERTATION - \clearpage -\else - \onecolumn -\fi -\begin{longtable}{ - >{\centering\arraybackslash}p{0.045\linewidth}| - >{\centering\arraybackslash}p{0.081\linewidth}| - >{\centering\arraybackslash}p{0.08\linewidth}| % hydro - >{\centering\arraybackslash}p{0.05\linewidth}| % drag - >{\centering\arraybackslash}p{0.035\linewidth}| % dynam - >{\centering\arraybackslash}p{0.08\linewidth}| % controls - >{\centering\arraybackslash}p{0.04\linewidth}| % mooring - >{\centering\arraybackslash}p{0.035\linewidth}| - >{\centering\arraybackslash}p{0.059\linewidth}| - >{\centering\arraybackslash}p{0.065\linewidth}| % econ - >{\centering\arraybackslash}p{0.072\linewidth}} -\rot{\textbf{Ref}} & \rot{\textbf{Type}} & \rot{\textbf{Hydro}}& \rot{\textbf{Drag}} & \rot{\textbf{Dyn}} & \rot{\textbf{Ctrl}} & \rot{\textbf{Moor}} & \rot{\textbf{PTO}} & \rot{\textbf{Struct}} & \rot{\textbf{Econ}} & \rot{\textbf{Waves}} \\ -\hline -This work & 2-PA & MEEM & DF & F+ & PI+ & - & DY & AN & STR, PTO & REG+, STO \\ - -\cite{mccabe_multidisciplinary_2022} & 2-PA & ALG & - & F+ & PI+ & - & EF & AN & STR & REG+, STO \\ - -\cite{khanal_multi-objective_2024} & 1-PA* & BEM & - & F & PI+ & - & - & - & GEO & REG \\ - -\cite{gaudin_single_2021} & 1-PA* & MEEM & - & F & PI & DE & - & AN & MOOR & IRR+, STO \\ - -\cite{edwards_optimisation_2022} & 1-PA & BEM & - & F & P & - & - & - & GEO & REG \\ - -\cite{garcia-teruel_reliability-based_2021}& 1-PA & BEM & LD & F+ & PI & - & - & LD & - & IRR+ \\ - -\cite{garcia-teruel_design_2022}& 1-PA & BEM & - & F+ & PI & - & - & - & - & IRR \\ - -\cite{cotten_multi-objective_2022} & ATN & BEM & - & F+ & P & - & - & LD & - & IRR+ \\ - -\cite{abdulkadir_control_2024} & 1-PA & MEEM & LD & T & NL P & - & - & - & - & IRR \\ - -\cite{housner_numerical_2024} & TRM & BEM & - & T & PI & DY & - & - & - & REG \\ - -\cite{al_shami_parameter_2019} & 2-PA & BEM & DF & F & PI & - & - & - & - & REG \\ - -\cite{RM3} & 2-PA & BEM & Q & T & P & DE & - & FEA & FULL & IRR+, STO \\ - -\cite{mi_multi-scale_2025} & OSW & BEM & Various & T & PI & DE & DY & FEA & FULL & IRR+, STO \\ - -\cite{an_optimal_2024} & OTP & CFD & CFD & T & - & - & - & FEA & - & STO \\ - -\cite{ambuhl_reliability-based_2014} & 1-PA & - & Q & - & - & DE & - & AN & STR & STO \\ - -\cite{nguyen_theoretical_2024} & OSW & MEEM & LD & F & PI & - & - & LD & - & IRR \\ - -\cite{ferri_balancing_2014} & 1-PA & BEM & LD & T & P,PI+, US & - & DY & LD & STR & IRR+ \\ - -\cite{rosati_control_2023} & OWC & BEM & - & T & NL P+ & - & EF & - & PTO & IRR+ \\ - -\cite{son_performance_2016} & 2-PA & MEEM & LD & F & P & - & EF & - & - & REG \\ - -\cite{gaebele_tpl_2025} & 2-PA & BEM & - & PS & PI & - & DY & - & PTO, GEO & IRR \\ - -\cite{devin_high-dimensional_2024} & 1-PA & FIT & - & PS & US & - & DY & - & - & REG \\ - -\cite{grasberger_control_2024} & OSW & BEM & - & PS & PI & - & DY & - & GEO & IRR \\ - -\cite{herber_dynamic_2014} & 1-PA & FIT & LD & T & US & - & - & - & - & IRR \\ - -\cite{lin_fast_2025} & 1-PA & BEM & LD & F+ & US & - & - & - & - & IRR+ \\ - -\cite{abdulkadir_optimal_2024} & 1-PA* & BEM & - & T & US & - & - & - & - & IRR \\ - -%\cite{zhang_performance_2024} & 3-PA & MEEM & LD & 3 & F & P & DY & - & - & - & IRR -%& BF & PWR & - & CTL & 2 \\ - -%\cite{wen_shape_2018} & 1-PA & BEM & - & 1 & F & P & - & - & - & - & IRR+ -%& DOE, LOC & PWR & STA & DIM, SHP & 3 \\ - -\hline -\caption{Comparison of the disciplinary scope and fidelity of optimization-relevant WEC models. -See \Cref{fig:model-taxonomy} for abbreviations. * in type column means array effects from multiple interacting devices are modeled.} -\label{tab:lit} -%\fillandplacepagenumber -\end{longtable} -\ifdefined\DISSERTATION -\else - \twocolumn -\fi - % Per-study paragraphs \ifdefined\DISSERTATION \paragraph{Notable large-scale optimization studies} @@ -263,6 +166,96 @@ \subsection{WEC System Modeling: State of the Art and Gaps} Existing studies therefore tend to sacrifice either disciplinary breadth or computational tractability: high-fidelity multidisciplinary models are generally too expensive for large-scale optimization, while optimization-oriented models often omit structures, powertrain dynamics, or economic assessment. This tradeoff motivates the development of a computationally efficient framework capable of concurrently capturing the major techno-economic couplings governing WEC performance. \fi +\begin{landscape} +\begingroup +\begin{figure}[htbp] +\centering +\includegraphics[width=1.05\linewidth]{figs/from-matlab/taxonomy.pdf} +\caption{WEC system modeling approaches, with MDOcean's capabilities outlined in black.}\label{fig:model-taxonomy} +\fillandplacepagenumber +\end{figure} +\endgroup +\end{landscape} +%Need to copy sources over from slide 27 and 35 of F24 reseach slides with original charts (with color) +\singleColMacro{ +\begin{longtable}{ + P{0.13\linewidth}| % Ref + P{0.055\linewidth}| % Type + P{0.08\linewidth}| % hydro + P{0.04\linewidth}| % drag + P{0.025\linewidth}| % dynam + P{0.05\linewidth}| % controls + P{0.04\linewidth}| % mooring + P{0.035\linewidth}| % PTO + P{0.05\linewidth}| % struct + P{0.06\linewidth}| % econ + P{0.065\linewidth}} % waves +\rot{\textbf{Ref}} & \rot{\textbf{Type}} & \rot{\textbf{Hydro}}& \rot{\textbf{Drag}} & \rot{\textbf{Dyn}} & \rot{\textbf{Ctrl}} & \rot{\textbf{Moor}} & \rot{\textbf{PTO}} & \rot{\textbf{Struct}} & \rot{\textbf{Econ}} & \rot{\textbf{Waves}} \\ +\hline +This work & 2-PA & MEEM, ALG, FIT & DF, M & F+ & PI+ & - & DY & AN & STR, PTO & REG+, STO \\ + +\cite{mccabe_multidisciplinary_2022} & 2-PA & ALG & - & F+ & PI+ & - & EF & AN & STR & REG+, STO \\ + +\cite{khanal_multi-objective_2024} & 1-PA* & BEM & - & F & PI+ & - & - & - & GEO & REG \\ + +\cite{gaudin_single_2021} & 1-PA* & MEEM & - & F & PI & DE & - & AN & MOOR & IRR+, STO \\ + +\cite{edwards_optimisation_2022} & 1-PA & BEM & - & F & P & - & - & - & GEO & REG \\ + +\cite{garcia-teruel_reliability-based_2021}& 1-PA & BEM & LD & F+ & PI & - & - & LD & - & IRR+ \\ + +\cite{garcia-teruel_design_2022}& 1-PA & BEM & - & F+ & PI & - & - & - & - & IRR \\ + +\cite{cotten_multi-objective_2022} & ATN & BEM & - & F+ & P & - & - & LD & - & IRR+ \\ + +\cite{abdulkadir_control_2024} & 1-PA & MEEM & LD & T & NL P & - & - & - & - & IRR \\ + +\cite{housner_numerical_2024} & TRM & BEM & - & T & PI & DY & - & - & - & REG \\ + +\cite{al_shami_parameter_2019} & 2-PA & BEM & DF & F & PI & - & - & - & - & REG \\ + +\cite{RM3} & 2-PA & BEM & Q & T & P & DE & - & FEA & FULL & IRR+, STO \\ + +\cite{mi_multi-scale_2025} & OSW & BEM & \small{Various} & T & PI & DE & DY & FEA & FULL & IRR+, STO \\ + +\cite{an_optimal_2024} & OTP & CFD & CFD & T & - & - & - & FEA & - & STO \\ + +\cite{ambuhl_reliability-based_2014} & 1-PA & - & Q & - & - & DE & - & AN & STR & STO \\ + +\cite{nguyen_theoretical_2024} & OSW & MEEM & LD & F & PI & - & - & LD & - & IRR \\ + +\cite{ferri_balancing_2014} & 1-PA & BEM & LD & T & P,PI+, US & - & DY & LD & STR & IRR+ \\ + +\cite{rosati_control_2023} & OWC & BEM & - & T & NL P+ & - & EF & - & PTO & IRR+ \\ + +\cite{son_performance_2016} & 2-PA & MEEM & LD & F & P & - & EF & - & - & REG \\ + +\cite{gaebele_tpl_2025} & 2-PA & BEM & - & PS & PI & - & DY & - & PTO, GEO & IRR \\ + +\cite{devin_high-dimensional_2024} & 1-PA & FIT & - & PS & US & - & DY & - & - & REG \\ + +\cite{grasberger_control_2024} & OSW & BEM & - & PS & PI & - & DY & - & GEO & IRR \\ + +\cite{herber_dynamic_2014} & 1-PA & FIT & LD & T & US & - & - & - & - & IRR \\ + +\cite{lin_fast_2025} & 1-PA & BEM & LD & F+ & US & - & - & - & - & IRR+ \\ + +\cite{abdulkadir_optimal_2024} & 1-PA* & BEM & - & T & US & - & - & - & - & IRR \\ + +%\cite{zhang_performance_2024} & 3-PA & MEEM & LD & 3 & F & P & DY & - & - & - & IRR +%& BF & PWR & - & CTL & 2 \\ + +%\cite{wen_shape_2018} & 1-PA & BEM & - & 1 & F & P & - & - & - & - & IRR+ +%& DOE, LOC & PWR & STA & DIM, SHP & 3 \\ + +\hline +\caption{Comparison of the disciplinary scope and fidelity of optimization-relevant WEC models. +See \Cref{fig:model-taxonomy} for abbreviations. +An asterisk * in type column means array effects from multiple interacting devices are modeled.} +\label{tab:lit} +%\fillandplacepagenumber +\end{longtable} +} \subsection{Relevance of Semi-Analytical Modeling} An optimization study can use any kind of simulation model to assess the design's performance, including simple algebraic models, theory-intensive semi-analytical models, and standard numerical models. @@ -272,8 +265,8 @@ \subsection{Relevance of Semi-Analytical Modeling} \centering \caption{Comparison of model types} \label{tab:model-types} -\begin{tabular}{c>{\centering\arraybackslash}p{0.25\linewidth}c>{\centering\arraybackslash}p{0.2\linewidth}} - \textbf{Model type}& \textbf{Expertise required}& \textbf{Accuracy}& \textbf{Computational cost}\\ \hline +\begin{tabular}{cM{0.14\linewidth}cM{0.13\linewidth}} + \textbf{Model type}& \textbf{Expertise required}& \textbf{Accuracy}& \textbf{Computation cost}\\ \hline \textbf{Algebraic}& Low& Low& Low\\ \textbf{Semi-analytical}& High& Med& Low\\ \textbf{Numerical}& Med& Med/High& High\\ @@ -311,16 +304,16 @@ \subsection{Paper Contributions and Roadmap} The modeling framework articulated here is scalable to more disciplines, and enables systematic WEC design optimization. \else This paper presents MDOcean \citep{mccabe_mdocean_2024}, an open-source semi-analytical WEC simulation framework addressing the gap identified above. - The contributions are: - \begin{itemize} - \item integration of semi-analytical models from hydrodynamics, dynamics, control, structures, and economics in a single framework suitable for optimization; - \item a linearized pseudo-spectral dynamics formulation that extends prior describing-function and constrained-control approaches with an analytical quadratically-constrained quadratic program (QCQP) for the constrained linear controller, with a geometric interpretation on the complex reflection-coefficient plane; - \item validation against WEC-Sim and the RM3 reference design across power, structural, and economic outputs; and - \item demonstration of computational performance (\resultsAOR[simRuntime] simulation runtime) enabling optimization studies that are otherwise infeasible. + The work provides the following contributions: + \begin{itemize}[leftmargin=*] + \item \textbf{Integrates semi-analytical models} from hydrodynamics, dynamics, control, structures, and economics in a single framework suitable for optimization + \item Formulates an underactuated optimal control problem with coupled nonlinear multi-port dynamics, and introduces the \textbf{linearized pseudo-spectral (LPS) method} to approximate time-domain constraints and nonlinearities + \item \textbf{Analytically solves} the LPS optimal control problem as a univariate quadratically-constrained quadratic program (QCQP) with a geometric interpretation on the complex reflection coefficient plane + \item \textbf{Validates} against WEC-Sim and the RM3 reference design across power, structural, and economic outputs + \item Demonstrates \textbf{computational performance} (\resultsAOR[simRuntime] simulation runtime) enabling optimization studies that are otherwise infeasible + \item Illustrates the \textbf{mathematical structure and effect magnitude} of subsystem interactions and techno-economic scaling relationships, providing intuition to facilitate early-stage design tradeoffs \end{itemize} - Beyond enabling optimization, the framework is intended to support physical interpretation of subsystem interactions, limiting behaviors, and techno-economic scaling relationships relevant to early-stage WEC design. \fi - \ifdefined\DISSERTATION In \Cref{sec:model-structure}, the research objective and design scope are first used to define a broad problem formulation and module decomposition. Next, \Cref{sec:modules} describes the construction of the simulation model. @@ -329,7 +322,7 @@ \subsection{Paper Contributions and Roadmap} Finally, \Cref{sec:discussion} showcases results from model sweeps and insights from the model structure. \else \Cref{sec:model-structure} introduces MDOcean's modular structure and the rationale behind its analysis architecture. - \Cref{sec:modules} develops each module -- geometry, hydrodynamics, dynamics and control, structures, and economics -- documenting the assumptions, analysis methods, and implementation choices. + \Cref{sec:modules} documents the assumptions, analysis methods, and implementation choices of each module -- geometry, hydrodynamics, dynamics and control, structures, and economics. \Cref{sec:validation-benchmarking} compares MDOcean against WEC-Sim and the RM3 reference design, demonstrates the runtime gain over established tools, and discusses where the model is reliable. \Cref{sec:discussion} presents results from variable sweeps and model-derived insights, along with limitations and future work. Extended derivations, validation studies, and implementation details are provided in the appendices to preserve readability of the primary narrative while maintaining reproducibility. diff --git a/pubs/applied-ocean-research-model/sections/model-overview.tex b/pubs/applied-ocean-research-model/sections/model-overview.tex index dd071973..6150e71a 100644 --- a/pubs/applied-ocean-research-model/sections/model-overview.tex +++ b/pubs/applied-ocean-research-model/sections/model-overview.tex @@ -13,11 +13,18 @@ \section{Model Structure} The module coupling structure is illustrated in the XDSM diagram in \Cref{fig:n2}. Inputs and outputs are organized according to standard MDO conventions \citep{lambe_extensions_2012}, and an optimizer would update design variables $x$ using objective $J$ and constraints $g$ until convergence to $x^*$. -\begin{figure}[htbp] + +\begin{figure*}[htbp] \centering \includegraphics[width=\linewidth]{figs/from-matlab/xdsm.pdf} \caption{Simplified XDSM diagram}\label{fig:n2} -\end{figure} +\end{figure*} +\begin{figure*}[b!] +\centering +\includegraphics[width=1\linewidth]{figs/from-matlab/control_analysis_flowcharts.pdf} +\caption{Optimization architectures organized by analysis method (MDF and SAND) and control method (explicit, implicit, simultaneous, and nested).} +\label{fig:control-arch} +\end{figure*} \ifdefined\DISSERTATION \paragraph{Absence of Feedback Coupling} @@ -102,17 +109,12 @@ \section{Model Structure} MDOcean uses an implicit-MDF architecture: the dynamics module iterates internally to converge the quasi-linear drag and saturation describing functions, while all other modules are explicit. The implicit formulation is preferred over explicit alternatives due to the structure of the derived optimal control conditions (\Cref{sec:mod-freq-domain}). It also avoids introducing control variables into the outer optimization while preserving optimal control behavior, minimizing total runtime given that hydrodynamics dominates computational cost (\Cref{sec:sim-runtime}). - Detailed architecture rationale, coupling analysis, and the relation to MDO literature are provided in [dissertation chapter]. + Detailed architecture rationale, coupling analysis, and the relation to MDO literature are provided in \cite{mccabe_dissertation_2026}. \Cref{fig:control-arch} summarizes the relationship between analysis and control architectures, with the implicit-MDF formulation used in this study highlighted. \fi -\begin{figure}[b!] -\centering -\includegraphics[width=1\linewidth]{figs/from-matlab/control_analysis_flowcharts.pdf} -\caption{Optimization architectures organized by analysis method (MDF and SAND) and control method (explicit, implicit, simultaneous, and nested).} -\label{fig:control-arch} -\end{figure} + diff --git a/pubs/applied-ocean-research-model/sections/module-details.tex b/pubs/applied-ocean-research-model/sections/module-details.tex index 47d89ea9..f2dbe7de 100644 --- a/pubs/applied-ocean-research-model/sections/module-details.tex +++ b/pubs/applied-ocean-research-model/sections/module-details.tex @@ -22,32 +22,32 @@ \subsection{Geometry}\label{sec:geom} The geometry module computes submerged volume $V_{\text{sub}}$, structural volume $V_{\text{struct}}$, and associated hydrostatic properties from bulk dimensions. \Cref{fig:dims} defines the principal geometric variables, where $T$ denotes drafts, $h$ heights, and $D$ diameters. Subscripts $f$, $s$, and $d$ refer to the float, spar, and damping plate, respectively. -Key points include the center of buoyancy $B$, center of gravity $G$, and metacenter $M$ for the combined system. +Key points include the center of buoyancy $B$, center of gravity $G$, the keel $K$, and metacenter $M$ for the combined system. Additional structural dimensions and thickness definitions are provided in \Cref{sec:structures}. \begin{figure}[htbp] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/dimensions.pdf} + \includegraphics[width=\linewidth]{figs/from-matlab/dimensions.pdf} \caption{Dimension labeling of system}\label{fig:dims} \end{figure} -Static stability is enforced by requiring $GM>0$, where +Static pitch stability is enforced by requiring: \begin{equation}\label{eq:GM} - GM = KB + BM - KG, + \overline{GM} = \overline{KB} + \overline{BM} - \overline{KG} > 0, \end{equation} -\noindent where $BM = \pi D_f^4/(64V_{\text{sub}})$ \cite{newman}. +\noindent where $\overline{BM} = \pi D_f^4/(64V_{\text{sub}})$ \citep{newman}. -Hydrostatic equilibrium is enforced by matching displaced water mass to total system mass and any mismatch is compensated using ballast. +Hydrostatic equilibrium is enforced by matching displaced water mass to total system mass, with any mismatch compensated using ballast with mass $m_{bal}$: \begin{equation} m_{bal} = \rho_w V_{sub} - \rho_M V_{struct}, \end{equation} where $\rho_w$ and $\rho_M$ are water and structural material densities. To ensure feasible buoyancy, $m_{bal}$ must be non-negative and physically storable within the hull. -Assuming seawater ballast, this imposes a volume feasibility constraint: +Assuming seawater ballast, this imposes a volume constraint: \begin{equation}\label{eq:vol-constraint} V_{struct} \le V_{sub} + V_{surf} - V_{pto}, \end{equation} -where $V_{\text{surf}}$ is the above-water structural volume and $V_{pto}$ is the (fixed) PTO volume. +where $V_{\text{surf}}$ (design-dependent) and $V_{pto}$ (fixed) are volumes of the above-water structures and PTO respectively. %TC:break Hydrodynamic Coefficients %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -59,52 +59,70 @@ \subsection{Hydrodynamic Coefficients}\label{sec:hydro} where, under monochromatic forcing, \begin{equation}\label{eq:hydro-forces-expanded} \begin{aligned} - \vec{\hat{F}}_{e} &= \vec{\gamma}\,\hat{\zeta}, - & \vec{\hat{F}}_{rad} &= -\mathbf{A}_h\vec{\hat{\ddot{\xi}}} - \mathbf{B}_h\vec{\hat{\dot{\xi}}}, \\ + \vec{\hat{F}}_{e} &= \vec{\gamma}\,\hat{\zeta},&~~ + \vec{\hat{F}}_{rad} &= -\mathbf{A}_h\vec{\hat{\ddot{\xi}}} - \mathbf{B}_h\vec{\hat{\dot{\xi}}},&~~ \vec{\hat{F}}_{res} &= -\mathbf{K}_h\vec{\hat{\xi}}. \end{aligned} \end{equation} -Here $\vec{\hat{\xi}}$ is the body displacement phasor vector, and the hydrodynamic coefficients $\vec{\gamma}$, $\mathbf{A}_h$, $\mathbf{B}_h$, and $\mathbf{K}_h$ are the wave excitation vector, added mass matrix, radiation damping matrix, and restoring stiffness matrix, respectively. -The origin is taken as the center of the float at the still waterline. +Here $\vec{\hat{\xi}}$ is the body displacement phasor vector from the origin at the center of the float at the still waterline. +The hydrodynamic coefficients $\vec{\gamma}$, $\mathbf{A}_h$, $\mathbf{B}_h$, and $\mathbf{K}_h$ are the wave excitation vector, added mass matrix, radiation damping matrix, and hydrostatic stiffness matrix, respectively. -For a two-body heaving point absorber under operational conditions ($op$), the hydrodynamic matrices are $2\times 2$, representing the float ($f$), spar ($s$), and their coupling ($c$): -\begin{equation}\label{eq:op-hydro-coeffs} +The displacement vector $\vec{\hat{\xi}}$ contains a state for each degree of freedom. +We define separate displacement and excitation vectors for the operational ($op$) and storm ($st$) cases. +A two-body heaving point absorber that rigidly locks the bodies in a storm following \cite{RM3} has vectors: +\begin{equation}\label{eq:define-vectors} \vec{\hat{\xi}}_{op} = \begin{bmatrix} \hat{\xi}_f \\ \hat{\xi}_s - \end{bmatrix},~ - \vec{\gamma}_{op} = \begin{bmatrix} + \end{bmatrix}, ~~ + \vec{\hat{\xi}}_{st}=\begin{bmatrix}\hat{\xi}_m + \end{bmatrix}, ~~ + \vec{\gamma}_{op} = \begin{bmatrix} \gamma_f \\ \gamma_s - \end{bmatrix},~ - \mathbf{A_{h,op}} = \begin{bmatrix} + \end{bmatrix},~~ + \vec{\gamma}_{st} = \begin{bmatrix} + \gamma_m + \end{bmatrix}. +\end{equation} +with subscripts for the float ($f$), spar ($s$), their coupling ($c$), and the merged float-spar ($m$). +Operational hydrodynamic matrices are symmetric $2\times 2$: +\begin{equation}\label{eq:op-hydro-coeffs} +\begin{aligned} + \mathbf{A}_{h,op} = \begin{bmatrix} A_f & A_c \\ A_c & A_s \end{bmatrix},~ - \mathbf{B_{h,op}} = \begin{bmatrix} + \mathbf{B}_{h,op} &= \begin{bmatrix} B_f & B_c \\ B_c & B_s - \end{bmatrix},~ - \mathbf{K_{h,op}} = \begin{bmatrix} + \end{bmatrix},\\ + \mathbf{K}_{h,op} &= \begin{bmatrix} K_f & 0\\ 0 & K_s \end{bmatrix} +\end{aligned} \end{equation} - -In storm conditions ($st$), the float and spar are assumed rigidly connected following \cite{RM3}, yielding a single-body representation ($m$): +while storm (merged-body) coefficients are scalar: \begin{equation}\label{eq:st-hydro-coeffs} -\vec{\hat{\xi}}_{st}, \vec{\gamma}_{st}, \mathbf{A}_{h,st}, \mathbf{B}_{h,st}, \mathbf{K}_{h,st}. +\begin{aligned} +\vec{\gamma}_{st}&=\gamma_{m}=\Sigma~\vec{\gamma}_{op},& +~\mathbf{A}_{h,st}&=A_m=\Sigma~\mathbf{A}_{h,op}, \\ +~\mathbf{B}_{h,st}&=B_m=\Sigma~\mathbf{B}_{h,op},& +~\mathbf{K}_{h,st}&=K_m=\Sigma~\mathbf{K}_{h,op}. +\end{aligned} \end{equation} - -The merged-body coefficients are obtained by summing the individual and coupling contributions, e.g., $A_m = A_f + A_s + 2A_c$, with analogous expressions for $B_m$, while hydrostatic stiffness is additive. +with $\Sigma$ denoting the sum of all elements in the matrix or vector (``grand sum''). The remaining task is to compute the float, spar, and coupling coefficients, with operational and storm implications discussed in \Cref{sec:design-load-cases}. -The hydrostatic stiffness terms ($K_f$, $K_s$) are frequency-independent and are computed directly from geometry (see \cref{sec:appendix-meem-details}), while the remaining coefficients require solving the frequency-domain radiation boundary value problem for the velocity potential. +The hydrostatic stiffness terms ($K_f$, $K_s$) are frequency-independent and are computed directly from geometry (see \Cref{eq:gamma-K} in \Cref{sec:appendix-meem-details}), while the remaining coefficients require solving the frequency-domain radiation boundary value problem for the velocity potential. Traditionally, this is performed using a BEM solver, which discretizes the body surfaces and solves a large linear system. In MDOcean, float coefficients are instead computed using the semi-analytical Matched Eigenfunction Expansion Method (MEEM), exploiting cylindrical symmetry for improved efficiency. Spar and coupling terms are obtained via algebraic approximations and interpolation of pre-computed BEM data, and could be extended to MEEM in future work. +\Cref{tab:hydro-methods} summarizes the methods used for each coefficient. %%%%%%%%%%%%%%% \subsubsection{Float Coefficients: Matched Eigenfunction Expansion Method}\label{sec:hydro-meem} MEEM solves the radiation and excitation problems by expanding the velocity potential in eigenfunctions within each fluid region and enforcing matching across region boundaries. -The MEEM radiation solution for two concentric heaving cylinders was first presented in the studies \cite{mavrakos_hydrodynamic_2004} and detailed in reference \cite{chau_inertia_2012}. +The MEEM radiation solution for two concentric heaving cylinders was first presented by \citet{mavrakos_hydrodynamic_2004} and detailed by \citet{chau_inertia_2012}. In this approach, the fluid domain is partitioned into regions with separable analytical solutions, and continuity conditions are enforced at region interfaces to determine the unknown expansion coefficients. -The open-source implementation used here, summarized in \cref{sec:appendix-meem-details}, follows the formulation derived in the study \cite{mccabe_open-source_2024}, with numerical properties detailed in the forthcoming companion paper \cite{bimali_matrix_2026,best_openflash_2026}. +The open-source implementation used here, summarized in \cref{sec:appendix-meem-details}, follows the formulation by \citet{mccabe_open-source_2024}. +Forthcoming companion papers \citet{bimali_matrix_2026,best_openflash_2026} detail the numerical properties and code implementation respectively of an expanded implementation. The primary advantage of MEEM is computational efficiency. Timing comparisons (\Cref{sec:sim-runtime}) show more than an order-of-magnitude speedup over Capytaine BEM at comparable convergence, significantly reducing the cost of multidisciplinary analysis. @@ -114,50 +132,44 @@ \subsubsection{Float Coefficients: Matched Eigenfunction Expansion Method}\label As a semi-analytical method, MEEM is restricted to simplified geometries; the formulation used here assumes a dual concentric cylinder (\Cref{fig:meem-geom}). Consequently, float coefficients neglect the damping plate and approximate the float bottom as flat rather than slanted. The damping plate enters through the spar coefficients (next subsection). -These approximations have minor impact but could be relaxed in future extensions using more general MEEM formulations \cite{olaya_hydrodynamic_2015,kokkinowrachos_behaviour_1986}. +These approximations have minor impact and could be relaxed in future extensions using more general MEEM formulations \citep{olaya_hydrodynamic_2015,kokkinowrachos_behaviour_1986,bimali_matrix_2026}. \begin{figure}[b!] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/MEEM_geometry.pdf} + \includegraphics[width=\linewidth]{figs/from-matlab/MEEM_geometry.pdf} \caption{MEEM geometry. Dashed lines represent geometry that is neglected in the MEEM calculations} \label{fig:meem-geom} \end{figure} %%%%%%%%%%%%%%% -\subsubsection{Other Coefficients: Approximations}\label{sec:hydro-other} +\subsubsection{Other coefficients: approximations}\label{sec:hydro-other} Because the damping plate affects the spar coefficients more strongly than the float, the former are obtained by scaling existing solutions that include the damping plate rather than using MEEM. -\Cref{tab:hydro-methods} summarizes the methods used for each coefficient. \begin{table}[htbp] \centering - \caption{Method of computing different hydrodynamic coefficients in MDOcean} + \caption{Method of computing hydrodynamic coefficients in MDOcean} \label{tab:hydro-methods} -\begin{tabular}{cc} + \begin{tabular} + {M{0.25\linewidth}M{0.63\linewidth}} Coefficient& Method\\ \hline - $A_f$, $B_f$, $\gamma_f$ & MEEM\\ - $A_s$& Approximate $\omega\rightarrow\infty$ solution (\cref{eq:As})\\ - $A_c$, $B_s$, $B_c$, $\gamma_s$& Nominal BEM solution, scaled with $D_d$ and $T_s$\\ + $A_f$, $B_f$, $\gamma_f$ & MEEM (\Cref{eq:hydro,eq:gamma-K} in \Cref{sec:appendix-meem-details})\\ + $A_s$& Approximate $\omega\rightarrow\infty$ solution (\Cref{eq:As} in \Cref{sec:appendix-additional-hydro})\\ + $A_c$, $B_s$, $B_c$, $\gamma_s$& Nominal BEM solution, scaled with $D_d$ and $T_s$ (\Cref{eq:scale-wamit} in \Cref{sec:appendix-additional-hydro})\\ \end{tabular} - \end{table} - +\end{table} %\hl{Describe how I did the WAMIT damping plate scaling for the spar excitation.} - %The MEEM implementation currently cannot take the damping plate into account, and the plate is expected to have a large impact on the spar added mass and damping. Therefore, algebraic approximations of the damping plate coefficients are used instead of the MEEM solution. -The spar added mass coefficient $A_s$ is approximated as a frequency-independent quantity based on the displaced water volume in the spherical projection of the damping plate \cite{philip_damping_2012}: -\begin{equation}\label{eq:As} - A_s =\frac{1}{3}\rho_w D_d^3\left(1 - \frac{3}{4}r^2\sqrt{1-r^2} - \frac{1}{4}(1-\sqrt{1-r^2})^2(2+\sqrt{1-r^2})\right), -\end{equation} -where $r = D_s/D_d$. -This expression smoothly transitions from a pure plate limit ($\frac{1}{3}\rho_w D_d^3$ at $r=0$) to a pure column limit ($\frac{1}{6}\rho_w D_d^3$ at $r=1$). +The spar added mass coefficient $A_s$ is approximated as a frequency-independent quantity based on the displaced water volume in the spherical projection of the damping plate, with the formula given in \Cref{sec:appendix-additional-hydro}. The spar damping ($B_s$), excitation ($\gamma_s$), and coupling terms ($A_c$, $B_c$) are obtained by interpolating a pre-computed WAMIT BEM solution over the nondimensional parameter $kD_d$ (see \cref{sec:appendix-additional-hydro}). -The excitation term is additionally scaled as $e^{-k(T_s - T_{s,nom})}$ to account for depth-dependent attenuation. +The excitation term is additionally scaled as $\exp(-k(T_s - T_{s,nom}))$ to account for depth-dependent attenuation. %\hl{I suspect this imposes a constraint on $D_d$ needing to be similar to the nominal design, see notebook 3/20/26 referencing analysis from 4/4/25.} Because float and spar coefficients are computed using different approximations, positive definiteness of the hydrodynamic matrices is not guaranteed (see \cref{sec:appendix-additional-hydro}). -To ensure physical consistency, the added mass and damping matrices are required to remain positive definite, and any violations are corrected by reducing the coupling terms $A_c$ and $B_c$. +To ensure physicality, the added mass and damping matrices are required to remain positive definite. +Any violations are corrected by reducing the coupling terms $A_c$ and $B_c$. % \subsubsection{Spar Dynamics} % The RM3 is a two-body system where the relative motion of the float and spar produces power. However, capturing these underactuated multibody dynamics would require more complicated modeling and control \cite{faedo_principle_2022}, especially in the context of force saturation%, as well as implementation of the excitation force phase in MEEM. @@ -198,42 +210,45 @@ \subsection{Dynamics and Control}\label{sec:dynamics} This subsection introduces the dynamics and control modeling assumptions before detailing the frequency-domain methods and control formulations. %%%%%%%%%%%%%%% -\subsubsection{linearized pseudo-spectral Techniques} +\subsubsection{Modified frequency-domain methods} \label{sec:mod-freq-domain} \ifdefined\DISSERTATION - Methods for WEC dynamics and control simulation in an optimization context were reviewed in \cref{sec:lit} and include time-stepping (nonlinear, numerical, high computation cost), frequency domain (linear, analytical, low computation cost), pseudo-spectral (nonlinear, numerical, medium computation cost), linearized pseudo-spectral techniques (approximate nonlinear, semi-analytical, low computational cost), and Pontryagin's maximum principle (nonlinear, analytical, low computational cost). + Methods for WEC dynamics and control simulation in an optimization context were reviewed in \cref{sec:lit} and include time-stepping (nonlinear, numerical, high computation cost), frequency domain (linear, analytical, low computation cost), pseudo-spectral (nonlinear, numerical, medium computation cost), quasi-linearized pseudo-spectral techniques (approximate nonlinear, semi-analytical, low computational cost), and Pontryagin's maximum principle (nonlinear, analytical, low computational cost). In wave energy, both spectral-domain numerical methods and analytical methods have generated recent research interest, with studies \cite{bacelli_numerical_2015,coe_initial_2020,tan_computationally_2026,nie_optimal_2016,da_silva_stochastic_2023} and studies \cite{abdulkadir_optimal_2024,lin_fast_2025} serving as illustrative examples of each case. - This section describes both the technical concept and the practical rationale for the choice of a linearized pseudo-spectral method, and contextualizes it with respect to other constrained optimal control methods in more detail than was possible in the literature review. + This section describes both the technical concept and the practical rationale for the choice of a quasi-linearized pseudo-spectral method, and contextualizes it with respect to other constrained optimal control methods in more detail than was possible in the literature review. MDOcean requires the ability to incorporate nonlinearities, at least approximately. Inherent nonlinearities such as drag significantly affect the amplitude and power production at resonance. Even with linearized intrinsic dynamics, maintaining optimality subject to dynamic constraints also requires a nonlinear controller. Dynamic constraints on amplitude, force, and peak power significantly influence the trade-off of power and cost, and amplitude limits are required to obtain a physically meaningful result. Besides nonlinearities, MDOcean also requires a low computation cost, since the dynamics will be evaluated many thousands of times during optimization. - Of the methods mentioned, Pontryagin's maximum principle (PMP) and linearized pseudo-spectral techniques meet these requirements. + Of the methods mentioned, Pontryagin's maximum principle (PMP) and quasi-linearized pseudo-spectral techniques meet these requirements. PMP requires an analytically difficult derivation unique to the specific nonlinearities in question because wave energy dynamics have a particular form requiring so-called singular control \cite{zou_optimal_2017}. - On the other hand, linearized pseudo-spectral techniques are simple and intuitive. - Therefore, MDOcean pursues a quasi-linear semi-analytical linearized pseudo-spectral method. + On the other hand, quasi-linearized pseudo-spectral techniques are simple and intuitive. + Therefore, MDOcean pursues a quasi-linear semi-analytical quasi-linearized pseudo-spectral method. \else - Standard linear frequency-domain WEC analysis assumes monochromatic wave forcing and linear dynamics, enabling closed-form impedance-matched optimal control. + Standard linear frequency-domain WEC analysis assumes unconstrained linear dynamics, enabling closed-form impedance-matched optimal control. However, real WECs are subject to (i) dynamic constraints (e.g., generator force, power, and stroke limits) and (ii) nonlinear forces, particularly viscous drag. Time-domain simulation handles both naturally but is computationally expensive, especially in an optimization loop. - Pseudo-spectral methods handle constraints by collocating dynamic equations within an optimization, but require an outer numerical optimizer for control synthesis. + Pseudo-spectral methods are nonlinear and handle constraints by collocating the dynamics at specific time-steps, but require an outer numerical optimizer for both simulation and control synthesis. \fi -MDOcean adopts a linearized pseudo-spectral approach that retains the speed of frequency-domain analysis while handling both constraints and drag nonlinearities. -This formulation seeks to preserve the computational advantages of classical frequency-domain analysis while extending its applicability to practical WEC operating regimes that include saturation, drag, and dynamic constraints. -Two techniques are combined: (a) describing functions, which quasi-linearize a nonlinear waveform by retaining only its fundamental harmonic, and (b) an analytical quadratically-constrained quadratic program (QCQP) for the constrained linear controller, exploiting the low dimensionality of the monochromatic problem. -The analytical structure of the monochromatic problem permits many constrained-control solutions to be obtained without iterative numerical optimization, substantially reducing computational cost in large parametric studies. -\Cref{fig:venn-diagram} positions these techniques relative to existing methods. +MDOcean adopts a quasi-linearized pseudo-spectral (QLPS) approach that retains the speed of frequency-domain analysis while handling both temporal constraints and drag nonlinearities. +This formulation seeks to preserve the computational advantages and simplicity of classical linear analysis while extending its applicability to practical WEC operating regimes that include saturation, drag, and dynamic constraints. +QLPS combines two techniques: (a) describing functions, which quasi-linearize a nonlinear waveform by retaining only its fundamental harmonic, +and (b) constrained spectral optimal control, which solves the optimal control problem as a quadratically-constrained quadratic program (QCQP) assuming linear dynamics and spectral-domain constraints. +\Cref{fig:venn-diagram} positions QLPS relative to existing methods. +Augmenting spectral optimal control with describing functions extends its scope to approximate time-domain constraints and nonlinear dynamics. +Furthermore, because QLPS assumptions under monochromatic forcing lead to a monochromatic problem, the QCQP becomes univariate for systems with a single control degree of freedom. +The optimal constrained linear controller can then be found analytically through a geometric interpretation that exploits the low dimensionality. +Its nonlinear counterpart can be reconstructed semi-analytically with inverse describing functions. +Avoiding iterative numerical optimization substantially reduces computational cost in large parametric studies and produces insight into the mathematical structure of dynamic tradeoffs to inform PTO sizing. \begin{figure}[htbp] \centering - \resizebox{0.5\linewidth}{!}{% - \includegraphics{figs/from-matlab/control_method_venn_diagram.pdf} - } - \caption{Venn diagram comparing different methods for control synthesis, with MDOcean's linearized pseudo-spectral approach at the intersection of frequency domain, nonlinear, and constrained optimal control. -See text for a description of each acronym.} + \includegraphics[width=.9\linewidth]{figs/from-matlab/control_method_venn_diagram.pdf} + \caption{Venn diagram comparing different methods for control synthesis, with MDOcean's quasi-linearized pseudo-spectral approach at the intersection of frequency domain, nonlinear, and constrained optimal control. + See text for acronym descriptions.} \label{fig:venn-diagram} \end{figure} @@ -248,80 +263,80 @@ \subsubsection{linearized pseudo-spectral Techniques} These methods are not optimal control methods and do not incorporate constraints, although it is possible to sub-optimally enforce constraints by applying a saturation nonlinearity to the constrained quantity. Finally, methods which combine frequency domain synthesis with constrained optimal control include constrained loop shaping (often enforcing performance constraints while minimizing the distance to an optimal loop shape, in this case the complex-conjugate impedance); linear spectral optimal control (OC), which poses a trajectory optimization problem for a linear system represented with global basis functions in the spectral domain; constrained impedance matching (IM), the narrow-band, frequency-domain equivalent of spectral optimal control; and convex control formulated using linear matrix inequalities (LMIs), typically in the context of robust control \cite{vanantwerp_tutorial_2000}. - At the intersection of all three categories lies nonlinear spectral OC, pseudo-spectral OC, and the method used in MDOcean, which we call ``linearized pseudo-spectral OC.'' + At the intersection of all three categories lies nonlinear spectral OC, pseudo-spectral OC, and the method used in MDOcean, which we call ``quasi-linearized pseudo-spectral OC.'' The pseudospectral method uses the same global basis function representation as spectral OC, but instead of enforcing the dynamics at all points in the trajectory, it does so only at specific collocation points, which enables the representation of time-domain constraints and nonlinearities as opposed to harmonic constraints and nonlinearities \cite{elnagar_pseudospectral_1995}. Pseudo-spectral methods often employ Legendre or Chebyshev polynomial bases, although in wave energy Fourier bases are more typical. In this case, the method is equivalent to a nonlinear trajectory optimization problem that enforces both spectral-domain harmonic balance and time-domain constraints. - The linearized pseudo-spectral method used in MDOcean borrows the pseudo-spectral concept of incorporating time-domain constraints in the spectral domain, but retains only the fundamental harmonic to collapse to the frequency domain. + The quasi-linearized pseudo-spectral method used in MDOcean borrows the pseudo-spectral concept of incorporating time-domain constraints in the spectral domain, but retains only the fundamental harmonic to collapse to the frequency domain. This is exactly analogous to the way that describing functions linearize the harmonic balance method. It can also be thought of as the combination of describing functions with constrained impedance matching. \fi -\Cref{fig:mod-freq-domain-synthesis,fig:mod-freq-domain-evaluation} distinguish two use cases for the linearized pseudo-spectral method: control synthesis (the design of a potentially nonlinear controller to obey constraints and maximize performance of the linearized system) and evaluation (simulation of an arbitrary controller's power performance, where the controller, the plant, or both are nonlinear). +\Cref{fig:mod-freq-domain-synthesis,fig:mod-freq-domain-evaluation} distinguish two application workflows for the QLPS method: control synthesis (the design of a potentially nonlinear controller to obey constraints and maximize performance of the quasi-linearized system) and evaluation (simulation of an arbitrary controller's power performance, where the controller, plant, or both are nonlinear). -MDOcean uses the linearized pseudo-spectral method for both control synthesis and evaluation. -In prior wave energy studies (see \Cref{tab:lit}), related linearized pseudo-spectral approaches have primarily been used for evaluation, often with simplified treatments of constraints. -In contrast, the present work applies the same framework to control synthesis, enabling nonlinear controller design within a frequency-domain formulation. +MDOcean uses QLPS for both control synthesis and evaluation. +In prior wave energy studies, related approaches have been used for evaluation, with simplified treatments of constraints (see F+ and PI+ entries in \Cref{tab:lit}). +We refer to these as the modified frequency-domain family of methods, of which QLPS is the most accurate and optimal. +In contrast, the present work applies QLPS and other modified frequency-domain methods to control synthesis, enabling nonlinear controller design within a frequency-domain formulation. This differs from linearization approaches such as stochastic linearization, which are generally limited to linear controller design. -\begin{figure}[htbp] +\begin{figure*}[htbp] \centering -\caption{Flowchart depicting the linearized pseudo-spectral method with describing function linearization (1)} -\begin{subfigure}[t]{1.1\linewidth} +\caption{Flowchart depicting the QLPS method (M1). The describing functions (DFs) for response/constraint, controller, and plant are shown in yellow, pink, and blue respectively.} +\begin{subfigure}[t]{.9\linewidth} \includegraphics[width=\linewidth]{figs/from-matlab/mod_freq_domain_ctrl_synthesis.pdf} - \caption{Control synthesis - applies to (1.1) only.} + \caption{Control synthesis - applies to (M1.2) only.} \label{fig:mod-freq-domain-synthesis} \end{subfigure} \vfill -\begin{subfigure}[t]{1.1\linewidth} +\begin{subfigure}[t]{.9\linewidth} \includegraphics[width=\linewidth]{figs/from-matlab/mod_freq_domain_ctrl_evaluation.pdf} - \caption{Control evaluation - applies to both (1.1) and (1.2).} + \caption{Control evaluation - applies to both (M1.1) and (M1.2).} \label{fig:mod-freq-domain-evaluation} \end{subfigure} -\end{figure} +\end{figure*} -When handling dynamic constraints with a linearized pseudo-spectral method, first the system is linearized and the unconstrained optimal controller (and the response signals corresponding to that controller) are determined in the frequency domain. +When handling dynamic constraints with a modified frequency-domain method, first the system is (quasi-) linearized and the unconstrained optimal controller (and the response signals corresponding to that controller) are determined in the frequency domain. Then the unconstrained response is checked for constraint violations. This check can be performed in the frequency domain if the time-domain constraint can be adequately linearized (via describing functions) as a harmonic constraint, as shown in the top of \Cref{fig:mod-freq-domain-synthesis}, or the response can be reconstructed in the time domain via the inverse describing function if the constraint is not easily representable via harmonics. If the response for a given sea state does not violate any constraints, then the unconstrained controller and corresponding power is used for that sea state as the frequency-domain linear optimal controller. Otherwise, a variety of approaches are possible to address the violated constraint, summarized in \Cref{tab:constraint-approaches} and described below. -Methods M1 and M2 can be used for both control synthesis and evaluation, while methods M4--M5 can only be used for evaluation since they do not attempt to find a controller that satisfies the constraint. +Methods M1--M2 can be used for both control synthesis and evaluation, while methods M4--M5 can only be used for evaluation since they do not attempt to find a controller that satisfies the constraint. Method M3 can, in general, only be used for evaluation. However, in the special case where the only time-domain constraint is on the control force, the method can be used for control synthesis. \begin{table}[htbp] \centering -\caption{Constraint Satisfaction Approaches for linearized pseudo-spectral Methods} +\caption{Constraint satisfaction approaches for modified frequency domain methods. Filled stars indicate better performance.} \label{tab:constraint-approaches} -\begin{tabular}{>{\raggedright\arraybackslash}p{0.6\linewidth}l l} -\textbf{Method} & \textbf{Accuracy} & \textbf{Optimality} \\ \hline -(M1) Nonlinear optimal (or near-optimal) controller, linearized via describing functions& High & High \\ -(M2) Linear optimal controller& High & Medium \\ -(M3) Saturate or zero the times with constraint violations, neglecting the effect on other signals& Medium & Medium-low\\ -(M4) Saturate or zero entire sea state& Low & Low \\ -(M5) Mark design as infeasible& N/A& Lowest \\ +\begin{tabular}{M{.03\linewidth} >{\raggedright\arraybackslash}m{0.45\linewidth}l l} +&\textbf{Method} & \textbf{Accuracy} & \textbf{Optimality} \\ \toprule +M1 & Nonlinear optimal (M1.1) or near-optimal (M1.2) controller, quasi-linearized via describing functions& \Stars{3} &\Stars{3} \\ \hline +M2 & Linear optimal controller& \Stars{3} & \Stars{2} \\ \hline +M3 & Saturate or zero the times with constraint violations, neglecting the effect on other signals& \Stars{2} & \Stars{1}\\ \hline +M4 & Saturate or zero entire sea state& \Stars{1} & \Stars{0.5} \\ \hline +M5 &Mark design as infeasible& N/A& \Stars{0} \\ \end{tabular} \end{table} % --- M1 --- -The linearized pseudo-spectral approach with the highest accuracy and optimality (M1) is to synthesize a new nonlinear optimal (or near-optimal) time-domain controller that obeys the constraint, taking into account the linearized effect of that nonlinear controller on the frequency domain response. +\paragraph{M1: QLPS with nonlinear control} +The modified frequency-domain approach with the highest accuracy and optimality (M1) is to synthesize a new nonlinear optimal (or near-optimal) time-domain controller that obeys the constraint, taking into account the linearized effect of that nonlinear controller on the frequency domain response. This synthesis can occur in a number of ways. +If a time-domain PMP control law is available for the constraint under consideration (M1.1), it can be used as the nonlinear optimal controller directly, essentially bypassing the use of QLPS in the control synthesis step (\cref{fig:mod-freq-domain-synthesis}) and using it only for evaluation (\cref{fig:mod-freq-domain-evaluation}). \ifdefined\DISSERTATION - If a time-domain PMP control law is available for the constraint under consideration (M1.1), it can be used as the nonlinear optimal controller directly, essentially bypassing the use of the linearized pseudo-spectral linearization in the control synthesis step (\cref{fig:mod-freq-domain-synthesis}) and using it only for evaluation (\cref{fig:mod-freq-domain-evaluation}). This differs from the standard PMP method in that the evaluation of the response and power corresponding to the time-domain PMP controller is determined in the frequency domain, rather than via numerical time-integration (as in the method of \cite{zou_optimal_2017}) or analytical time-integration (as in the method of \cite{lin_fast_2025}). This is advantageous if the computational cost or mathematical labor of the numerical or analytical PMP evaluation methods, respectively, are of concern. - For situations without an explicit PMP control law (M1.2), a nonlinear near-optimal controller can be constructed by first deriving the optimal linear controller via constrained optimization, and then adjusting it to be nonlinear using insights from signal saturation and filtering, which is the approach pursued here for the force limit. - This approach (M1.2) should be considered the primary formulation of the linearized pseudo-spectral method, and the other approaches (M2)--(M5) in \Cref{tab:constraint-approaches} act as approximations of this method with varying levels of accuracy and optimality. - In either case (M1.1 and M1.2), the nonlinear controller is linearized via the describing function method. -\else - For situations without an explicit PMP control law (M1.2), a nonlinear near-optimal controller can be constructed by first deriving the optimal linear controller via constrained optimization - and then adjusting it to be nonlinear using insights from signal saturation and filtering. - This is the approach pursued here for the force limit. - The nonlinear controller is linearized via the describing function method. \fi +For situations without an explicit PMP control law (M1.2), a nonlinear near-optimal controller can be constructed by first deriving the optimal linear controller via constrained optimization +(M2) and then adjusting it to be nonlinear using insights from signal saturation and filtering. +This is the approach pursued here for the force limit. +M1.2 should be considered the primary formulation of QLPS, and the other modified frequency-domain methods (M2)--(M5) in \Cref{tab:constraint-approaches} act as approximations of QLPS with varying levels of accuracy and optimality. +M1.1 and M1.2 both linearize the nonlinear controller via the describing function method. In addition to the control nonlinearity, a describing function is also utilized to linearize the plant drag nonlinearity. % --- M2 --- +\paragraph{M2: QLPS with linear control} If a nonlinear controller or the corresponding describing function to linearize it are not available, the most optimal linear controller that obeys the constraint can be used instead (M2). This approach accurately evaluates the maximum power that a linear controller could produce, although this power will be lower (less optimal) than that of the optimal nonlinear controller. MDOcean pursues this approach for amplitude limits in the operational design load case. @@ -336,14 +351,15 @@ \subsubsection{linearized pseudo-spectral Techniques} \fi % --- M3, M4, M5 --- +\paragraph{M3--M5: Constrained frequency-domain estimates} If neither a numerical nor an analytical implementation is available to identify the optimal constrained linear controller, then the violating portion of the time-domain signals can be saturated or zeroed for the sake of power calculation (M3). This approach neglects the corresponding effect of the saturation/zeroing on other quantities in the system, which lowers the accuracy, and the controller does not have the opportunity to even approximately consider the constraint, resulting in a power production potentially much lower than the constrained optimal. \ifdefined\DISSERTATION This simple approach does not require any derivations to implement. - It has been used in several previous WEC optimization studies for various constraints \cite{garcia-teruel_design_2022,garcia-teruel_reliability-based_2021,cotten_multi-objective_2022,mccabe_constrained_2013} and is evidently the most common linearized pseudo-spectral method in the marine energy field. + It has been used in several previous WEC optimization studies for various constraints \cite{garcia-teruel_design_2022,garcia-teruel_reliability-based_2021,cotten_multi-objective_2022,mccabe_constrained_2013} and is evidently the most common quasi-linearized pseudo-spectral method in the marine energy field. \else This simple approach has been used in several previous WEC optimization studies for various constraints in the studies \cite{garcia-teruel_design_2022,garcia-teruel_reliability-based_2021,cotten_multi-objective_2022,mccabe_constrained_2013} - and is the most common linearized pseudo-spectral method in the marine energy field. + and is the most common quasi-linearized pseudo-spectral method in the marine energy field. \fi It is utilized here for the power limit. @@ -389,7 +405,7 @@ \subsubsection{linearized pseudo-spectral Techniques} The former encodes the describing function formula directly into the simulation (see derivation in \cref{sec:drag}), while the latter implements it indirectly using knowledge of the limit cases of describing functions for multiple nonlinear control laws (see discussion in \cref{sec:appendix-force-sat}). \ifdefined\DISSERTATION - This is possible because the describing function linearized pseudo-spectral method is used for both controller synthesis and evaluation, so the time-domain nonlinear optimal controller that is output from the right of \cref{fig:mod-freq-domain-synthesis} is fed into the top left of \cref{fig:mod-freq-domain-evaluation}, resulting in a chaining of the inverse control describing function with the forward control describing function, which cancel out. + This is possible because the describing function quasi-linearized pseudo-spectral method is used for both controller synthesis and evaluation, so the time-domain nonlinear optimal controller that is output from the right of \cref{fig:mod-freq-domain-synthesis} is fed into the top left of \cref{fig:mod-freq-domain-evaluation}, resulting in a chaining of the inverse control describing function with the forward control describing function, which cancel out. Unless intentionally made different to study robustness, the linearized plant model is identical for the synthesis and evaluation steps. In this case, two workflows essentially become one, with no further computation required after the optimal control step. Reconstruction of the nonlinear controller is not necessary for the evaluation of the power corresponding to that controller. @@ -405,7 +421,7 @@ \subsubsection{linearized pseudo-spectral Techniques} %%%%%%%%%%%%%%% \subsubsection{Equation of Motion}\label{sec:eom} \ifdefined\DISSERTATION - Before applying the linearized pseudo-spectral method described above, we first present the equations of motion for a standard linear WEC model, starting in the time domain and quickly moving to the frequency domain. + Before applying the quasi-linearized pseudo-spectral method described above, we first present the equations of motion for a standard linear WEC model, starting in the time domain and quickly moving to the frequency domain. Specifically, we combine the bi-conjugate network model of \cite{coe_co-design_2025} with the underactuated multi-degree of freedom model of \cite{faedo_principle_2022} to arrive at a formulation suitable for the impedance matching of a multibody WEC subject to wave forcing, powertrain kinematics and dynamics, and drag. \else The MDOcean dynamic model combines the bi-conjugate network model of \cite{coe_co-design_2025} with the underactuated multi-degree of freedom model of \cite{faedo_principle_2022} to express the equations of motion for a multibody WEC subject to wave forcing, powertrain kinematics and dynamics, and drag. @@ -475,7 +491,7 @@ \subsubsection{Equation of Motion}\label{sec:eom} \begin{figure}[htbp] \centering - \includegraphics{figs/from-matlab/circuit_intrinsic.pdf} + \includegraphics[width=\linewidth]{figs/from-matlab/circuit_intrinsic.pdf} \caption{Multiport circuit for two hydrodynamic degrees of freedom with intrinsic and powertrain impedances $\mathbf{Z}_i$ and $\mathbf{Z}_p$} \label{fig:multiport-circuit-intrinsic} \end{figure} @@ -507,10 +523,10 @@ \subsubsection{Power Production}\label{sec:power} In the frequency domain assuming sinusoidal waveforms, the active power $P$, reactive power $Q$, complex power $S$, and apparent power $|S|$ from electrical engineering \cite{saadat_power_1999} are \begin{equation}\label{eq:power-PQS} \begin{aligned} - P &= \tfrac{1}{2}\Re(\hat{e}\hat{q}^*) = \tfrac{1}{2}\Re(Z)|\hat{q}|^2, - & Q &= \tfrac{1}{2}\Im(\hat{e}\hat{q}^*) = \tfrac{1}{2}\Im(Z)|\hat{q}|^2, \\ - S &= P + iQ = \tfrac{1}{2}\hat{e}\hat{q}^* = \tfrac{1}{2}Z|\hat{q}|^2, - & |S| &= \tfrac{1}{2}|Z||\hat{q}|^2, + P &= \tfrac{1}{2}\Re(\hat{e}\hat{q}^*) = \tfrac{1}{2}\Re(Z)|\hat{q}|^2,\\ + Q &= \tfrac{1}{2}\Im(\hat{e}\hat{q}^*) = \tfrac{1}{2}\Im(Z)|\hat{q}|^2, \\ + S &= P + iQ = \tfrac{1}{2}\hat{e}\hat{q}^* = \tfrac{1}{2}Z|\hat{q}|^2, \\ + |S| &= \tfrac{1}{2}|Z||\hat{q}|^2, \end{aligned} \end{equation} where the port is terminated with complex impedance $Z$ and $(\,)^*$ denotes the conjugate transpose. @@ -592,11 +608,11 @@ \subsubsection{Optimal Control}\label{sec:optimal-control} \end{bmatrix} \end{equation} The matrices for various metrics are given in \cref{tab:power-metrics}. -\begin{table}[htbp] +\begin{table}[bhp] \centering \caption{Matrices for calculating various metrics at a port via \cref{eq:metrics}} \label{tab:power-metrics} -\begin{tabular}{l|ccccc} +\begin{tabular}{>{\raggedright\arraybackslash}p{0.1\linewidth}|ccccc} Metric $M$ & $p_{\text{avg}}=P$ & $Q$ & $S$ & $|\hat{e}|^2$ & $|\hat{q}|^2$ \\ \hline Matrix $\mathbf{A}_M$ & $\begin{bmatrix} @@ -716,15 +732,14 @@ \subsubsection{Optimal Control}\label{sec:optimal-control} \centering \begin{subfigure}[t]{0.48\linewidth} \includegraphics[width=\linewidth]{figs/from-matlab/saturation_desc_fcn.pdf} - \caption{Saturated sin vs. time (red) and its fundamental amplitude (blue).}\label{fig:sat} + \caption{Generator torque saturation}\label{fig:sat} \end{subfigure} \hfill \begin{subfigure}[t]{0.48\linewidth} \includegraphics[width=\linewidth]{figs/from-matlab/drag_desc_fcn.pdf} - \caption{Conceptual demonstration of the drag describing function in the time domain. -The nonlinear signal (red) is approximated by its fundamental amplitude (blue).}\label{fig:drag-df} + \caption{Quadratic drag}\label{fig:drag-df} \end{subfigure} -\caption{Describing function approximations for force saturation (a) and drag (b).}\label{fig:desc-fcns} +\caption{Conceptual demonstration of describing function approximations. Plots show force versus time. Nonlinear signals in red; fundamental amplitudes in blue.}\label{fig:desc-fcns} \end{figure} Details of the calculation are provided in \cref{sec:appendix-force-sat}. The result is updated effective values for control gains $B_l$ and $K_l$ that represent the nonlinear saturated controller with a torque that obeys the constraint and is used in place of \cref{eq:constrained-qp-solution}. @@ -749,13 +764,18 @@ \subsubsection{Drag}\label{sec:drag} To improve computational efficiency and allow solution in the frequency domain, the nonlinearity is quasi-linearized with a describing function, following \cite{quartier_influence_2021}. The pressure along the bottom surface of each body due to drag $\vec{p}_d(y,t)$ is modeled with empirical drag coefficient vector $\vec{C}_d$: \begin{equation}\label{eq:drag-pressure} + \vec{p}_{d}(y,t) = \frac{1}{2}\rho_w \vec{C}_d ~\vec{v}_{rel}(y,t) |\vec{v}_{rel}(y,t)| +\end{equation} +\Cref{eq:drag-pressure-approx} expands the velocity term then applies the describing function approximation: +\begin{equation}\label{eq:drag-pressure-approx} \begin{aligned} - \vec{p}_{d}(y,t) &= \frac{1}{2}\rho_w \vec{C}_d ~\vec{v}_{rel}(y,t) |\vec{v}_{rel}(y,t)| \\ - &= \frac{1}{2}\rho_w \vec{C}_d ~|\vec{\hat{V}}_{rel}(y)|^2 \cos(\omega t+\angle \vec{\hat{V}}_{rel}(y)) \left| \cos(\omega t+\angle \vec{\hat{V}}_{rel}(y))\right|\\ - &\approx \frac{1}{2}\rho_w \vec{C}_d ~|\vec{\hat{V}}_{rel}(y)|^2\frac{8}{3\pi} \cos(\omega t+\angle \vec{\hat{V}}_{rel}(y)) + \left(\vec{v}_{rel}(y,t)~\cdot\right.\\ + \left.\left|\vec{v}_{rel}(y,t)\right|\right) + &= |\vec{\hat{V}}_{rel}(y)|^2 \cos(\omega t+\angle \vec{\hat{V}}_{rel}(y)) + \left| \cos(\omega t+\angle \vec{\hat{V}}_{rel}(y))\right|\\ + &\approx~|\vec{\hat{V}}_{rel}(y)|^2\frac{8}{3\pi} \cos(\omega t+\angle \vec{\hat{V}}_{rel}(y)) \end{aligned} \end{equation} -where the last line of \cref{eq:drag-pressure} applies the describing function approximation. \ifdefined\DISSERTATION Computation of the relative velocity uses the following expression for the incident wave velocity phasor in finite depth water evaluated at the body draft, with $\vec{T}=[T_{f,2},T_s]$ collecting the bottom draft of the float and spar: \begin{equation}\label{eq:wave-velocity-phasor} @@ -813,11 +833,11 @@ \subsubsection{Energy Production} To find the long-term average power production $\overline{P}_{\text{elec}}$ in a location with a given distribution of sea states, the power matrix $\mathbf{P}^{H,T}_{\text{elec}}$ for each sea state is weighted by that sea state's probability using a Joint Probability Density (JPD) matrix and then summed. This method is illustrated in \cref{fig:JPD-multiply}. -The JPD used in this analysis represents wave conditions in Humboldt Bay, CA and is taken from the study \cite{janzou_system_2020}. +The JPD used in this analysis represents wave conditions in Humboldt Bay, CA and is taken from the study by \cite{janzou_system_2020}. \begin{figure}[b!] \centering -\includegraphics[width=0.75\linewidth]{figs/manual/jpd_mult.pdf} +\includegraphics[width=\linewidth]{figs/manual/jpd_mult.pdf} \caption{Power matrix multiplication} \label{fig:JPD-multiply} \end{figure} @@ -943,20 +963,23 @@ \subsubsection{Dynamic Limits and Design Load Cases} Incorporating semi-analytical hydrodynamic coefficients for surge is a possible area of future work. Maximum acceptable heave %and surge forces will be evaluated with structural factors of safety in \cref{sec:structures}. -The structures model currently lacks the ability to account for surge forces, so the surge force calculation is not used in the optimization. +The structures model currently lacks the ability to account for surge forces, so the surge force calculation is not used in the simulation. The table also indicates the maximum permissible amplitudes. In operational sea states, the geometric clearance between the float and spar is enforced to prevent overtravel. The permissible heights for upward and downward motion of the float relative to the spar are: \begin{equation}\label{eq:h-fs-up-down} - h_{fs,\text{up}} = h_s - T_s - (h_f- T_{f,2}), \qquad h_{fs,\text{down}} = T_s - T_{f,2} - h_d +\begin{aligned} + h_{fs,\text{up}} &= h_s - T_s - (h_f- T_{f,2}) \\ + h_{fs,\text{down}} &= T_s - T_{f,2} - h_d +\end{aligned} \end{equation} Additionally, a limit of $h_{fs,\text{clear}}$ (see definition in \cref{fig:dims}) is required to ensure clearance of the tubular structure connecting the float and spar. The conditions for linear wave-body interaction are enforced to maintain compatibility with the linear potential flow theory approach of \cref{sec:hydro}. This maximum amplitude $\xi_{\text{linear}}$ for the float and spar respectively is \begin{equation} - \xi_{f,\text{linear}} = (h-T_{f,2})/10, \qquad \xi_{s,\text{linear}} = (h-T_s)/10 + \xi_{f,\text{linear}} = \frac{1}{10}(h-T_{f,2}), \quad \xi_{s,\text{linear}} = \frac{1}{10}(h-T_s) \end{equation} \ifdefined\DISSERTATION which is derived by requiring that the height of the water column below the body (the relevant dimension determining the hydrodynamic coefficients in MEEM) does not change by more than 10\% from its equilibrium value when the body is at its peak displacement. @@ -1025,7 +1048,7 @@ \subsubsection{Irregular Waves}\label{sec:irregular-waves} % Approximation accuracy discussion. Trim closing clause for AOR. In the absence of dynamic constraints, the regular wave assumption does not significantly affect power production if the ultimate design utilizes a high-order linear controller capable of impedance matching over the bandwidth of each sea state. If the ultimate design instead utilizes a simple PI controller, the assumption slightly overestimates the power produced in a broadband wave environment. -A previous study using a Bretschneider spectrum found the irregular wave power generation for a point absorber with PI control to be around 92\% of the perfectly matched power \cite{coe_practical_2021}, which is sufficiently close to justify the approximation for early design phases without worrying about the complexity tradeoff of controller order. +\cite{coe_practical_2021} used a Bretschneider spectrum and found the irregular wave power generation for a point absorber with PI control to be around 92\% of the perfectly matched power, which is sufficiently close to justify the approximation for early design phases without worrying about the complexity tradeoff of controller order. \ifdefined\DISSERTATION Nonetheless, irregular time-domain transient peaks and the variation of hydrodynamic coefficients over the spectral width of each individual sea state are not considered. @@ -1033,7 +1056,7 @@ \subsubsection{Irregular Waves}\label{sec:irregular-waves} \Cref{sec:discussion} will discuss possible future extensions to incorporate irregular wave effects into the analysis. \else Transient peaks and spectral coefficient variation within sea states are not captured, so MDOcean likely underestimates sensitivity to peak force, power, and amplitude constraints. - Stochastic linearization is a candidate future extension; see \cref{sec:unmodeled-effects}. + Stochastic linearization is a candidate for future extension; see \cref{sec:unmodeled-effects}. \fi %TC:break Structures @@ -1053,7 +1076,7 @@ \subsection{Structures}\label{sec:structures} the endurance limit assesses high-cycle fatigue. \fi -The factor of safety is directly calculated for the most heavily loaded subcomponent in each of the three components: the bottom plate in the float, the cylindrical shell of the spar, and the annular plate of the damping plate. +The factor of safety is directly calculated for the most heavily loaded subcomponent in each of the three components: the bottom plate in the float, the cylindrical shell of the spar, and the annular disc of the damping plate. Thicknesses of other subcomponents are either held constant or scaled with the dimensions of subcomponents that are directly assessed, as outlined in \cref{tab:struct-calc-vs-scale}. \ifdefined\DISSERTATION In the future, extending the structures module to directly calculate the factor of safety for each subcomponent would allow for more realistic structural design. @@ -1063,16 +1086,16 @@ \subsection{Structures}\label{sec:structures} \centering \caption{Structural analysis methodology by subcomponent} \label{tab:struct-calc-vs-scale} -\begin{tabular}{ccc} +\begin{tabular}{M{0.15\linewidth} M{0.4\linewidth} M{0.3\linewidth}} \textbf{Component}& \textbf{Structural Subcomponent}& \textbf{Method of Analysis}\\ \hline - \multirow{5}{*}{\centering Float}& Bottom trapezoidal plate& Directly calculated\\ - & Top trapezoidal plate& Scaled from bottom plate\\ - & Radial rectangular plates& Scaled from bottom plate\\ - & Circumferential rectangular plates& Scaled from bottom plate\\ - & PTO connection tubes& Held constant\\ \hline + \multirow[c]{5}{*}{\centering Float}& \begin{itemize}\item Bottom trapezoidal plate\end{itemize}& Directly calculated\\ + & \begin{itemize}\item Top trapezoidal plate\end{itemize}& Scaled from bottom plate\\ + & \begin{itemize}\item Radial rectangular plate\end{itemize}& Scaled from bottom plate\\ + & \begin{itemize}\item Circumferential rectangular plates\end{itemize}& Scaled from bottom plate\\ + & \begin{itemize}\item PTO connection tubes\end{itemize}& Held constant\\ \hline Spar& Cylindrical shell& Directly calculated\\ \hline - \multirow{2}{*}{\centering Damping plate}& Annular plate& Directly calculated\\ - & Spar connection tubes& Held constant\\ + \multirow[c]{2}{*}{\centering Damping plate}& \begin{itemize}\item Annular plate\end{itemize}& Directly calculated\\ + & \begin{itemize}\item Spar connection tubes\end{itemize}& Held constant\\ \end{tabular} \end{table} @@ -1080,7 +1103,7 @@ \subsection{Structures}\label{sec:structures} \begin{figure}[t!] \centering -\includegraphics[width=0.5\linewidth]{figs/from-matlab/FBD.pdf} +\includegraphics[width=.9\linewidth]{figs/from-matlab/FBD.pdf} \caption{Applied loads and fixed points of each structure}\label{fig:FBDs} \end{figure} @@ -1177,4 +1200,4 @@ \subsection{Economics}\label{sec:econ} The per-WEC unit costs and prices $C_{\text{fixed}}$, $C_{pto,\text{constant}}$, $p_{P}$, $p_{F}$, and $p_{s}$ all decrease with the number of devices $N_{WEC}$ via a power law (\cref{eq:cost-power-law} in \cref{sec:appendix-econ}), with curve-fit parameters given in \cref{tab:econ-model-values}. The total $CAPEX$ and $OPEX$ are found by multiplying the unit costs by the number of devices $N_{WEC}$. -For consistency with prior references \cite{neary_reference_2014,RM3}, all costs are intentionally kept in units of 2012 USD without adjusting for inflation. +For consistency with prior references \cite{neary_reference_2014,RM3}, all costs are intentionally kept in units of 2012 USD without adjusting for inflation. \ No newline at end of file diff --git a/pubs/applied-ocean-research-model/sections/other-appendices.tex b/pubs/applied-ocean-research-model/sections/other-appendices.tex index 1cb3f5ac..2ff1128a 100644 --- a/pubs/applied-ocean-research-model/sections/other-appendices.tex +++ b/pubs/applied-ocean-research-model/sections/other-appendices.tex @@ -2,27 +2,40 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Additional Hydrodynamics Details} \label{sec:appendix-additional-hydro} -The hydrodynamic coefficients for the spar are calculated by fitting and interpolating data from WAMIT and reference \cite{olaya_hydrodynamic_2015}. -The data show the excitation and damping coefficients experiencing ``rebound'' effects where rather than slowly asymptoting to zero like those of a single cylinder, they go quickly to zero, increase to a peak, then slowly asymptote down to zero again. -The coefficients have the following form: -\begin{equation} +\subsection{Spar Coefficients} +As \Cref{sec:hydro-other} introduces, because of the importance of the damping plate, the spar hydrodynamic coefficients do utilize other methods besides MEEM. +We approximate the spar added mass coefficient $A_s$ as a frequency-independent quantity using a formula proposed by \cite{philip_damping_2012} based on the displaced water volume in the spherical projection of the damping plate: +\begin{equation}\label{eq:As} + A_s =\frac{1}{3}\rho_w D_d^3\left(1 - \frac{3}{4}r^2\sqrt{1-r^2} - \frac{1}{4}(1-\sqrt{1-r^2})^2(2+\sqrt{1-r^2})\right), +\end{equation} +where $r = D_s/D_d$. +This expression smoothly transitions from a pure plate limit ($\frac{1}{3}\rho_w D_d^3$ at $r=0$) to a pure column limit ($\frac{1}{6}\rho_w D_d^3$ at $r=1$). + +The other spar coefficients are calculated by fitting and interpolating data from WAMIT for the nominal design, then adjusting for known dimensional scalings based on differences from the nominal design. +To derive the dimensional scalings, we use the numerical results that \cite{olaya_hydrodynamic_2015} present for a variant of MEEM that incorporates the damping plate. +By plotting the data for a variety of geometries on a single plot and scaling the axes by dimensionless geometric ratios until the points begin to overlap as a single trend, +we empirically arrive at the following form for the excitation coefficient magnitude: +\begin{equation}\label{eq:scale-wamit} |\gamma_s| = \rho g \pi D_s^2/4 |H_0(kR_x)| \exp(-k e_1) kh \alpha \max(1,\alpha) \sqrt{\frac{f(\frac{kh R_x \alpha^3}{\beta R_p})} {\beta(1+\frac{\beta}{\alpha})}} \end{equation} -\Cref{fig:olaya-hydro-data} shows the nondimensionalized data used for interpolation. +where $f(\cdot)$ indicates a nonlinear relationship to be interpolated from data, and $\alpha$, $\beta$, $e_1$, $R_p$, and $R_x$ are dimensions defined in \citet{olaya_hydrodynamic_2015}. +\Cref{fig:olaya-hydro-data} shows this nonlinear relationship, which has collapsed remarkably for the variety of $\alpha$ and $\beta$ values. \begin{figure}[htbp] \centering - \includegraphics{figs/from-matlab/case4_v1_auto_semilogx_khRxRpa4b_fv1.pdf} + \includegraphics[width=.6\linewidth]{figs/from-matlab/case4_v1_auto_semilogx_khRxRpa4b_fv1.pdf} \caption{Nondimensionalized hydrodynamic data used to interpolate spar coefficients, showing the rebound effect.} \label{fig:olaya-hydro-data} \end{figure} +The dip at mid frequencies shows the excitation and damping coefficients experiencing ``rebound'' effects where rather than slowly asymptoting to zero like those of a single cylinder, they go quickly to zero, increase to a peak, then slowly asymptote down to zero again. +\subsection{Positive-definiteness} Computing the float coefficients in a different way than the spar and coupling coefficients introduces the possibility of violating physical requirements of the solution. In particular, the damping matrix must be positive definite because it represents a dissipative force where outgoing waves remove energy from the system, and the added mass matrix must be positive definite because it represents the kinetic energy of the fluid which must be non-negative. \begin{equation}\label{eq:open-loop-stabilty} - \det(\mathbf{A_{h,op}})=A_f A_s - A_c^2>0 \,, \quad - \det(\mathbf{B_{h,op}})=B_f B_s - B_c^2>0 + \det(\mathbf{A}_{h,op})=A_f A_s - A_c^2>0 \,, \quad + \det(\mathbf{B}_{h,op})=B_f B_s - B_c^2>0 \end{equation} This positive-definiteness may be violated for certain combinations of coefficients and can cause the system dynamics to become unstable. Even when other forces such as body inertia, drag, friction, and control stabilize the system, violation of the conditions of \cref{eq:open-loop-stabilty} can artificially inflate power production. @@ -91,7 +104,7 @@ \subsubsection{Multiport Circuit Conventions} \subsubsection{Time-Domain Expansion of Instantaneous Power} \label{sec:appendix-power-time} -Equation~\eqref{eq:power-PQS} can be used to rewrite the instantaneous power $p(t)$ at any port as \cite{saadat_power_1999} +Equation~\eqref{eq:power-PQS} can be used to rewrite the instantaneous power $p(t)$ at any port as \citep{saadat_power_1999} \begin{equation}\label{eq:power-time-PQS} p(t) = P + |S|\cos\!\left(2(\omega t + \angle\hat{e}) - \tan^{-1}\!\left(\frac{Q}{P}\right)\right), \end{equation} @@ -171,7 +184,7 @@ \subsubsection{PTO Kinematics and Underactuation} In the circuit representation with two hydrodynamic and one PTO degree of freedom, the hybrid matrix \cref{eq:h-matrix-kinematics} can be represented as a three-port element, as shown in \cref{fig:multiport-circuit-kinematics}. \begin{figure}[htbp] \centering - \includegraphics{figs/from-matlab/circuit_kinematics.pdf} + \includegraphics[width=.7\linewidth]{figs/from-matlab/circuit_kinematics.pdf} \caption{Multiport circuit with powertrain kinematics represented as a hybrid matrix, followed by the Th\'{e}venin equivalent of the intrinsic dynamics and kinematics} \label{fig:multiport-circuit-kinematics} \end{figure} @@ -191,7 +204,9 @@ \subsubsection{Combining Intrinsic Dynamics with PTO Kinematics} Although, as previously mentioned, a true cascade matrix must have the same number of input and output ports, it is possible to formulate a cascade-like matrix that can transform from the PTO port to the body ports. However, it cannot transform in the other direction due to underactuation. -While this matrix cannot transform to the body ports directly due to the wave excitation, body velocities can be nonzero even when the PTO force and velocity are both zero. + +Even in the favorable direction, this matrix cannot transform to the body ports directly. +Due to the wave excitation, body velocities can be nonzero even when the PTO force and velocity are both zero. Hence, there is no way to express $\hat{\dot{\xi}}$ as a weighted sum of $\hat{F}_{PTO}$ and $\dot{\hat{X}}_{PTO}$. $\hat{\dot{\xi}}$ could be expressed as a function of the PTO variables by substituting $\mathbf{Z}_p=\mathbf{T}_{kin}^T (\vec{\hat{F}}_{PTO}/\vec{\hat{\dot{X}}}_{PTO}) \mathbf{T}_{kin}$ into \cref{eq:eom-freq-domain}, but this yields a nonlinear expression of limited utility. Instead, we can use cascade-like matrices to solve for a modified variable $\vec{\dot{\xi}}_\Delta$ defined as the change in the body velocities caused by the PTO, in other words the difference between the actual body velocities $\vec{\dot{\xi}}$ and their values if the PTO were disconnected with $F_{PTO}=0$, $\vec{\dot{\xi}}_{0}$: @@ -338,7 +353,7 @@ \subsubsection{PTO Dynamics} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Linear Solution Procedure} \label{sec:appendix-dynamics-solution} -\Cref{sec:eom,sec:pto-dynamics-overview} introduced linear models for the rigid body dynamics, PTO kinematics, and PTO dynamics using impedance, hybrid, and cascade matrices respectively. +\Cref{sec:eom,sec:pto-dynamics-overview} introduce linear models for the rigid body dynamics, PTO kinematics, and PTO dynamics using impedance, hybrid, and cascade matrices respectively. Here we unify the formulation and present a clear system-level solution procedure for all port variables: \begin{enumerate} \item Use \cref{eq:thevenin-electrical} to assemble the Th\'{e}venin source voltage and impedance $\hat{V}_{s,th}$ and $Z_{s,th}$ as seen from the port on which power should be maximized (in this case, the electrical generator port $\{V,I\}$). @@ -357,7 +372,7 @@ \subsection{Indirect Optimal Control Solution Procedure} \label{sec:appendix-qp-numerical} The indirect optimal control approach applies optimality conditions and then numerically solves the resulting equations. The optimality conditions depend on the active set of constraints. -In the most straightforward case where just one constraint exists, the constraint is active if and only if the unconstrained solution violates that constraint. +In the most straightforward case where simply one constraint exists, the constraint is active if and only if the unconstrained solution violates that constraint. For example, if only a maximum force constraint exists, the constrained solution has a force $\hat{F}_{p,\text{constr}}$ that is the minimum of the unconstrained force $\hat{F}_{p,\text{unconstrained}}$ and the maximum allowed force fundamental $F_{\text{max}}$: \begin{equation} \hat{F}_{p,constr} = \min\left(\hat{F}_{p,unconstrained}, F_{max}\right) @@ -370,13 +385,13 @@ \subsection{Indirect Optimal Control Solution Procedure} When multiple constraints exist, the situation is more complex, since it is possible for a given constraint to be violated by the unconstrained solution but not active in the constrained solution, or for multiple constraints to be active simultaneously. In the special case where the unconstrained solution violates no more than one constraint, that constraint is active and the solver can proceed as described above. When the unconstrained solution violates multiple constraints, the solver must try all subsets of 1-2 of those constraint(s) as the active set, find which of these active sets produces the largest power while satisfying all constraints, and then set the error as the deviation of the guessed solution from the constrained solution for that active set. -Since the active set cannot be calculated a-priori as in the simpler case, this approach ends up looking more like a numerical optimization than a numerical root-finding, although the optimization is over a finite set of possibilities instead of a continuous domain. -The main advantage of the indirect approach is its ability to handle nonlinear constraints, although for the quadratic constraints considered here, the analytical approach described in \cref{sec:appendix-qp-solution} is more efficient. +Since the active set cannot be calculated a-priori as in the simpler case, this approach ends up looking more like a numerical optimization than a numerical root-finding, although the optimization is over a finite set of active set possibilities instead of a continuous control domain as it would be in a direct optimal control method. +The main advantage of the indirect approach is its ability to handle arbitrary nonlinear\footnote{Nonlinear meaning greater than second order in an optimization sense, not in a dynamic sense, because the dynamics are already quasi-linearized here} constraints, although for the quadratic constraints considered here, the analytical approach described in \cref{sec:appendix-qp-solution} is more efficient. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Analytical Quadratic Program for Constrained Power Maximization} \label{sec:appendix-qp-solution} -This section describes the solution process for the constrained optimal control problem \cref{eq:opt-problem} introduced in \cref{sec:optimal-control}. +This section describes the analytical solution process for the constrained optimal control problem \cref{eq:opt-problem} introduced in \cref{sec:optimal-control}. %%%%%%%%%%%%%%% \subsubsection{Quadratic Constraint Coefficients} @@ -495,7 +510,7 @@ \subsubsection{Reduced Dimension Formulation} Since \cref{eq:obj-gamma} shows that power decreases with increasing $|\Gamma|$, the power maximization problem \cref{eq:opt-problem} can be recast as minimizing the norm of $\Gamma$: \begin{subequations}\label{eq:opt-problem-gamma} \begin{alignat}{2} -& \min_{\Gamma \in \mathbb{C}} & |\Gamma|& \\[3pt] +& \min_{\Gamma \in \mathbb{C}} & |\Gamma|^2& \\[3pt] & \text{s.t.} & S_{\mu}|\Gamma - \Gamma_{c,\mu}|^2 &\le S_{\mu} r_{\mu}^2, \quad \forall \mu \label{eq:disc-constraint}\\ @@ -588,14 +603,14 @@ \subsubsection{Solution} $a_{\mu\nu} = \frac{1}{2}(r_{\mu}^2-r_{\nu}^2+d_{\mu\nu}^2)/d_{\mu\nu}$ is the distance from the center of the $\mu$th circle to the midpoint of the two intersection points, and $1i$ is the imaginary unit. -The closest-point candidates $\mathcal{C}_{cl}$ match \cref*{eq:circle-intersection}~4 of \cite{mccabe_force-limited_2024}, +The closest-point candidates $\mathcal{C}_{cl}$ match \equationautorefname~4 of \cite{mccabe_force-limited_2024}, which provides the solution for problems with only one constraint. Importantly, we need not compute the entire feasible space $\mathcal{F}$, only compute $\mathcal{C}$ and then check whether each candidate point is feasible. \Cref{fig:qp-solution-geometry} illustrates the various cases for feasibility and optimality for a toy problem with $r_{\mu}=4$, $S_{\mu}=1$, and various $\Gamma_{c,\mu}$. \begin{figure}[htbp] \centering - \includegraphics[width=\textwidth]{figs/from-matlab/qp_circles.pdf} + \includegraphics[width=.8\textwidth]{figs/from-matlab/qp_circles.pdf} \caption{Visualization of the optimization problem \cref{eq:opt-problem-gamma} in the complex plane of $\Gamma$. Toy problems with various constraint circle centers $\Gamma_{c,\mu}$, all with radius $r_{\mu}=4$ and $S_{\mu}=1$, demonstrate the various solution cases in \cref{eq:gamma-opt-quadratic}.} \label{fig:qp-solution-geometry} @@ -609,7 +624,7 @@ \subsubsection{Solution} \subsubsection{Quartic Apparent Power Constraints} \label{sec:appendix-mag-S-constraints} While constraints on effort and flow magnitudes and real and reactive power are quadratic in $\vec{x}$, constraints on the magnitude of apparent power $|S|<|S_{\text{max}}|$, and therefore the maximum and minimum instantaneous power, are fourth-order in $\vec{x}$ and thus $\Gamma$. -This results in a ``Cassini oval'' curve on the complex plane of $\Gamma$, which consists of two closed curves when $|S_{\text{max}}|<|S_0|$ and one curve when $|S_{\text{max}}|\geq|S_0|$, where $|S_0|$ is the apparent power of the unconstrained optimum \cite{weisstein_cassini_2026}. +This results in a ``Cassini oval'' curve \citep{weisstein_cassini_2026} on the complex plane of $\Gamma$, which consists of two closed curves when $|S_{\text{max}}|<|S_0|$ and one curve when $|S_{\text{max}}|\geq|S_0|$, where $|S_0|$ is the apparent power of the unconstrained optimum. Though the optimization is no longer a quadratically constrained quadratic program, the same geometric solution procedure applies, adding circle-oval intersections $\mathcal{C}_{int,\text{oval}-\text{circ}}$, oval-oval intersections, and oval closest points $\mathcal{C}_{cl,\text{oval}}$ to the set of candidate points $\mathcal{C}$. The new candidate points are: \begin{equation} @@ -629,8 +644,9 @@ \subsubsection{Quartic Apparent Power Constraints} \end{aligned} \end{equation} where $\mu$ indexes circle constraints and $k$ indexes oval constraints. -While the earlier circle-circle intersections \cref{eq:candidate-points} result in a quadratic equation that is solved analytically for the coordinates of $\Gamma$, oval-circle intersections lead to quartic equations that are solved numerically with standard polynomial root-finding routines. -We currently assume there is only one apparent power constraint, so oval-oval intersections are not derived, but if present would likewise yield a fourth-order polynomial in $\Gamma$. +While the earlier circle-circle intersections of \Cref{eq:candidate-points} result in a quadratic equation that is solved analytically for the coordinates of $\Gamma$, oval-circle intersections lead to quartic equations that are solved numerically with standard polynomial root-finding routines. +Typically there is only one apparent power constraint, so oval-oval intersections are not derived, but if present would likewise yield a fourth-order polynomial in $\Gamma$. +Therefore, the presence of apparent power constraints makes the optimal control method semi-analytical rather than analytical. %%%%%%%%%%%%%%% \subsubsection{Extension to Higher Dimensions and Convexity} @@ -645,15 +661,15 @@ \subsubsection{Extension to Higher Dimensions and Convexity} The key practical implication is that any locally optimal solution is also globally optimal, so gradient-based solvers converge to the true optimum without exhaustive search, and solution time scales polynomially with the number of decision variables rather than exponentially. After reducing to the $\Gamma$ parameterization, the power objective $|\Gamma|^2$ is a sum of squares and therefore convex. -Even though the original power matrix $\mathbf{A}_P$ (see \cref{tab:power-metrics}) has eigenvalues $\{-1,1\}$ and is indefinite, \cite{bacelli_numerical_2015} shows that substituting the linear dynamics equality constraint into the objective---which is exactly what the reduction to $\Gamma$ in \cref{eq:obj-gamma} accomplishes---renders the objective convex for both regular and irregular wave forcing. +Even though the original power matrix $\mathbf{A}_P$ (see \cref{tab:power-metrics}) has eigenvalues $\{-1,1\}$ and is indefinite, \citet{bacelli_numerical_2015} show that substituting the linear dynamics equality constraint into the objective---which is exactly what the reduction to $\Gamma$ in \cref{eq:obj-gamma} accomplishes---renders the objective convex for both regular and irregular wave forcing. Nonlinear dynamics can reintroduce non-convexity and should be checked on a case-by-case basis. The inequality constraints are more problematic. A constraint is convex if the set of points satisfying it forms a convex region. In the $\Gamma$ plane, a convex constraint corresponds to a circle whose feasible region is the \emph{interior} ($S_\mu=1$ in \cref{eq:disc-constraint}), while a non-convex constraint corresponds to a circle whose feasible region is the \emph{exterior} ($S_\mu=-1$). -Constraints on the effort and flow magnitude (e.g., force limits $|\hat{F}|\leq F_{\text{max}}$) can be convex or non-convex, depending on the complex angle of the Th\'{e}venin equivalent impedance of the relevant source impedance and on the ratio of the constraint limit to its unconstrained value. -Reference \cite{mccabe_force-limited_2024} visualizes the feasible regions for families of such constraints on the complex $\Gamma$ plane. -It shows that the constraints are always convex for purely real impedances; that for reactive plants, nonconvexities occur when the value of the constraint limit exceeds some threshold; and that this threshold lowers as the reactivity increases. +Constraints on the effort and flow magnitude (e.g., force limits $|\hat{F}|\leq F_{\text{max}}$) can be convex or non-convex, depending on the complex angle of the Th\'{e}venin equivalent of the relevant source impedance and on the ratio of the constraint limit to its unconstrained value. +\citet{mccabe_force-limited_2024} visualize the feasible regions for families of such constraints on the complex $\Gamma$ plane. +They show that the constraints are always convex for purely real impedances; that for reactive plants, nonconvexities occur when the value of the constraint limit exceeds some threshold; and that this threshold lowers as the reactivity increases. The quadratic matrices $\mathbf{Q}$ in \cref{tab:qp-constraints} confirm this: $\mathbf{A}_{|\hat{q}|^2}$ and $\mathbf{A}_{|\hat{e}|^2}$ are only positive \emph{semi}definite, indicating constraints that define half-spaces rather than bounded ellipsoids. For the 2D problem studied here, non-convexity is handled by the exhaustive enumeration of candidate points in \cref{eq:gamma-opt-quadratic}, which is efficient in two dimensions but does not scale to higher dimensions. Future work should investigate semidefinite relaxations of these non-convex constraints using the S-procedure \cite{vanantwerp_tutorial_2000}, which replaces each non-convex constraint with a conservative but convex outer approximation. @@ -662,13 +678,13 @@ \subsubsection{Extension to Higher Dimensions and Convexity} %%%%%%%%%%%%%%% \subsubsection{Comparison to Other WEC Constrained Control Formulations} Notably, $\Gamma$ is equal to the reflection coefficient at the generator electrical port, which is commonly used in microwave circuit design to quantify impedance mismatch. -References \cite{mccabe_force-limited_2024,coe_co-design_2025} further discuss the reflection coefficient for WECs. +\citet{mccabe_force-limited_2024,coe_co-design_2025} further discuss the reflection coefficient for WECs. Similar graphical interpretations of constraints for wave energy converters are presented in the studies \cite{mccabe_force-limited_2024,merigaud_geometrical_2023,bacelli_geometric_2013}, with all three visualizing quadratic constraints as circles on the 2D plane. -Reference \cite{mccabe_force-limited_2024} focuses on force/current limits but discusses amplitude, phase voltage, and apparent power constraints as well. -Reference \cite{merigaud_geometrical_2023} considers exclusively amplitude limits while Ref.~\cite{bacelli_geometric_2013} covers amplitude and force limits. -Reference \cite{mccabe_force-limited_2024} uses Smith charts to show constraints directly in $\Gamma$ space as we do here, while Ref.~\cite{merigaud_geometrical_2023} works in the complex plane of the hydrodynamic transmission coefficient and Ref.~\cite{bacelli_geometric_2013} in that of the PTO force. -Reference \cite{bacelli_geometric_2013} also extends the formulation to capture higher harmonics (nonsinusoidal forces), where the exact constraint cannot be easily represented geometrically but sufficient conditions for constraint violation and satisfaction in the higher dimensional space can be represented as hyperspheres and hyperellipsoids. +\citet{mccabe_force-limited_2024} focus on force/current limits but discuss amplitude, phase voltage, and apparent power constraints as well. +\citet{merigaud_geometrical_2023} consider exclusively amplitude limits while \citet{bacelli_geometric_2013} cover amplitude and force limits. +\cite{mccabe_force-limited_2024} use Smith charts to show constraints directly in $\Gamma$ space as we do here, while \citet{merigaud_geometrical_2023} work in the complex plane of the hydrodynamic transmission coefficient and \citet{bacelli_geometric_2013} in that of the PTO force. +\citet{bacelli_geometric_2013} also extend the formulation to capture higher harmonics (nonsinusoidal forces), where the exact constraint cannot be easily represented geometrically but sufficient conditions for constraint violation and satisfaction in the higher dimensional space can be represented as hyperspheres and hyperellipsoids. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -702,7 +718,7 @@ \subsection{Constraint Sensitivity Analysis via Parametric Programming} In our case, we have a semi-analytical solution for the optimal reflection coefficient at each sea state in terms of the constraint centers and radii at that sea state, $\Gamma_{opt}(r_{\mu}, \Gamma_{c,\mu})$, so we can compute the Lagrange multipliers as follows: -\begin{equation} +\begin{equation}\label{eq:lagrange} \begin{aligned} \lambda_\mu &\equiv @@ -725,6 +741,7 @@ \subsection{Constraint Sensitivity Analysis via Parametric Programming} \right) \end{aligned} \end{equation} +The computation of each of the five partials in the final expression of \Cref{eq:lagrange} is discussed next. %%%%%%%%%%%%%%% \subsubsection{Effect of Constraint Limits on Circle Radii} @@ -749,7 +766,7 @@ \subsubsection{Effect of Constraint Limits on Circle Radii} , \quad \frac{\partial r_\mu}{\partial b_\mu} = \frac{1}{2\sqrt{d_\mu(b_\mu + c_\mu)}} \end{equation} -Since $r_\mu\in\mathbb{R}$, then $b_\mu + c_\mu$ and $d_\mu$ must have the same sign, so $\frac{\partial r_\mu}{\partial b_\mu}$ is always positive. +Since $r_\mu\in\mathbb{R}$ (circle radii cannot have an imaginary part), then $b_\mu + c_\mu$ and $d_\mu$ must have the same sign, so $\frac{\partial r_\mu}{\partial b_\mu}$ is always positive. In other words, increasing $b_\mu$ always enlarges the circle radius, regardless of whether the constraint represents the interior or exterior of the circle. Meanwhile, \cref{eq:gamma-disc-center} shows that the center $\Gamma_{c,\mu}$ does not depend on $b_\mu$, so $\frac{\partial \Gamma_{c,\mu}}{\partial b_\mu}=0$. Thus, the constraint limit $L_\mu$ does not affect the circle center. @@ -808,8 +825,8 @@ \subsubsection{Optimal Reflection Coefficient as a Function of Constraint Limits This function can be constructed as $f(\Gamma) = \Gamma + \sum_\mu \mathbb{I}_\mu(\Gamma)$, where $\mathbb{I}_\mu(\Gamma)$ is an indicator function that returns $0$ if $\Gamma$ satisfies the $\mu$th constraint and $\infty$ otherwise. The penalty function is convex in $\Gamma$ if and only if all constraints are convex, i.e., if $S_\mu=1 ~\forall \mu$. -Then, the optimal reflection coefficient of each candidate type can be written as the minimum of the candidate points passed through the penalty function, and the overall optimal reflection coefficient is the minimum of those: -\begin{equation}\label{eq:gamma-opt-min} +Then, the optimal reflection coefficient of each candidate type can be written as the minimum of the candidate points passed through the penalty function. +\begin{equation}\label{eq:gamma-opt-min-per-candidate} \begin{aligned} |\Gamma_{cl,\mu^*}| &= \min_\mu \left\{ @@ -818,7 +835,12 @@ \subsubsection{Optimal Reflection Coefficient as a Function of Constraint Limits |\Gamma_{int,\mu^*,\nu^*}| = \min_{\mu,\nu} \left\{ | f\left(\Gamma_{\mu,\nu}(b_\mu, b_\nu)\right) | - \right\} \\ + \right\} +\end{aligned} +\end{equation} +Finally, the overall optimal reflection coefficient $|\Gamma_{opt}|$ is the minimum of all the per-candidate-type minima: +\begin{equation}\label{eq:gamma-opt-min} +\begin{aligned} |\Gamma_{opt}| &= \min \left\{ |\Gamma_{cl,\mu^*}(b_\mu)|, |\Gamma_{int,\mu^*,\nu^*}(b_\mu, b_\nu)| @@ -851,7 +873,6 @@ \subsubsection{Optimal Reflection Coefficient as a Function of Constraint Limits If that were the case, then by the same argument as for the closest points, the minimum of feasible $|\Gamma_{int}|$ would be concave if and only if all constraints are concave. The optimum $|\Gamma_{opt}|$ is the minimum of the feasible $|\Gamma_{cl}|$ and $|\Gamma_{int}|$, so if both of those are concave, then it is also concave. Therefore, the criteria for the concavity of $|\Gamma_{opt}|$ with respect to $b_\mu$ is that all constraints are convex (interior of circles), and the curvature of $\left(\textrm{aff}(b_{\mu}) + \textrm{aff}(b_{\nu})\right) \sqrt{\textrm{quad}(b_{\mu}, b_{\nu})}$ is concave. - As future work, the above curvature analysis should be verified with a disciplined convex programming package such as CVXPY. To facilitate the analysis of the next section, we observe that although the squared optimal reflection coefficient magnitude, $|\Gamma_{opt}|^2$, @@ -931,7 +952,7 @@ \subsubsection{Average Power as a Function of Constraint Limits} \right] \end{aligned} \end{equation} -Expanding under the assumption that each constraint is sometimes active by itself and sometimes active with another constraint, +Expanding under the assumption that each constraint is active by itself in some sea states and active with another constraint in others, we find that $\overline{P}_{\text{elec}}$ is a piecewise function of each constraint limit $b_\mu$, with a quadratic leading order term as well as affine, square root of affine, and affine times square root of quadratic terms. The two Kronecker delta terms are combined with the weighting term $w_\beta$ to create revised weighting terms $\tilde{w}_{1\mu\nu\beta}$ and $\tilde{w}_{2\mu\nu\beta}$ that are non-zero only when the corresponding constraint(s) are active in that sea state: \begin{equation}\label{eq:power-triple-sum} @@ -964,9 +985,7 @@ \subsubsection{Average Power as a Function of Constraint Limits} &\left[ \sqrt{\textrm{aff}_{\mu\nu\beta}(b_{\mu})} + \textrm{quad}_{\mu\nu\beta}(b_{\mu}, b_{\nu}) - \right.+ - \\ - &\left. + + \left(\textrm{aff}_{\mu\nu\beta}(b_{\mu}) + \textrm{aff}_{\mu\nu\beta}(b_{\nu})\right) \sqrt{\textrm{quad}_{\mu\nu\beta}(b_{\mu}, b_{\nu})} \right] @@ -1083,7 +1102,7 @@ \subsection{Drag Model As described in \cref{sec:drag}, a describing function approximation is used to model the form drag. The relative-velocity computation in \cref{eq:drag-pressure} uses the standard finite-depth linear-wave velocity phasor evaluated at body draft, with $\vec{T} = [T_{f,2}, T_s]$ collecting the bottom drafts of the float and spar: \begin{equation}\label{eq:wave-velocity-phasor} - \vec{\hat{V}}_{wave}(y) = \frac{H}{2}\,\frac{gk}{\omega}\, e^{-k\vec{T}}\, e^{-iky}. + \vec{\hat{V}}_{wave}(y) = \frac{H}{2}\,\frac{gk}{\omega}\, \exp\{-k\vec{T}-iky\}. \end{equation} Strip theory integrates the spatially-dependent infinitesimal drag force vector $d\vec{F}_d$ over the wetted surface area $A_{wet}$. @@ -1118,7 +1137,7 @@ \subsection{Drag Model where we use radii vectors $\vec{R}=[D_f/2,D_d/2]$ and $\vec{R}_{in}=[D_{f,in},0]$ to represent the float and spar. $\vec{\alpha}_v(y)$ is the amplitude ratio of the relative velocity and WEC velocity: \begin{equation}\label{eq:drag-alpha-definitions} - \vec{\alpha}_v(y) = \frac{\vec{|\hat{V}}_{rel}(y)|}{|\vec{\dot{\xi}}|} = \sqrt{1 + \vec{r}^{\, 2} e^{-2ky} - 2 \vec{r} e^{-ky} \sin\angle \vec{\hat{\dot{\xi}}}} + \vec{\alpha}_v(y) = \frac{\vec{|\hat{V}}_{rel}(y)|}{|\vec{\dot{\xi}}|} = \sqrt{1 + \vec{r}^{\, 2} \exp\{-2ky\} - 2 \vec{r} \exp\{-ky\} \sin\angle \vec{\hat{\dot{\xi}}}} \end{equation} and $\vec{r} = |\hat{V}_{\text{wave}}|/ |\vec{\dot{\xi}}|$ is the amplitude ratio of the incident wave velocity and WEC velocity. @@ -1231,7 +1250,7 @@ \subsubsection{Comparison to Other Drag Formulations} This means that Ref.~\cite{quartier_influence_2021} captures only the Froude-Krylov-like component of the drag excitation, neglecting the diffraction-like component. Early MDOcean simulations that used the formulation of \cite{quartier_influence_2021} found this approximation to be unacceptable, producing unstable dynamics (negative drag damping coefficients whose absolute value exceeds the radiation damping) at high frequencies and large wave heights in the operational regime. The strip-theory approach presented in \cref{sec:drag} avoids this issue. -Comparisons against WEC-Sim confirm that even without strip theory, assuming the relative velocity is perfectly in phase with the body velocity yields more accurate results than the drag formulation used in the study \cite{quartier_influence_2021}. +Comparisons against WEC-Sim confirm that even without strip theory, assuming that the relative velocity is perfectly in phase with the body velocity (i.e. a drag force of purely damping and in-phase excitation, with no stiffness or out-of-phase excitation) yields more accurate results than the drag formulation used in the study \cite{quartier_influence_2021}. %%%%%%%%%%%%%%% \subsubsection{Iterative Solution} @@ -1240,7 +1259,7 @@ \subsubsection{Iterative Solution} Iteration is chosen here, the same approach used for a frequency-domain drag simulation of floating offshore wind turbines \cite{hall_open-source_2022}. MDOcean provides users with two options: to employ a typical nonlinear root-finding algorithm, or to simply use the equation of motion \eqref{eq:eom} to obtain subsequent iterates. The latter is a form of fixed point iteration where $\vec{\hat{\dot{\xi}}}_k = g(\vec{\hat{\dot{\xi}}}_{k-1})$, where $g(\vec{\hat{\dot{\xi}}})$ is the equation of motion \cref{eq:eom-freq-domain} and $k$ is the iteration index. -It is therefore guaranteed to converge as long as the dynamics are contracting at the solution, which has the following criteria \cite{chicone_contraction_2006}: +It is therefore guaranteed to converge as long as the dynamics are contracting at the solution, which \cite{chicone_contraction_2006} shows has the following criteria : \begin{equation} \left| \frac{\partial} @@ -1264,8 +1283,8 @@ \subsubsection{Optimal Control Condition} \begin{equation} \begin{aligned} 0 &= \frac{\partial }{\partial B_l} \left[p_{avg,VI}\right] -= \frac{\partial}{\partial B_l} \left[\frac{1}{2} B_l |\hat{I}|^2\right] \\ -&= \frac{\partial }{\partial B_l} \left[\frac{1}{2} B_l \left|\frac{\hat{V}_{s,th}}{Z_{s,th}+B_l+K_l/s}\right|^2 \right] += \frac{\partial}{\partial B_l} \left[\frac{1}{2} B_l |\hat{I}|^2\right] += \frac{\partial }{\partial B_l} \left[\frac{1}{2} B_l \left|\frac{\hat{V}_{s,th}}{Z_{s,th}+B_l+K_l/s}\right|^2 \right] \end{aligned} \end{equation} and the typical simplification $\frac{\partial \hat{V}_{s,th}}{\partial B_l}=\frac{\partial Z_{s,th}}{\partial B_l}=0$ must be replaced with the implicit dependence of the Th\'{e}venin equivalent parameters on $B_l$ because $B_l$ affects $\hat{\dot\xi}$ and therefore $\mathbf{B}_{d}$ and $\vec{\gamma}_{d}$ via \eqref{eq:drag-coeffs}. @@ -1344,8 +1363,8 @@ \subsubsection{Case 2: Large Wave Amplitudes} Once \cref{eq:D-slamming} is satisfied, the set of acceptable $\angle\hat{\xi}$ that avoid negative values for both $\xi_{\text{max},\text{slam}}$ and the right hand side of \cref{eq:slamming-squared} are: \begin{equation}\label{eq:slamming-allowed-angles} \begin{aligned} - \angle\hat{\xi} \in &\left[0,~~\arcsin\left(\frac{\Delta z_{slam}}{H/2}\right)-\frac{kD}{2}\right] \\ - \cup & \left[2\pi-\arcsin\left(\frac{\Delta z_{slam}}{H/2}\right)+\frac{kD}{2}, ~~2\pi\right] + \angle\hat{\xi} \in &\left[0,~~\arcsin\left(\frac{\Delta z_{slam}}{H/2}\right)-\frac{kD}{2}\right] + \cup \left[2\pi-\arcsin\left(\frac{\Delta z_{slam}}{H/2}\right)+\frac{kD}{2}, ~~2\pi\right] \end{aligned} \end{equation} At the edge-case $\angle\hat{\xi}$ on the inner boundaries of this interval, only a single amplitude is permitted because the values for $\xi_{\text{max},\text{slam}}$ and $\xi_{\text{min},\text{slam}}$ coincide at $\sqrt{(H/2)^2-\Delta z_{\text{slam}}^2}$. @@ -1360,7 +1379,7 @@ \subsubsection{Implementation} If $\Delta z_{\text{slam}}{\centering\arraybackslash}m{.5\linewidth} | >{\centering\arraybackslash}m{.5\linewidth} } +\begin{tabular}{c m{1em} | M{.47\linewidth} | M{.47\linewidth} } && \multicolumn{2}{c}{Drag Coefficient}\\ && $C_d=0$& $C_d=1$\\ \hline \multirow{2}{*}[-3em]{\rot{Hydro Coeffs}} &\rot{Identical}& @@ -1484,7 +1504,7 @@ \subsubsection{WEC-Sim} \begingroup \begin{figure}[htbp] \centering -\begin{tabular}{c m{1em} | >{\centering\arraybackslash}m{.5\linewidth} | >{\centering\arraybackslash}m{.5\linewidth} } +\begin{tabular}{c m{1em} | M{.47\linewidth} | M{.47\linewidth} } && \multicolumn{2}{c}{Drag Force}\\ && Magnitude& Phase\\ \hline \multirow{2}{*}[-3em]{\rot{Body}} &\rot{Float}& @@ -1500,7 +1520,7 @@ \subsubsection{WEC-Sim} \end{figure} \endgroup \end{landscape} -The reason for this discrepancy remains unclear, since the expected error of the describing function approximation due to higher harmonics would increase with frequency, whereas the observed error is worst at mid-range frequencies. +The reason for this discrepancy remains unclear, since the expected error of the describing function approximation due to higher harmonics would increase with frequency, whereas the observed error is worse at mid-range frequencies. Evidently, amplitude and power results are quite sensitive to the presence of drag, but not very sensitive to the exact magnitude and phase of the drag force. The WEC-Sim simulation uses drag-only Morison elements with a maximum width of one tenth the incident wavelength to ensure accuracy of strip theory, and uses a Fourier transform to extract the fundamental amplitude and phase of the nonlinear drag waveform. Future work should investigate the source of this error and whether it can be reduced by adjusting the drag model. @@ -1509,19 +1529,19 @@ \subsubsection{WEC-Sim} The general application of the describing function technique requires that the fundamental frequency $\omega$ contains the vast majority of the energy, resulting in body displacements of the form $\xi(t) \sim \cos(\omega t+\phi)$. \Cref{fig:wecsim-thd} shows the total harmonic distortion of the float displacement for each sea state in the $C_d=1$ case, which peaks at 1\% in the worst-case sea state and is below 0.5\% for the majority of sea states. This confirms that higher harmonics have a negligible impact on the system and that the describing function method is a valid approach for modeling drag in this context. -The paricular describing function relation used for drag assumes that the drag force has the form $F_d(t) \sim |\cos(\omega t+\phi)|\cos(\omega t+\phi)$. -\Cref{fig:wecsim-drag-waveform} shows the assumed and actual drag force waveform for four representative sea states corresonding to the four corners of the JPD. +The particular describing function relation used for drag assumes that the drag force has the form $F_d(t) \sim |\cos(\omega t+\phi)|\cos(\omega t+\phi)$. +\Cref{fig:wecsim-drag-waveform} shows the assumed and actual drag force waveform for four representative sea states corresponding to the four corners of the JPD. The match is excellent in all four cases, confirming that the drag force is indeed well approximated by this form, even in high-amplitude, high-frequency conditions where the nonlinearity is expected to be strongest. \begin{figure}[htbp] \centering - \includegraphics[width=1.1\linewidth]{figs/from-matlab/position_THD_contour.pdf} + \includegraphics[width=.6\linewidth]{figs/from-matlab/position_THD_contour.pdf} \caption{Total Harmonic Distortion of Float Displacement in WEC-Sim Simulations with Drag} \label{fig:wecsim-thd} \end{figure} \begin{figure}[htbp] \centering - \includegraphics[width=1.1\linewidth]{figs/from-matlab/drag_force_desc_fcn.pdf} + \includegraphics[width=.6\linewidth]{figs/from-matlab/drag_force_desc_fcn.pdf} \caption{Comparison of Assumed and Actual Drag Force Signal Shape} \label{fig:wecsim-drag-waveform} \end{figure} @@ -1532,7 +1552,7 @@ \subsubsection{Reference Model Report} \Cref{fig:report-power-validation} compares these values against those predicted by WEC-Sim and MDOcean. \begin{figure}[htbp] \centering - \includegraphics[width=1.1\linewidth]{figs/from-matlab/wecsim_rpt_multi_drag_on_meem_on__power_mech_unsat.pdf} + \includegraphics[width=1.05\linewidth]{figs/from-matlab/wecsim_rpt_multi_drag_on_meem_on__power_mech_unsat.pdf} \caption{Comparison of mechanical power for RM3 report \cite{RM3} and MDOcean} \label{fig:report-power-validation} \end{figure} @@ -1590,7 +1610,7 @@ \subsection{RM3 Cost Breakdown Structure} \caption{Cost Breakdown Structure} \label{tab:CBS} -\begin{tabular}{>{\raggedright\arraybackslash}p{0.45\linewidth}>{\centering\arraybackslash}p{0.2\linewidth}>{\centering\arraybackslash}p{0.2\linewidth}>{\centering\arraybackslash}p{0.15\linewidth}} +\begin{tabular}{>{\raggedright\arraybackslash}p{0.3\linewidth}>{\centering\arraybackslash}p{0.2\linewidth}>{\centering\arraybackslash}p{0.2\linewidth}>{\centering\arraybackslash}p{0.15\linewidth}} CBS Category& Nominal \% for $N_{WEC}=1$& Nominal \% for $N_{WEC}=100$&Scales with design?\\\hline 1.1 - Development& 26\%& 3\%&No\\ 1.2 - Infrastructure& 6\%& 4\%&No\\ @@ -1631,10 +1651,10 @@ \subsection{Power-Law Scaling with Number of Devices} \begin{tabular}{ >{\raggedright\arraybackslash}p{0.33\linewidth} >{\raggedright\arraybackslash}p{0.1\linewidth} - >{\raggedright\arraybackslash}p{0.25\linewidth} - >{\raggedright\arraybackslash}p{0.25\linewidth}} + >{\raggedright\arraybackslash}p{0.18\linewidth} + >{\raggedright\arraybackslash}p{0.18\linewidth}} CBS Category& Exponent $\beta$&Unit cost at scale $C_{\infty}$&Single-unit cost \newline $C_1$\\\hline - \textbf{1.1-1.3, 1.6-1.9} - Design-Independent CAPEX& 0.741& 1.24 \$M&13.92 \$M\\ + \textbf{1.1-1.3, 1.6-1.9} - Design-Indep. CAPEX& 0.741& 1.24 \$M&13.92 \$M\\ \textbf{1.4} - Structural CAPEX & 0.481& 2.387 \$/kg&4.294 \$/kg\\ \textbf{1.5} - Power Take Off CAPEX& 0.206& Constant: 92.59 \$k Power: 0.4454 \$/W Force: 0.0086 \$/N&Constant: 93.64 \$k Power: 1.355 \$/W Force: 0.0204 \$/N\\ \textbf{2.1-2.6} - OPEX&0.557& 0&1.193 \$M\\ @@ -1675,7 +1695,7 @@ \section{Structures Module Details} It outputs a factor of safety for each limit case. First, it obtains equivalent section properties based on the plate and stiffener dimensions. Second, it relates applied loads to stresses and deflection using analytical solutions to structural boundary value problems obtained from the well-known Roark's handbook \cite{young_roarks_2001} and the references therein. -Finally, it utilizes design standards from organizations like Det Norske Veritas, the American Bureau of Shipping, and the American Iron and Steel Institute to develop limit state expressions for each possible failure mode of each major structural element under each design load case. +Finally, it utilizes design standards from organizations like Det Norske Veritas, the American Bureau of Shipping, and the American Iron and Steel Institute to develop limit-state expressions for each possible failure mode of each major structural element under each design load case. For consistency with the reference model report \cite{RM3}, all structural elements are assumed to be welded, although future work should consider the use of pinned joints for certain interfaces to reduce reaction moments and enhance structural efficiency. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -1708,7 +1728,7 @@ \subsection{Equivalent Section Properties for Stiffened Plates} where $M$ is the maximum moment per unit length and $z_{\text{max}}$ is the maximum distance from the neutral axis of the stiffened section. This parallels the stress for an unstiffened plate of thickness $t$, which is $\sigma=6M/t^2 = 12 M(t/2)/t^3$. -Note that this equivalent-thickness method of accounting for stiffeners assumes that as a whole, the stiffener-plate system deflects like a single plate element, rather than as set of multiple plate elements separated by stiffeners. +Note that this equivalent-thickness method of accounting for stiffeners assumes that as a whole, the stiffener-plate system deflects like a single plate element, rather than as a set of multiple plate elements separated by stiffeners. This is a reasonable assumption when the stiffeners are small and densely spaced with respect to the plate, but breaks down if the stiffeners come to dominate the system. The so-called effective breadth ratio $\rho$ is used to quantify the validity of this assumption, where the equivalent-thickness method requires $\rho=1$: \begin{equation} @@ -1764,12 +1784,12 @@ \subsection{Float} \begin{figure}[htbp] \centering - \includegraphics[width=0.75\linewidth]{figs/from-matlab/trapezoid.pdf} + \includegraphics[width=0.5\linewidth]{figs/from-matlab/trapezoid.pdf} \caption{Trapezoidal float plate with inscribed circle used to determine $x_{\text{min}}$} \label{fig:trapezoid} \end{figure} -To approximately capture the dimensions of a trapezoid which effectively set the deflection curvature and therefore stress, $x_{\text{min}}$ is set to the diameter of the largest circle that can be inscribed in the trapezoid tangent to the shorter base $b_1$. +To approximately capture the dimensions of a trapezoid, which effectively set the deflection curvature and therefore stress, $x_{\text{min}}$ is set to the diameter of the largest circle that can be inscribed in the trapezoid tangent to the shorter base $b_1$. This inscribed circle scheme satisfies all limit cases mentioned above and can be calculated as \begin{equation} x_{\text{min}} = \begin{cases}b_1 \left( m + \sqrt{1+m^2}\right), & h_0 \geq \sqrt{b_1b_2} \\ h_0, & h_0 < \sqrt{b_1b_2}\end{cases} @@ -1805,11 +1825,11 @@ \subsubsection{Float Top Plate} using $\beta_1$ and $\beta_2$ tabulated as a function of $x_{\text{max}}/x_{\text{min}}$ in reference \cite{young_roarks_2001}. The bending moment per unit length is then found as $M_{\text{top}}=\sigma_{\text{edge}}t_{f,\text{top},eq}^2/6$. However, unlike \cref{eq:bottom-plate-moment}, it is found that \cref{eq:top-plate-stress} predicts stresses that do not match well with the finite element results for the nominal RM3 design in references \cite{RM3,previsic_reference_2011}. -It is unknown whether this is an implementation or modeling issue and remains future work to address. +It is unknown whether this is an implementation or modeling issue, and remains as future work to address. In the meantime, the top plate thickness is set by scaling the bottom plate thickness, as described in \cref{tab:struct-calc-vs-scale}. \subsubsection{Other Float Plates} -The other float plates besides the top and bottom ones experience only the edge reaction moments from the top and bottom plates, with no additional external loads in heave. +The other float plates besides the top and bottom ones experience only the edge reaction moments that the top and bottom plates apply at their interface, with no additional external loads in heave. Assuming that the maximum allowable stress is the same in all float plates, the required side plate thickness can be found by simply scaling the input float bottom plate thickness using the float thickness ratios in the nominal RM3 design. This prevents the need for a separate structural analysis of the side plates. @@ -1876,8 +1896,8 @@ \subsection{Column} where $\zeta$ and $\omega$ are defined as follows: \begin{equation} \begin{aligned} - \zeta &= 1 - P_r(1 - P_r)~\frac{\sigma_Y}{F_{\text{crit}}/A} - \frac{\sigma_\theta}{\sigma_Y} \\ - \omega &= \frac{\sigma_\theta}{2\sigma_Y} \left(1 - \frac{\sigma_\theta}{2\sigma_Y}\right) + \zeta &= 1 - P_r(1 - P_r)~\frac{\sigma_Y}{F_{\text{crit}}/A} - \frac{\sigma_\theta}{\sigma_Y}, \qquad + \omega = \frac{\sigma_\theta}{2\sigma_Y} \left(1 - \frac{\sigma_\theta}{2\sigma_Y}\right) \end{aligned} \end{equation} and the hoop stress $\sigma_\theta$ is diff --git a/pubs/shared/els-cas/cas-common.sty b/pubs/shared/els-cas/cas-common.sty index e3e4154b..155ebc61 100644 --- a/pubs/shared/els-cas/cas-common.sty +++ b/pubs/shared/els-cas/cas-common.sty @@ -1986,7 +1986,7 @@ \else \__short_authors: :~ \fi - { \rmfamily \itshape Preprint~ submitted ~to ~Elsevier } + { \rmfamily \itshape Preprint~ submitted ~to ~\journal } \group_end: } diff --git a/pubs/shared/references.bib b/pubs/shared/references.bib index 7c7f0aec..20c8b9e8 100644 --- a/pubs/shared/references.bib +++ b/pubs/shared/references.bib @@ -1,3 +1,13 @@ +@phdthesis{mccabe_dissertation_2026, + address = {Ithaca, NY}, + type = {{PhD} dissertation}, + title = {Leveraging {Semi}-{Analytical} {Modeling}, {Multidisciplinary} {Design} {Optimization}, and {System} {Value} {Metrics} to {Advance} {Wave} {Energy} {Converter} {Viability}}, + url = {https://calkit.io/symbiotic-engineering/mdocean/publications?path=pubs%2Fdissertation%2FsampleThesis.pdf}, + school = {Cornell University}, + author = {McCabe, Rebecca}, + month = aug, + year = {2026}, +} @book{hamming_numerical_1962, address = {New York},