Fourth-order Predictor-corrector Method Research Articles

Computational study on aerodynamically coupled piezoelectric harvesters

In this work, the authors present a two-dimensional computational model for predicting the aeroelastic response as well as the output power of vertically arranged harvesters by taking into account all aerodynamic interactions. The piezo-aeroelastic framework consists of the following: (1) an aerodynamic model based on the unsteady vortex-lattice method to compute the aerodynamic forces; (2) a discrete parameter model for each harvester with 3 degrees of freedom (plunge motion, pitch motion, and the voltage generated by the piezoelectric effect); (3) an inter-model connection to exchange information between models at each time step; and (4) a numerical scheme based on the Hamming’s fourth-order predictor–corrector method to integrate all the governing equations in the time domain. The results obtained allow us to infer new insights into the flutter onset as well as the post-critical behavior of harvester arrangements. An interesting finding is that the flutter speed is significantly decreased as the distance between the harvesters is reduced. The results suggest the strong possibility of effective energy extraction at low flow speeds using properly distributed harvester arrangements. However, in post-critical conditions, the output power is significantly enhanced as the free-stream speed is increased.

Numerical study on the performance of semicircular and rectangular submerged breakwaters

A systematic numerical comparative study of the performance of semicircular and rectangular submerged breakwaters interacting with solitary waves is the basis of this paper. To accomplish this task, Nwogu\'s extended Boussinesq model equations are employed to simulate the interaction of the wave with breakwaters. The finite difference technique has been used to discretize the spatial terms while a fourth-order predictor-corrector method is employed for time discretization in our numerical model. The proposed computational scheme uses a staggered-grid system where the first-order spatial derivatives have been discretized with fourth-order accuracy. For validation purposes, five test cases are considered and numerical results have been successfully compared with the existing analytical and experimental results. The performances of the rectangular and semicircular breakwaters have been examined in terms of the wave reflection, transmission, and dissipation coefficients (RTD coefficients) denoted by K_R, K_T, K_D. The latter coefficient K_D emerges due to the non-energy conserving K_R and K_T,. Our computational results and graphical illustrations show that the rectangular breakwater has higher reflection coefficients than semicircular breakwater for a fixed crest height, but as the wave height increases, the two reflection coefficients approach each other. On the other hand, the rectangular breakwater has larger dissipation coefficients compared to that of the semicircular breakwater and the difference between them increases as the height of the crest increases. However, the transmission coefficient for the semicircular breakwater is greater than that of the rectangular breakwater and the difference in their transmission coefficients increases with the crest height. Quantitatively, for rectangular breakwaters the reflection coefficients K_R are 5-15% higher while the diffusion coefficients K_D are 3-23% higher than that for the semicircular breakwaters, respectively. The transmission coefficients K_T for rectangular breakwater shows the better performance up to 2.47% than that for the semicircular breakwaters. Based on our computational results, one may conclude that the rectangular breakwater has a better overall performance than the semicircular breakwater. Although the model equations are non-dissipative, the non-energy conserving transmission and reflection coefficients due to wave-breakwater interactions lead to dissipation type contribution.

A co-simulation methodology to simulate the nonlinear aeroelastic behavior of a folding-wing concept in different flight configurations

A methodology to simulate the unsteady, nonlinear aeroelastic behavior of a folding-wing concept in multiple flight configurations is presented. It is based on a partitioned or co-simulation scheme, which divides the dynamical system into interacting aerodynamic and structural models. The aerodynamic model that predicts the loads is based on the unsteady vortex-lattice method. The structural model that predicts the motion of the folding wing is based on the finite-element method. The grids in the two models are non-matching. The method for simultaneously integrating the combined set of equations is based on the fourth-order predictor-corrector method developed by Hamming. The technique for transferring information between the non-matching grids is the radial basis interpolation method. This system is partially validated by showing that the predicted loads here agree very closely with numerical results based on the Euler equations for a wing with prescribed unsteady twisting and pitching motion and with analytical solutions available in the literature. Finally, a series of numerical simulations related to the aeroelastic behavior of a folding-wing concept inspired by gull wings provide new insights into flutter boundaries as a function of the dihedral angles of the inner and outer wings. The findings in this paper strongly suggest that the present numerical aeroelastic model will be a valuable computational tool for further studies of aircraft with morphing wings.

Computational Dynamics of Flapping Wings in Hover Flight: A Co-Simulation Strategy

A co-simulation strategy for modeling the unsteady dynamics of flying insects and small birds as well as biologically inspired flapping-wing micro-air-vehicles is developed in this work. In particular, the dynamic system under study is partitioned in two subsystems (the structural model and the aerodynamic model) that exchange information in a strong way. The vehicle or insect system is modeled as a collection of rigid bodies and lifting surfaces that can undergo deformations such as spanwise twisting, in-plane bending, out-of-plane bending, and an arbitrary combination of these deformation mechanisms. To account for the loads associated with the airflow, an aerodynamic model based on an extended version of the unsteady vortex-lattice method is used. The motion equations are integrated by using a fourth-order predictor-corrector method along with a procedure to stabilize the solution of the resulting differential-algebraic equations. The numerical results obtained for the unsteady lift and dynamics of a fruit fly in free hover flight are found to be in close agreement with prior experimental results reported in the literature. Furthermore, the inclusion of an adequate wing deformation pattern results in an increase of the lift force compared with that of a rigid wing surface, pointing to the importance of wing flexibility on aerodynamic performance. From the findings reported in this paper, it is believed that the numerical simulation framework presented here could serve as a computational tool for further studies of flying insects and micro-air-vehicles.

Calculation of fluid–structure interaction: methods, refinements, applications

This paper starts by reviewing the most significant previous studies on fluid–structure interaction of flows past stationary and moving circular cylinders, including those undergoing fluid dynamic excitation in the form of vortex-induced vibrations. It shows the importance of flow past a circular cylinder as a convenient validation test case for new numerical methods and demonstrates the advantages of the immersed boundary method. A new method based on the immersed boundary method to accurately compute fluid–structure interaction and vortex-induced vibrations is then presented and applied. The fluid solver is strongly coupled with the equation of motion through Hamming's fourth-order predictor–corrector method, which allows computation of the interaction of a fluid with a moving structure. The accuracy, validity and efficiency of the proposed method are demonstrated by a series of validation cases, including flow past stationary and moving cylinders, both prescribed and freely vibrating. The numerical method reproduces fairly well the vibration of an elastically mounted cylinder subjected to a fluid flow. In particular, vibration amplitude as well as vibration frequency are computed quite accurately. In addition, the predicted vortex-induced cylinder vibration lock-in regime is in good agreement with laboratory experiments.

Numerical solution for hybrid fuzzy systems by Milne's fourth order predictor-corrector method

In this paper we study numerical methods for hybrid fuzzy differential equations by an appllication of the Milne’s fourth order predictor- corrector method for fuzzy differential equations. We prove a convergence result and give numerical examples to illustrate the theory.

Nonlinear analysis and enhancement of wing-based piezoaeroelastic energy harvesters

We investigate the level of harvested power from aeroelastic vibrations for an elastically mounted wing supported by nonlinear springs. The energy is harvested by attaching a piezoelectric transducer to the plunge degree of freedom. The considered wing has a low-aspect ratio and hence three dimensional aerodynamic effects cannot be neglected. To this end, the three dimensional unsteady vortex lattice method for the prediction of the unsteady aerodynamic loads is developed. A strong coupling scheme that is based on Hamming's fourth-order predictor–corrector method and accounts for the interaction between the aerodynamic loads and the motion of the wing is employed. The effects of the electrical load resistance, nonlinear torsional spring and eccentricity between the elastic axis and the gravity axis on the level of the harvested power, pitch and plunge amplitudes are investigated for a range of operating wind speeds. The results show that there is a specific wind speed beyond which the pitch motion does not pick any further energy from the incident flow. As such, the displacement in the plunge direction grows significantly and causes enhanced energy harvesting. The results also show that the nonlinear torsional spring plays an important role in enhancing the level of the harvested power. Furthermore, the harvested power can be increased by an order of magnitude by properly choosing the eccentricity and the load resistance. This analysis is helpful in designing piezoaeroelastic energy harvesters that can operate optimally at specific wind speeds.

Modeling and performance analysis of cambered wing-based piezoaeroelastic energy harvesters

We investigate the effects of aerodynamic loads on the performance of wing-based piezoaeroelastic energy harvesters. The rigid airfoil consists of pitch and plunge degrees of freedom supported by flexural and torsional springs with a piezoelectric coupling attached to the plunge degree of freedom. The effects of aerodynamic loads are investigated by considering a camber in the airfoil. A two-dimensional unsteady vortex-lattice method (UVLM) is used to model the unsteady aerodynamic loads. An iterative scheme based on Hamming’s fourth-order predictor–corrector method is employed to solve the governing equations simultaneously and interactively. The effects of varying the camber, its location, and the nonlinear torsional spring coefficient are determined. The results show that, for small values of the camber location, the flutter speed changes greatly on increasing the camber of the airfoil. On the other hand, for large values of the camber location, the variation of the flutter speed when changing the camber is very negligible. We demonstrate that the symmetric airfoil case is the best configuration to design enhanced wing-based piezoaeroelastic energy harvesters. Furthermore, the results show that an increase in the camber results in a decrease in the level of the harvested power. For cambered airfoils, we demonstrate that an increase in the camber location leads to an increase in the level of the harvested power. The results show that an increase in the airfoil camber delays the appearance of a secondary Hopf bifurcation.

A FINITE DIFFERENCE METHOD FOR THE ONE-DIMENSIONAL VARIATIONAL BOUSSINESQ EQUATIONS

The variational Boussinesq equations derived by Klopman et. al. (2005) con-verse mass, momentum and positive-definite energy. Moreover, they were shown to have significantly improved frequency dispersion characteristics, making it suitable for wave simulation from relatively deep to shallow water. In this paper we develop a numerica lcode for the variational Boussinesq equations. This code uses a fourth-order predictor-corrector method for time derivatives and fourth-order finite difference method for the first-order spatial derivatives. The numerical method is validated against experimen-tal data for one-dimensional nonlinear wave transformation problems. Furthermore, the method is used to illustrate the dispersive effects on tsunami-type of wave propagation.DOI : http://dx.doi.org/10.22342/jims.14.1.57.1-11

Camber Effects on the Power Harvesting from Piezoaeroelastic Systems

We investigate the effects of the aerodynamic loads on the performance of piezoaeroelastic energy harvesters. The harvester consists of a rigid airfoil having a pitch and plunge degrees of freedom with a piezoelectric coupling attached to the plunge degree of freedom. The Unsteady Vortex Lattice Method is used to model the unsteady flow and predict the loads. An iterative scheme based on Humming’s fourth order predictor-corrector method is employed to solve simultaneously and interactively the governing equations. The effects of varying the airfoil camber coeffcient are determined. We demonstrate that increasing the camber does not necessarily increase the level of the harvested power.

A new form of generalized Boussinesq equations for varying water depth

A new set of equations of motion for wave propagation in water with varying depth is derived in this study. The equations expressed by the velocity potentials and the wave surface elevations include first-order non-linearity of waves and have the same dispersion characteristic to the extended Boussinesq equations. Compared to the extended Boussinesq equations, the equations have only two unknown scalars and do not contain spatial derivatives with an order higher than 2. The wave equations are solved by a finite element method. Fourth-order predictor–corrector method is applied in the time integration and a damping layer is applied at the open boundary for absorbing the outgoing waves. The model is applied to several examples of wave propagation in variable water depth. The computational results are compared with experimental data and other numerical results available in literature. The comparison demonstrates that the new form of the equations is capable of calculating wave transformation from relative deep water to shallow water.

Gravitational Fall Velocity of Sphere in Viscous Fluid

The gravitational transient fall velocity of a rigid sphere in an otherwise quiescent viscous fluid is studied through examining the physical reasonableness of simulated results and comparing the results to Moorman's 1955 free-fall experiments. The published sphere dynamic equations from low-to-moderate sphere Reynolds number are solved numerically by a fourth order predictor-corrector method and iterations on sphere velocity and acceleration. Among the published sphere dynamic expressions, Mei and Adrian's 1992 expression has the best agreement with Moorman's data. The relative importance of the steady and unsteady drags along the gravitational transient process is also discussed. It is found that neglecting the unsteady drag indeed simplifies the computational procedure, but the accuracy on the time-varying fall velocity is significantly compromised. The added mass and history terms are of great importance at early stages of gravitational transient falling for low sphere-to-fluid density ratio and low sphere Reynolds number.

Fast spin dynamics algorithms for classical spin systems

We have proposed new algorithms for the numerical integration of the equations of motion for classical spin systems. In close analogy to symplectic integrators for Hamiltonian equations of motion used in Molecular Dynamics, these algorithms are based on the Suzuki-Trotter decomposition of exponential operators and unlike more commonly used algorithms exactly conserve spin length and, in special cases, energy. Using higher order decompositions we investigate integration schemes of up to fourth order and compare them to a well-established fourth order predictor-corrector method. We demonstrate that these methods can be used with much larger time steps than the predictor-corrector method and thus may lead to a substantial speedup of computer simulations of the dynamical behavior of magnetic materials.

Generalized Gear's method for computing the flow of a viscoelastic fluid

A fourth-order predictor-corrector method is developed for obtaining the numerical solution of a class of singular boundary value problems in which the highest derivative is further multiplied by a parameter which can take arbitrarily small values. First, the method is tested on a prototype linear differential equation. It is then used to compute the two-dimensional stagnation point flow of a viscoelastic fluid. Finally, the flow over a stretching sheet is computed for which an exact solution exists. The comparison of the results shows that the method gives highly accurate results for a moderately sized integration step.

Critical dynamics in the two-dimensional classical XY model: A spin-dynamics study.

Using spin-dynamics techniques we have performed large-scale computer simulations of the dynamic behavior of the classical three component XY model (i.e., the anisotropic limit of an easy-plane Heisenberg ferromagnet), on square lattices of size up to ${192}^{2}$, for several temperatures below, at, and above ${\mathit{T}}_{\mathrm{KT}}$. The temporal evolution of spin configurations was determined numerically from coupled equations of motion for individual spins by a fourth-order predictor-corrector method, with initial spin configurations generated by a hybrid Monte Carlo algorithm. The neutron-scattering function S(q,\ensuremath{\omega}) was calculated from the resultant space-time displaced spin-spin correlation function. Pronounced spin-wave peaks were found both in the in-plane and the out-of-plane scattering function over a wide range of temperatures. The in-plane scattering function ${\mathit{S}}^{\mathit{xx}}$ also has a large number of clear but weak additional peaks, which we interpret to come from two-spin-wave scattering. In addition, we observed a small central peak in ${\mathit{S}}^{\mathit{xx}}$, even at temperatures well below the phase transition. We used dynamic finite-size scaling theory to extract the dynamic critical exponent z. We find z=1.00(4) for all T\ensuremath{\le}${\mathit{T}}_{\mathrm{KT}}$, in excellent agreement with theoretical predictions, although the shape of S(q,\ensuremath{\omega}) is not well described by current theory. \textcopyright{} 1996 The American Physical Society.

Critical dynamics of the body-centered-cubic classical Heisenberg antiferromagnet.

We have used spin dynamics techniques to perform large-scale simulations of the dynamic behavior of the L\ifmmode\times\else\texttimes\fi{}L\ifmmode\times\else\texttimes\fi{}L body-centered-cubic classical Heisenberg antiferromagnet with L\ensuremath{\le}48 at a range of temperatures above and below as well as at the critical point ${\mathit{T}}_{\mathit{c}}$. The temporal evolutions of the spin configurations were determined numerically from coupled equations of motion for individual spins by a fourth-order predictor-corrector method, with initial spin configurations generated by Monte Carlo simulations. The neutron scattering function S(q,\ensuremath{\omega}) was calculated from the space- and time-displaced spin-spin correlation function. We used a previously developed dynamic finite-size scaling theory to extract the dynamic critical exponent z from S(q,\ensuremath{\omega}) at ${\mathit{T}}_{\mathit{c}}$. Our results are in agreement with the theoretical prediction of z=1.5 and with experimental results; however, we find that the asymptotic regime was only entered at L\ensuremath{\approxeq}30. In the analysis of the form of the transverse and longitudinal components of S(q,\ensuremath{\omega}) we found that a central diffusion peak appears below ${\mathit{T}}_{\mathit{c}}$ predominantly in the longitudinal component and remains present through and above ${\mathit{T}}_{\mathit{c}}$. The transverse component of the spin-wave peak is Lorentzian below ${\mathit{T}}_{\mathit{c}}$ but for T\ensuremath{\geqslant}${\mathit{T}}_{\mathit{c}}$ is described best by a more complex functional form. Below ${\mathit{T}}_{\mathit{c}}$ we see evidence of multiple spin-wave peaks in the longitudinal component. \textcopyright{} 1996 The American Physical Society.

SPIN DYNAMICS SIMULATIONS OF CLASSICAL, THREE-DIMENSIONAL HEISENBERG MAGNETS

Spin dynamics simulations have been used to study dynamic critical behavior of classical Heisenberg magnets. The temporal evolutions of the spin configurations were determined numerically from coupled equations of motion by a fourth-order predictor corrector method, with initial spin configurations generated by Monte-Carlo simulations. The neutron scattering function, S (q, ω), was calculated from the space and time displaced spin-spin correlation function and the dynamic critical exponent z was extracted using dynamic finite size scaling theory. For both ferromagnetic and antiferromagnetic models we find good agreement with theoretical predictions and experimental results.

Time-Dependent Numerical Code for Extended Boussinesq Equations

The extended Boussinesq equations derived by Nwogu (1993) significantly improve the linear dispersive properties of long-wave models in intermediate water depths, making it suitable to simulate wave propagation from relatively deep to shallow water. In this study, a numerical code based on Nwogu's equations is developed. The model uses a fourth-order predictor-corrector method to advance in time, and discretizes first-order spatial derivatives to fourth-order accuracy, thus reducing all truncation errors to a level smaller than the dispersive terms retained by the model. The basic numerical scheme and associated boundary conditions are described. The model is applied to several examples of wave propagation in variable depth, and computed solutions are compared with experimental data. These initial results indicate that the model is capable of simulating wave transformation from relatively deep water to shallow water, giving accurate predictions of the height and shape of shoaled waves in both regular and irregular wave experiments.

Spin-dynamics study of the dynamic critical behavior of the three-dimensional classical Heisenberg ferromagnet.

Using spin-dynamics techniques we have performed large-scale computer simulations of the dynamic behavior of the L\ifmmode\times\else\texttimes\fi{}L\ifmmode\times\else\texttimes\fi{}L body-centered-cubic classical Heisenberg ferromagnet with L\ensuremath{\le}40 in the vicinity of the critical point ${\mathit{T}}_{\mathit{c}}$. The temporal evolutions of spin configurations were determined numerically from coupled equations of motion for individual spins by a fourth-order predictor-corrector method with initial spin configurations generated by Monte Carlo methods. The space- and time-displaced spin-spin correlation functions and their space-time Fourier transforms were calculated to determine the neutron-scattering functions. We developed a dynamic finite-size scaling theory for the neutron-scattering function at ${\mathit{T}}_{\mathit{c}}$ and used it to extract the dynamic critical exponent z. Within our resolution limit, the value of z was estimated to be 2.478(28), in excellent agreement with the dynamic scaling prediction, and dynamic scaling was found to be valid for momentum transfer q up to 0.4\ensuremath{\pi} in the (1,0,0) direction.

DYNAMICAL INTEGRATION

In this letter we show how to use a new form of integration, called dynamical integration, that utilizes the dynamics of a system defined by an ODE to construct a map that is in effect a one-step integrator. This method contrasts sharply with classical numerical methods that utilize polynomial or rational function approximations to construct integrators. The advantages of this integrator is that it uses only one step while preserving important dynamical properties of the solution of the ODE: First, if the ODE is conservative, then the one-step integrator is measure preserving. This is significant for a system having a highly nonlinear component. Second, the one-step integrator is actually a one-parameter family of one-step maps and is derived from a continuous transformation group as is the set of solutions of the ODE. If each element of the continuous transformation group of the ODE is topologically conjugate to its inverse, then so is each member of the one-parameter family of one-step integrators. If the solutions of the ODE are elliptic, then for sufficiently small values of the parameter, the one-step integrator is also elliptic. In the limit as the parameter of the one-step family of maps goes to zero, the one-step integrator satisfies the ODE exactly. Further, it can be experimentally verified that if the ODE is chaotic, then so is the one-step integrator. In effect, the one-step integrator retains the dynamical characteristics of the solutions of the ODE, even with relatively large step sizes, while in the limit as the parameter goes to zero, it solves the ODE exactly. We illustrate the dynamical, in contrast to numerical, accuracy of this integrator with two distinctly different examples: First we use it to integrate the unforced Van der Pol equation for large ∊, ∊≥10 which corresponds to an almost continuous square-wave solution. Second, we use it to obtain the Poincaré map for two different versions of the periodically forced Duffing equation for parameter values where the solutions are chaotic. The dynamical accuracy of the integrator is illustrated by the reproduction of well-known strange attractors. The production of these attractors is eleven times longer when using a conventional fourth-order predictor-corrector method. The theory presented here extends to higher dimensions and will be discussed in detail in a forthcoming paper. However, we caution that the theory we present here is not intended as a line of research in numerical methods for ODEs.