Mechanics2DRigidGivenLaw.cpp

Go to the documentation of this file.

/*--------------------------------*- VM2D -*-----------------*---------------*\
| ##  ## ##   ##  ####  #####   |                            | Version 1.14   |
| ##  ## ### ### ##  ## ##  ##  |  VM2D: Vortex Method       | 2026/03/06     |
| ##  ## ## # ##    ##  ##  ##  |  for 2D Flow Simulation    *----------------*
|  ####  ##   ##   ##   ##  ##  |  Open Source Code                           |
|   ##   ##   ## ###### #####   |  https://www.github.com/vortexmethods/VM2D  |
|                                                                             |
| Copyright (C) 2017-2026 I. Marchevsky, K. Sokol, E. Ryatina, A. Kolganova   |
*-----------------------------------------------------------------------------*
| File name: Mechanics2DRigidGivenLaw.cpp                                     |
| Info: Source code of VM2D                                                   |
|                                                                             |
| This file is part of VM2D.                                                  |
| VM2D is free software: you can redistribute it and/or modify it             |
| under the terms of the GNU General Public License as published by           |
| the Free Software Foundation, either version 3 of the License, or           |
| (at your option) any later version.                                         |
|                                                                             |
| VM2D is distributed in the hope that it will be useful, but WITHOUT         |
| ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       |
| FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License       |
| for more details.                                                           |
|                                                                             |
| You should have received a copy of the GNU General Public License           |
| along with VM2D.  If not, see <http://www.gnu.org/licenses/>.               |
\*---------------------------------------------------------------------------*/
 
 
#include <algorithm>
 
#include "Mechanics2DRigidGivenLaw.h"
 
#include "Airfoil2D.h"
#include "Boundary2D.h"
 
 
#include "MeasureVP2D.h"
#include "StreamParser.h"
#include "Velocity2D.h"
#include "Wake2D.h"
#include "World2D.h"
#include "Gmres2D.h"
 
using namespace VM2D;
 
MechanicsRigidGivenLaw::MechanicsRigidGivenLaw(const World2D& W_, size_t numberInPassport_)
    : Mechanics(W_, numberInPassport_, true, false)
{
    MechanicsRigidGivenLaw::ReadSpecificParametersFromDictionary();
    Initialize(VelocityOfCenterOfMass(0.0), PositionOfCenterOfMass(0.0), AngularVelocity(0.0), RotationAngle(0.0));
};
 
 
//Вычисление гидродинамической силы, действующей на профиль
void MechanicsRigidGivenLaw::GetHydroDynamForce()
{
    W.getTimers().start("Force");
 
    const double& dt = W.getPassport().timeDiscretizationProperties.dt;
 
    hydroDynamForce = { 0.0, 0.0 };
    hydroDynamMoment = 0.0;
 
    viscousForce = { 0.0, 0.0 };
    viscousMoment = 0.0;
 
    Point2D hDFGam = { 0.0, 0.0 };      //гидродинамические силы, обусловленные присоед.завихренностью
    Point2D hDFdelta = { 0.0, 0.0 };    //гидродинамические силы, обусловленные приростом завихренности
    Point2D hDFQ = { 0.0, 0.0 };        //гидродинамические силы, обусловленные присоед.источниками
    
    double hDMGam = 0.0;            //гидродинамический момент, обусловленный присоед.завихренностью
    double hDMdelta = 0.0;          //гидродинамический момент, обусловленный приростом завихренности
    double hDMQ = 0.0;          //гидродинамический момент, обусловленный присоед.источниками
 
        
    for (size_t i = 0; i < afl.getNumberOfPanels(); ++i)
    {
        /*Point2D*/ Vcm = VeloOfAirfoilRcm(W.getCurrentStep() * dt);
        /*Point2D*/ VcmOld = VeloOfAirfoilRcm((std::max(W.getCurrentStep(), (size_t)1) - 1) * dt);
        
        /*double*/ Wcm = AngularVelocityOfAirfoil(W.getCurrentStep() * dt);
        /*double*/ WcmOld = AngularVelocityOfAirfoil((std::max(W.getCurrentStep(), (size_t)1) - 1) * dt);
 
        Point2D rK = 0.5 * (afl.getR(i + 1) + afl.getR(i));
        
        Point2D dr = rK - afl.rcm;
 
        Point2D velK = { -Wcm * dr[1], Wcm * dr[0] };
        velK += Vcm;
 
        double gAtt = Wcm * (dr ^ afl.tau[i]) + (Vcm & afl.tau[i]);
        double gAttOld = WcmOld * (dr ^ afl.tau[i]) + (VcmOld & afl.tau[i]);
        double deltaGAtt = gAtt - gAttOld;
 
        double qAtt = Wcm * (dr ^ afl.nrm[i]) + (Vcm & afl.nrm[i]);
    
        double deltaK = boundary.sheets.freeVortexSheet(i, 0) * afl.len[i] - afl.gammaThrough[i] + deltaGAtt * afl.len[i];
        
        /*1*/
        hDFdelta += deltaK * Point2D({ -dr[1], dr[0] });
        hDMdelta += 0.5 * deltaK * dr.length2();
 
        /*2*/
        hDFGam += 0.25 * (afl.getV(i) + afl.getV(i+1)).kcross() * gAtt * afl.len[i];
        hDMGam += 0.25 * dr ^ (afl.getV(i) + afl.getV(i + 1)).kcross() * gAtt * afl.len[i];
 
        /*3*/
        hDFQ -= 0.25 * (afl.getV(i) + afl.getV(i + 1)) * qAtt * afl.len[i];
        hDMQ -= 0.25 * dr ^ (afl.getV(i) + afl.getV(i + 1)) * qAtt * afl.len[i];        
    }
 
    const double rho = W.getPassport().physicalProperties.rho;
 
    hydroDynamForce = rho * (hDFGam + hDFdelta * (1.0 / dt) + hDFQ);
    hydroDynamMoment = rho * (hDMGam + hDMdelta / dt + hDMQ);
 
    if ((W.getPassport().physicalProperties.nu > 0.0) && (W.getCurrentStep() > 0))
        for (size_t i = 0; i < afl.getNumberOfPanels(); ++i)
        {
            Point2D rK = 0.5 * (afl.getR(i + 1) + afl.getR(i)) - afl.rcm;
            viscousForce += rho * (afl.viscousStress[i] * afl.tau[i]);
            viscousMoment += rho * (afl.viscousStress[i] * afl.tau[i]) & rK;
        }
 
    W.getTimers().stop("Force");
}// GetHydroDynamForce()
 
 
// Вычисление скорости центра масс
Point2D MechanicsRigidGivenLaw::VeloOfAirfoilRcm(double currTime)
{
    if (W.getPassport().physicalProperties.typeAccel.second == 3)
    {
        switch ((int)(W.getPassport().physicalProperties.timeAccel))
        {
        case 0:
            return { -1.0, 0.0 };
        case 1:
            return { 0.0, -1.0 };
        case 2:
            return { 0.0, 0.0 };
        default:
            return { 0.0, 0.0 };
        }
    }
 
    return VelocityOfCenterOfMass(currTime);
}//VeloOfAirfoilRcm(...)
 
 
// Вычисление положения центра масс
Point2D MechanicsRigidGivenLaw::PositionOfAirfoilRcm(double currTime)
{   
    return PositionOfCenterOfMass(currTime);
}//PositionOfAirfoilRcm(...)
 
// Вычисление угловой скорости профиля вокруг центра масс
double MechanicsRigidGivenLaw::AngularVelocityOfAirfoil(double currTime)
{
    if (W.getPassport().physicalProperties.typeAccel.second == 3)
    {
        switch ((int)(W.getPassport().physicalProperties.timeAccel))
        {
        case 0:
            return 0.0;
        case 1:
            return 0.0;
        case 2:
            return -1.0;
        default:
            return 0.0;
        }
    }
 
    return AngularVelocity(currTime);
}//AngularVelocityOfAirfoil(...)
 
 //Вычисление угла поворота профиля
double MechanicsRigidGivenLaw::AngleOfAirfoil(double currTime)
{
    return RotationAngle(currTime);
}//AngleOfAirfoil(...)
 
 
 
 
 
 
 
// Вычисление скоростей начал панелей
void MechanicsRigidGivenLaw::VeloOfAirfoilPanels(double currTime)
{
    Point2D veloRcm = VeloOfAirfoilRcm(currTime);
    double angVelo = AngularVelocityOfAirfoil(currTime);
 
    std::vector<Point2D> vel(afl.getNumberOfPanels(), { 0.0, 0.0 });
    for (size_t i = 0; i < afl.getNumberOfPanels(); ++i)
    {
        vel[i][0] = (afl.getR(i) - afl.rcm).kcross()[0] * angVelo;
        vel[i][1] = (afl.getR(i) - afl.rcm).kcross()[1] * angVelo;
        vel[i] += veloRcm;
    }
    
    afl.setV(vel);
 
    circulation = 2.0 * afl.area * angVelo;
    circulationOld = 2.0 * afl.area * WcmOld;
}//VeloOfAirfoilPanels(...)
 
void MechanicsRigidGivenLaw::Move()
{
    double t = W.getCurrentTime();
    double dt = W.getPassport().timeDiscretizationProperties.dt;
 
    //Point2D airfoilVelo = VeloOfAirfoilRcm(t);
    Point2D aflRcmOld = afl.rcm;
    double aflPhiOld = afl.phiAfl;
    afl.Move(PositionOfAirfoilRcm(t + dt) - aflRcmOld);
    afl.Rotate(AngleOfAirfoil(t + dt) - aflPhiOld);
 
    Vcm = VeloOfAirfoilRcm(t + dt);
    Phi = AngleOfAirfoil(t + dt);
    Wcm = AngularVelocityOfAirfoil(t + dt);
}//Move()
 
 
void MechanicsRigidGivenLaw::ReadSpecificParametersFromDictionary()
{
    /*
    mechParamsParser->get("timeAccel", timeAccel);
    W.getInfo('i') << "time of accelerated movement: " << "timeAccel = " << timeAccel << std::endl;
 
    mechParamsParser->get("initPosition", initPosition);
    W.getInfo('i') << "initial position: " << "initPosition = " << initPosition << std::endl;
 
    //mechParamsParser->get("targetVelocity", targetVelocity);
    //W.getInfo('i') << "target Velocity: " << "targetVelocity = " << targetVelocity << std::endl;
 
    mechParamsParser->get("targetAmplitude", targetAmplitude);
    W.getInfo('i') << "target Amplitude: " << "targetAmplitude = " << targetAmplitude << std::endl;
    */
 
   /*
    mechParamsParser->get("Comega", Comega);
    
    W.getInfo('i') << "frequency " << "Comega = " << Comega << std::endl;
    
    mechParamsParser->get("CA", CA);
        
    W.getInfo('i') << "Amplitude CA = " << CA << std::endl;
    */
}//ReadSpecificParametersFromDictionary()