Numerical Robustness Research Articles

An accelerated, high-order accurate direct solver for the Lippmann–Schwinger equation for acoustic scattering in the plane

An efficient direct solver for solving the Lippmann–Schwinger integral equation modeling acoustic scattering in the plane is presented. For a problem with N degrees of freedom, the solver constructs an approximate inverse in $\mathcal {O}(N^{3/2})$ operations and then, given an incident field, can compute the scattered field in $\mathcal {O}(N \log N)$ operations. The solver is based on a previously published direct solver for integral equations that relies on rank-deficiencies in the off-diagonal blocks; specifically, the so-called Hierarchically Block Separable format is used. The particular solver described here has been reformulated in a way that improves numerical stability and robustness, and exploits the particular structure of the kernel in the Lippmann–Schwinger equation to accelerate the computation of an approximate inverse. The solver is coupled with a Nyström discretization on a regular square grid, using a quadrature method developed by Ran Duan and Vladimir Rokhlin that attains high-order accuracy despite the singularity in the kernel of the integral equation. A particularly efficient solver is obtained when the direct solver is run at four digits of accuracy, and is used as a preconditioner to GMRES, with each forwards application of the integral operators accelerated by the FFT. Extensive numerical experiments are presented that illustrate the high performance of the method in challenging environments. Using the 10th-order accurate version of the Duan–Rokhlin quadrature rule, the scheme is capable of solving problems on domains that are over 500 wavelengths wide to relative error below 10− 10 in a couple of hours on a workstation, using 26M degrees of freedom.

A multi-level nonlinear elimination-based JFNK method for multi-scale multi-physics coupling problem in pebble-bed HTR

In the pebble-bed high-temperature gas-cooled reactor, various physical phenomena occur at different spatial scales and are coupled together tightly in a nonlinear manner. More specifically, the local temperature profiles inside the fuel sphere and TRISO particle, which are essential in reactor design and safety analysis, are coupled with the global physical fields, such as the pebble bed temperature, helium temperature, as well as the fission heat source. Therefore, the multi-physics calculation combined with multi-scale model should be employed to describe the reactor behavior. The computational efficiency is the critical issue for the numerical solution of such global and local coupled systems. Jacobian-free Newton Krylov (JFNK) algorithm is an excellent candidate for multi-physics coupling issues due to its high convergence rate and superior computational performance, but it doesn’t consider the physical features at different scales. In this work, a multi-level nonlinear elimination is proposed within the JFNK framework to deal with the multi-scale coupled system in high-temperature gas-cooled reactor. The numerous local multi-scale temperature variables, including the fuel sphere and TRISO particle temperatures, are eliminated efficiently based on the Schur complement factorization due to its well sparse matrix structure. As a result, the number of unknowns is significantly reduced and the computational performance is improved. A steady-state pebble-bed HTR model is used in this work, and the numerical results show that multi-level nonlinear elimination-based JFNK method has high computational efficiency and numerical robustness, whose performance is 4 ∼ 6 times higher than that of the original JFNK algorithm.

A local to global (L2G) finite element method for efficient and robust analysis of arbitrary cracking in 2D solids

This paper presents and validates a new local to global (L2G) FEM approach that can analyze multiple, interactive fracture processes in 2D solids with improved numerical efficiency and robustness. The method features: 1) forming local problems for individual and interactive cracks; and 2) parallel solving local problems and returning local solutions as part of the trial solution for global iteration. It has been demonstrated analytically (through a simple 1D problem) and numerically (through several benchmarking examples) that, the proposed method can substantially improve the robustness of the global solution process and significantly reduce the costly global iteration for convergence. The demonstrated improvement in numerical efficiency is up to 20∼40% for mildly unstable problems. For problems with severely unstable crack initiation and propagation, the improvement can be more significant. This new method is readily applicable to other popular methods such as the extended FEM (X-FEM), Augmented FEM (A-FEM) and Phantom-node method (PNM).

Isogeometric topology optimization of strain gradient materials

In this paper optimal topologies of isotropic linear elastic strain gradient materials are investigated by means of isogeometric topology optimization. The employment of strain gradient theory allows not only to capture the microstructural effects of materials but also to regularize stress/strain concentration phenomena and to address the so-called wedge forces. Isogeometric analysis is used in order to meet the requirements of higher-order continuity of strain gradient elasticity theory. For the purpose of determining the constitutive parameters of strain gradient materials, a novel experiment is conceived by using the Digital Image Correlation (DIC) technique. The so-called SIMP (Solid Isotropic Material with Penalization) method is applied to interpolate the material stiffness tensors (including the fourth-order and the sixth-order stiffness tensor) by power law with penalty. Benchmark numerical experiments are conducted to illustrate the computational effectiveness and numerical robustness of the proposed method.

Stress Optimization of Vent Holes with Different Shapes Using Efficient Switching Delayed PSO Algorithm

An effective integrated design optimization method is developed to reduce the maximum von Mises stress around vent holes of a high-pressure turbine sealing disk. It mainly includes four different shape designs (circular, elliptical, race-track, and four-arc) for holes, an updated self-developed modelling and meshing tool, an APDL-based strength analysis, and a self-proposed efficient switching delayed particle swarm optimization (SDPSO) algorithm. The main idea of SDPSO is: (1) by evaluating an evolutionary factor and utilizing a probability transition matrix, a non-homogeneous Markov chain is determined and auto-updated in each generation; (2) the evolutionary factor and the Markov chain are used to adaptively select the inertia weight, acceleration coefficients, and delayed information to adjust the particle’s velocity. The performance of SDPSO is evaluated through two benchmark optimization problems with constraints. The results show that SDPSO is superior to two well-known PSO algorithms in optimization capability, numerical robustness, and convergence speed. Furthermore, SDPSO is used for the stress optimization of vent holes with four different shapes. The results show that: (1) SDPSO is suitable and valuable for practical engineering optimization problems with constraints; (2) the developed integrated design optimization method is effective and advanced for reducing the maximum von Mises stress around the vent holes; and (3) the four-arc hole has more tremendous advantages in reducing the maximum von Mises stress, followed by the elliptical hole, the race-track hole, and the circular hole.

A high-fidelity solver based on hybrid numerical methods on unstructured grids for incompressible multiphase flows

A high-fidelity solver is developed on the unstructured grids to simulate incompressible multiphase flows with free surface based on volume of fluid (VOF) method. In this framework, the velocity is defined by both volume integrated average (VIA) and point value (PV) as computational variables which are updated simultaneously by volume-average/point-value multi-moment (VPM) method at each time step. Compared with our previous work, a node-based finite element (FE) method is used instead of cell-centered finite volume (FV) method for the discretization of pressure Poisson equation, which provides a new pattern that can naturally treat with the velocity and pressure coupling. It does not require Rhie–Chow interpolation or the use of a Petrov–Galerkin finite element formulation with tunable parameters to circumvent the well-known instability conditions. Moreover, we propose a novel balanced-force formulation with FE discretization to accommodate the surface tension force, which effectively eliminates the large numerical errors associated with the non-orthogonality of unstructured grids. In VOF method, the volume fraction is defined as VIA to identify the moving interface which is updated by solving advection equation through the FV formulation. The resulting model that combines hybrid discretizations, rigorously preserves the mass and momentum conservation, substantially improves the numerical accuracy and robustness of pressure field, and significantly suppresses the spurious velocity to nearly machine error even for multiphase flows with large density ratios. Various numerical examples have been presented to demonstrate the excellent capabilities of present solver which can provide high-fidelity solutions for practical applications in presence of largely deformed interface and highly complex geometry.

An aerial ammunition ad hoc network collaborative localization algorithm based on relative ranging and velocity measurement in a highly-dynamic topographic structure

In the internet of battlefield things, ammunition is becoming networked and intelligent, which depends on location information. Therefore, this paper focuses on the self-organized network collaborative localization of munitions with an aerial three-dimensional (3D) highly-dynamic topographic structure under a satellite denied environment. As for aerial networked munitions, the measurement of munitions is objectively incomplete due to the degenerated and interrupted link of munitions. For this reason, a cluster-oriented collaborative localization method is put forward in this paper. Multidimensional scaling (MDS) was first integrated with a trilateration localization method (TLM) to construct a relative localization algorithm for determining the relative location of a mobile cluster network. The information related to relative velocity was then combined into a collaborative localization framework to devise a TLM-vMDS algorithm. Finally, an iterative refinement algorithm based on scaling by majorizing a complicated function (SMACOF) was employed to effectively eliminate the influence of incomplete link observation on localization accuracy. Compared with the currently available advanced algorithms, the proposed TLM-vMDS algorithm achieves higher localization accuracy and faster convergence for a cluster of extensively networked munitions, and also offers better numerical stability and robustness for high-speed motion models.

A stabilized mixed-FE scheme for frictional contact and shear failure analyses in deformable fractured media

Simulation of fracture contact mechanics in deformable fractured media is of paramount important in computational mechanics. Previous studies have revealed that compressive loading may produce mode II fractures, which is quite different from mode I fractures induced by tensile loading. Furthermore, fractures can cross each other. This will increase the complexity of their network deformation under different loading types significantly. In this work, a stabilized mixed-finite element (FE) scheme with Lagrange multipliers is proposed in the framework of variational formulation, which is able to simulate frictional contact, shear failure (mode II) and opening (mode I) of multiple crossing fractures. A novel treatment is devised to guarantee physical solutions at the intersection of crossing fractures. A preconditioner is introduced to re-scale the saddle-point algebraic system and to preserve the numerical robustness. Then, a solution strategy is designed to calculate the unknowns, displacement and Lagrange multipliers, in one algebraic system. Later, numerical tests are conducted to study mechanical behaviors of fractured media. Benchmark study is performed to verify the presented mixed-FE scheme. A deformable medium with crossing fractures is simulated under mixed-mode loading types. The characteristics of fracture contact, surface sliding, opening and variation of stress intensity factor are analyzed. Simulation results show that the curve of slippage induced by compression, as well as the opening induced by internal fluid pressure, along the fracture length holds a parabolic shape. The diagonal contact point, at the intersecting position of the crossing fractures, is studied in detail, specially under different stress states. Finally, the impact of intersecting fractures on frictional contact mechanics is investigated for different loading conditions.

Robust hybrid/mixed finite elements for rubber-like materials under severe compression

A new family of hybrid/mixed finite elements optimized for numerical stability is introduced. It comprises a linear hexahedral and quadratic hexahedral and tetrahedral elements. The element formulation is derived from a consistent linearization of a well-known three-field functional and related to Simo–Taylor–Pister (STP) elements. For the quadratic hexahedral and tetrahedral elements we derive (static reduced) discontinuous hybrid elements, as well as continuous mixed finite elements with additional primary unknowns for the hydrostatic pressure and the dilation, whereas the linear hexahedral element is of the discontinuous type. The elements can readily be used in combination with any isotropic, invariant-based hyperelastic material model and can be considered as being locking-free. In a representative numerical benchmark test the elements numerical stability is assessed and compared to STP-elements and the family of discontinuous hybrid elements implemented in the commercial finite element code Abaqus/Standard. The new elements show a significant advantage concerning the numerical robustness.

System-level, input–output and new parameterizations of stabilizing controllers, and their numerical computation

It is known that the set of internally stabilizing controller Cstab is non-convex, but it admits convex characterizations using certain closed-loop maps: a classical result is the Youla parameterization, and two recent notions are the system-level parameterization (SLP) and the input–output parameterization (IOP). In this paper, we address the existence of new convex parameterizations and discuss potential tradeoffs of each parameterization in different scenarios. Our main contributions are: (1) We reveal that only four groups of stable closed-loop transfer matrices are equivalent to internal stability: one of them is used in the SLP, another one is used in the IOP, and the other two are new, leading to two new convex parameterizations of Cstab. (2) We investigate the properties of these parameterizations after imposing the finite impulse response (FIR) approximation, revealing that the IOP has the best ability of approximating Cstab given FIR constraints. (3) These four parameterizations require no a priori doubly-coprime factorization of the plant, but impose a set of equality constraints. However, these equality constraints will never be satisfied exactly in floating-point arithmetic computation and/or implementation. We prove that the IOP is numerically robust for open-loop stable plants, in the sense that small mismatches in the equality constraints do not compromise the closed-loop stability; but a direct IOP implementation will fail to stabilize open-loop unstable systems in practice. The SLP is known to enjoy numerical robustness in the state feedback case; here, we show that numerical robustness of the four-block SLP controller requires case-by-case analysis even when the plant is open-loop stable.

Software update: The ORCA program system—Version 5.0

AbstractVersion 5.0 of the ORCA quantum chemistry program suite was released in July 2021. ORCA 5.0 represents a major improvement over all previous versions of ORCA and features (1) highly improved performance, (2) increased numerical robustness, (3) a host of new functionality, and (4) greatly improved user friendliness. The article describes the most salient features of the program.This article is categorized under: Electronic Structure Theory > Ab Initio Electronic Structure Methods Data Science > Computer Algorithms and Programming Software > Quantum Chemistry

A computationally efficient scroll compressor model for both single-phase and two-phase compression considering scroll wrap temperature distribution

A general scroll compressor model with high accuracy, fast simulation speed and good numerical robustness is proposed in this study. The governing ordinary differential equations derived from mass and energy conservation use pressure and specific enthalpy as state variables, which can be applied to both single-phase compression and two-phase compression. Both the governing equations for inner compression and the nonlinear equations for internal leakages are explicit in pressure and specific enthalpy to avoid severe numerical stiff problems. In addition, the detailed scroll wrap temperature distribution is calculated based on the proposed second order differential equation considering heat conduction and periodic heat convection. Validation against experimental data shows that cooling capacity errors are less than 3.5% and power consumption errors are less than 2.8% in non-injection and liquid injection conditions. Evaluation of model speed and robustness shows that the model has a convergence rate of 100% and an average calculation speed of 12.4 s per case in various conditions.

Multiple-precision sparse matrix–vector multiplication on GPUs

Sparse matrix–vector multiplication (SpMV) appears in many application domains, and performance is the key consideration when implementing SpMV kernels. At the same time, accuracy is also important because rounding errors can drastically change the computed result. Multiple-precision arithmetic is a common approach to improve the accuracy of results. In this paper, we implement and evaluate multiple-precision SpMV kernels in the CSR, JAD, ELLPACK, and DIA matrix storage formats for graphics processing units. In the proposed implementation, the matrix is represented in double precision, while the input and output vectors are in multiple precision and internal computations are also performed in multiple precision. Our underlying floating-point arithmetic algorithms are based on the residue number system, which is attractive because of its carry-free nature and provides an arbitrary level of precision determined by the set of moduli that comprise the base of the system. In particular, we apply sets from 8 to 64 moduli, reaching precision values from 106 to 848 bits, and demonstrate how higher precision reduces the rounding errors in the SpMV kernel. We also conduct a thorough analysis of the CSR kernel in terms of roofline performance, occupancy and memory bandwidth, and evaluate the performance and memory consumption for all the proposed kernels. Numerical experiments on matrices from real-world applications show that ELLPACK offers the highest performance in many examples, while JAD is a good trade-off between performance and space requirements. Furthermore, our proposed kernels in general substantially improve the efficiency of sparse matrix–vector multiplication compared to implementations built on top of existing multiple-precision CUDA libraries. Finally, we integrate the multiple-precision SpMV into a preconditioned conjugate gradient linear solver and identify test cases where our implementation exhibits superior convergence and numerical robustness.

An adaptive method of moving asymptotes for topology optimization based on the trust region

A novel gradient-based optimizer, which makes use of Trust Region based Moving Asymptotes and is abbreviated as TRMA method, is proposed for the topology optimization. Compared with the most frequently used optimizer, the method of moving asymptotes (MMA), the step-based subproblem adopted in the TRMA method shows better accuracy in approximating the response surface near the current iterative point. Also, a criterion to judge the quality of the trial step is introduced in the trust region scheme which is regarded as an adaptive process for automatically adjusting the expansion parameters of the moving asymptotes. The above features endow the TRMA method with inherent numerical robustness and efficiency. Three benchmark problems (including the compliance minimization problem and the stress constraint problem) are considered to investigate the performance of the TRMA method. Although those problems are well addressed by many other creative strategies, which always involve modification of the optimization formulations, few studies focus on the development of the algorithms directly. Hence, the inherent complexity of those formulations (like the high nonlinearity introduced by aggregation and the large-scale constraints when the aggregation-free formulation is adopted) is still a great difficulty for the existing optimizers. Results show that the TRMA method is a powerful optimizer for handling problems with high nonlinearity and large-scale constraints and may promote the exploration and application of the topology optimization methods, even when the traditional MMA becomes inefficient.

Revisit to the THINC/QQ scheme: Recent progress to improve accuracy and robustness

AbstractIn this article, we revisit to the THINC/QQ (tangent of hyperbola interface capturing method with quadratic surface representation and Gaussian quadrature) scheme to make thorough investigations on its numerical property. Furthermore, it presents some recent progress made extensively on THINC/QQ algorithm to improve numerical accuracy and robustness which include (1) construct the interface equation on the global coordinate so as to apply THINC/QQ on polygonal/polyhedral elements, which ensures the boundedness of numerical solutions as well; (2) determine the coefficients of reconstruction polynomial from the unit normal vector of interface to reduce computational cost and alleviate numerical diffusion; (3) represent the interface with a cubic polynomial rather than quadratic representation used in THINC/QQ to obtain better geometric fidelity; (4) provide an appropriate choice considering different combinations, such as the number of quadrature points and temporal discretization methods. The verification benchmark tests are carried out on different types of grids, which convinces that THINC/QQ scheme with normalization on global coordinate delivers higher accuracy and robustness than the original one, and demonstrates effective competition compared with other existing interface capturing methods.

Recent Upgrades in a 2D Turbulent Transport Solver Based on a Hybrid Discontinuous Galerkin Method for the Simulation of Fusion Plasma in Tokamak

The simulation of fusion plasmas in realistic magnetic configurations and tokamak geometries still requires the development of advanced numerical algorithms owing to the complexity of the problem. In this context, we propose a Hybrid Discontinuous Galerkin (HDG) method to solve 2D transport fluid equations in realistic magnetic and tokamak wall geometries. This high-order solver can handle magnetic equilibrium free structured and unstructured meshes allowing a much more accurate discretization of the plasma facing components than current solvers based on magnetic field aligned methods associated with finite-differences (volumes) discretization. In addition, the method allows for handling realistic magnetic equilibrium, eventually non steady, a critical point in the modeling of full discharges including ramp up and ramp down phases. In this paper, we introduce the HDG algorithm with a special focus on recent developments related to the treatment of the cross-field diffusive terms, and to an adaptive mesh refinement technique improving the numerical efficiency and robustness of the scheme. The updated solver is verified with a manufactured solution method, and numerical tests are provided to illustrate the new capabilities of the code.

A Modified Fifth-Order WENO-HLLC Riemann Solver for Transonic Viscous Flows around Helicopter Rotor in Hover

This paper presents a modified fifth-order WENO-HLLC Riemann solver for the transonic compressible viscous flows around the helicopter rotor in hover. The HLLC approximate Riemann solver is proposed to discrete the convection term involving grid velocity of the Navier-Stokes equations. In order to solve the interface flow accurately, a modified fifth-order WENO scheme is presented by designing the smoothness indicators based on L 1 norm measurement. The improved WENO scheme can provide the optimal approximation order even at critical points. Numerical accuracy and robustness are validated by several benchmark inviscid flow problems. Then the numerical properties of the WENO-HLLC solver in conjunction with the implicit LU-SGS time integration method with high efficiency are further validated by simulating transonic viscous flows over RAE2822 airfoil and ONERA-M6 wing. The results show that the accuracy of calculating shock, discontinuity, and the vortex is significantly improved. Finally, the method is developed to compute the transonic vortex flow around the helicopter rotor with a domain discretized by overset grids. The results indicate that the proposed method is very robust and effective in acquiring high resolution for vortex wake.

Preliminary study on Reynolds stress model based on <i>ν</i><sub><i>t</i></sub>-scale equation

Reynolds stress model has always been the frontier and challenging problem in turbulence model theory research, where improving numerical robustness is the key to its wide application in engineering. Referring to the classical <inline-formula><tex-math id="M8">\begin{document}$k$\end{document}</tex-math><alternatives><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M8.jpg"/><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M8.png"/></alternatives></inline-formula>-<inline-formula><tex-math id="M9">\begin{document}$kL$\end{document}</tex-math><alternatives><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M9.jpg"/><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M9.png"/></alternatives></inline-formula> turbulence model, a new <inline-formula><tex-math id="M10">\begin{document}${\nu_t}$\end{document}</tex-math><alternatives><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M10.jpg"/><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M10.png"/></alternatives></inline-formula>-scale equation is constructed and used to couple the SSG/LRR model to form a so-called SSG/LRR-<inline-formula><tex-math id="M11">\begin{document}${\nu_t}$\end{document}</tex-math><alternatives><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M11.jpg"/><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M11.png"/></alternatives></inline-formula> Reynolds stress model. Four benchmark cases, including zero pressure gradient turbulent plate boundary layer, airfoil wake flow, supersonic square duck flow and separated flow over NACA0012 airfoil at 45 degree angle of attack, are carried out to test the new turbulence model. At the same time, high-order numerical schemes are used to discretize the turbulence equations in order to assess its numerical robustness. The results are compared with those of SA eddy viscosity model and SSG/LRR-<inline-formula><tex-math id="M12">\begin{document}$\omega$\end{document}</tex-math><alternatives><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M12.jpg"/><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M12.png"/></alternatives></inline-formula> Reynolds stress model. It is shown that the <inline-formula><tex-math id="M13">\begin{document}${\nu_t} $\end{document}</tex-math><alternatives><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M13.jpg"/><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M13.png"/></alternatives></inline-formula>-scale equation is strictly equal to zero at the viscous wall boundary. Compared with the traditional <inline-formula><tex-math id="M14">\begin{document}$\omega $\end{document}</tex-math><alternatives><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M14.jpg"/><graphic xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="16-20220417_M14.png"/></alternatives></inline-formula>-scale, it has better numerical robustness. Along with this, the new model can be matched with the high-order numerical schemes and obtain a better efficiency in the mesh convergence. Moreover, the new model has the inherent advantage of Reynolds stress model in simulating the corner flow and has the potential in scale adaptive simulation of unsteady separated flow.

Discrete-Time Differential Dynamic Programming on SO(3) With Pose Constraints

The motion of many robotics systems, such as the rotation of unmanned aerial vehicles, can be modeled by SO(3). However, the difficulties in the parameterization of the SO(3) makes it hard to implement the state-of-the-art collision avoidance algorithm. In this paper, we present a new method for control on SO(3) via the combination of the nonlinear state constraints and geometric control technique. We first define the nonlinear constraints for trajectory optimization on SO(3). Then we solve the trajectory optimization problem by constrained Differential Dynamic Programming (DDP) in a Riemannian geometry framework. The first and second-order expansion of the cost function and dynamics are derived in a coordinate-free manner. The safety condition represented by nonlinear constraints is incorporated into the DDP via an active set method, where we find the active set in the backward path to obtain the optimal control policy. We validated our methods on motion planning of rigid model on SO(3) manifold where pose constraints are imposed. Our methods outperform the baseline methods in terms of convergence speed and numerical robustness to disturbance. The numerical robustness is essential when the system is initialized far from the local optimum.

Inequality Constrained Optimal Control for Rope-Assisted ASV Docking Maneuvers

In difficult whether conditions the complexity of docking maneuvers of autonomous surface vessels (ASVs) increases. Docking maneuvers take place in proximity to piers and other vessels. The slow speed at which the maneuvers take place decreases the maneuverability of the vessel and environmental disturbances such as wind and waves increase in effect. In these scenarios ropes can be utilized to support the maneuver and compensate for the compromised actuation capabilities. The modeling of the modified mechanics plays a center role. It has been shown in previous work that this setup, connecting the vessel to the pier with a constant distance (e.g. with a rod) can be modeled using the concept of differential-algebraic equations (DAEs). With the adapted vessel dynamics a trajectory generation problem can be solved using optimal control theory. Here, a direct numerical method, in this case a full discretization of the optimal control problem (OCP), is used. The contribution in this paper allows for slackness in the rope. This is done by using a smooth penalty function, as well as solving a complementary problem. Both methods lead to feasible solutions. The penalty function shows are more smooth behavior but suffers from numerical robustness in comparison to the linear complementary solution.