7.10. 사용자 정의 물리 모델 (UDD, User-Defined Dynamics)

Flow360은 임의 물리 모델 입력을 지원합니다. 본 예제에서는 Proportional-Integral 제어기를 모사하여 Aero-Structural Interaction(ASI) 계산을 수행합니다.

7.10.1. 양력 계수 제어를 위한 받음각 PI 제어기

본 예제에서 om6Wing 예제에 제어기를 추가합니다. 이 예제에서 솔버 설정 파일을 기준으로 해석에서 얻은 양력 계수, CL은 약 0.26입니다. 해석 목표는 이 예제에 제어기를 추가하여 CL의 주어진 목표 값(예를 들어, 0.4)에 대해 필요한 받음각, alphaAngle을 추정하는 것 입니다. 솔버 설정 파일 내에 이를 위한 PI 제어기 정의는 다음과 같습니다.

"userDefinedDynamics": [
    {
        "dynamicsName": "alphaController",
        "inputVars": [
            "CL"
        ],
        "constants": {
            "CLTarget": 0.4,
            "Kp": 0.2,
            "Ki": 0.002
        },
        "outputVars": {
            "alphaAngle": "if (pseudoStep > 500) state[0]; else alphaAngle;"
        },
        "stateVarsInitialValue": [
            "alphaAngle",
            "0.0"
        ],
        "updateLaw": [
            "if (pseudoStep > 500) state[0] + Kp * (CLTarget - CL) + Ki * state[1]; else state[0];",
            "if (pseudoStep > 500) state[1] + (CLTarget - CL); else state[1];"
        ],
        "inputBoundaryPatches": [
            "1"
        ]
    }
]

사용자 정의 물리 모델의 매개변수에 대한 자세한 설명은 다음과 같습니다.

dynamicsName: 이 항목은 사용자가 지정한 물리 모델 이름입니다. 사용자 정의 물리 모델의 결과는 “udd_dynamicsName_v2.csv”에 저장됩니다.
constants: 여기에는 UDD와 관련된 모든 상수가 포함됩니다. 설정된 상수는 updateLaw 및 outputVars에서 사용됩니다. 본 예제에서 사용한 CLTarget, Kp 및 Ki는 각각 양력 계수의 목표 값, 제어기의 비례(Propotional) 이득 및 적분(Integral) 이득을 나타냅니다.
outputVars: 각 outputVars 는 변수 이름과 변수와 상태 변수 간의 관계로 구성됩니다. 본 예제에서 출력 변수는 alphaAngle이며, 처음 500 번의 pseudoSteps(즉, state[0]) 동안 초기 받음각(state[0])이 유지되다가, 이후 PI 제어기가 계산한 alphaAngle 결과가 해석의 alphaAngle 변수에 지정되도록 설정되어 있습니다.
stateVarsInitialValue: 이 제어기에는 두 개의 상태 변수가 있습니다. 첫 번째, state[0]는 제어기가 계산한 받음각이고, 두 번째, state[1]은 반복 계산 과정에서 도출된 오류(CLTarget-CL)의 합계입니다. state[0] 및 state[1]의 초기 값은 각각 alphaAngle 및 0로 설정됩니다.

제어 시스템에 대한 방정식은 다음과 같습니다.

$\alpha^n = \alpha^{n-1} + [K_p \cdot (C_{L,target}-C_L^n)+K_i \cdot e^n]$
$e^n = e^{n-1} +C_{L,target} - C_L^n$

여기서 상첨자 $n$ 은 pseudo step을 나타내고, alphaAngle ( $\alpha$ )는 state[0], 오차의 적분 ( $e$ )는 state [1]을 각각 나타냅니다.

updateLaw: 상태 변수에 대한 표현식(expression)이 여기에 지정됩니다. 첫 번째 및 두 번째 표현식은 각각 state[0] 및 state [1]에 대한 방정식에 해당됩니다. 조건문은 제어기가 처음 500 번의 pseudoSteps 동안 실행되지 않도록 강제합니다. 이렇게 하면 제어기에 의한 동적 변화가 반영되기 전에 유동장을 수렴시킬 수 있습니다.
inputBoundaryPatches: 사용자 정의 물리 모델에 사용되는 경계 조건이 여기에 지정됩니다. 본 예제에서 CL은 noSlipWall 경계에서 계산됩니다. 솔버 설정 파일에서 볼 수 있듯이 No-Slip 경계면의 지정된 경계 이름은 "1" 입니다.

다음 그림은 pseudoStep에 따른 양력 계수, $CL$ 과 받음각 alphaAngle에 대한 결과를 보여줍니다. 솔버 설정 파일에서 지정한 대로 받음각의 초기 값은 3.06입니다. 입력한 사용자 물리 모델에 의해 목표 양력 계수에 도달할 때 까지 alphaAngle은 자동으로 업데이트 됩니다. pseudoStep 2000에서 목표 양력 계수인 $CL=0.4$ 에 도달하였으며, 이 때 받음각, alphaAngle은 4.61로 산출 됩니다.

Fig. 7.10.1 Results of alphaController for om6Wing case study.

7.10.2. 구조-공기 역학적 하중을 사용한 동적 격자 회전

이번 예제에서는 평판이 공기역학적 하중과 회전 스프링 모멘트, 감쇠에 의한 영향을 동시에 받을 때, 힘의 균형 및 불균형에 의한 동적 회전 상태를 모사합니다.

아래 그림에 나타나 있듯이, 평판 모델은 원형 슬라이딩 인터페이스로 분리된 2개의 구역으로 분리되어 있습니다. 안쪽 구역은 회전 구역으로 설정되고, 바깥쪽 구역은 고정된 것으로 설정합니다. 이 방법은 회전익기에서 로터의 회전을 모사하는 방법과 동일합니다.

Fig. 7.10.2 Mesh of flat plate showing sliding interface for rotation.

아래 동영상은 평판의 flutter 거동을 보여줍니다.

Fig. 7.10.3 Flutter motion of the plate under aerodynamic and structural forces.

솔버 설정 파일(JSON) 내에서 사용자 정의 물리 모델(userDefinedDynamics)와 슬라이딩 인터페이스(slidingInterfaces)는 아래와 같이 설정됩니다.

{
    "slidingInterfaces": [
        {
            "stationaryPatches": [
                "farFieldBlock/rotIntf"
            ],
            "rotatingPatches": [
                "plateBlock/rotIntf"
            ],
            "axisOfRotation": [
                0,
                1,
                0
            ],
            "centerOfRotation": [
                0,
                0,
                0
            ],
            "volumeName": "plateBlock",
            "theta": [
                "0",
                "0",
                "0"
            ]
        }
    ],
    "userDefinedDynamics": [
        {
            "dynamicsName": "dynamicTheta",
            "inputVars": [
                "momentY"
            ],
            "constants": {
                "I": 0.443768309310345,
                "zeta": 4.0,
                "K": 0.0161227107,
                "omegaN": 0.190607889,
                "theta0": 0.0872664626
            },
            "outputVars": {
                "omegaDot": "state[0];",
                "omega": "state[1];",
                "theta": "state[2];"
            },
            "stateVarsInitialValue": [
                "-1.82621384e-02",
                "0.0",
                "1.39626340e-01"
            ],
            "updateLaw": [
                "if (pseudoStep == 0) (momentY - K * ( state[2] - theta0 ) - 2 * zeta * omegaN * I *state[1] ) / I; else state[0];",
                "if (pseudoStep == 0) state[1] + state[0] * timeStepSize; else state[1];",
                "if (pseudoStep == 0) state[2] + state[1] * timeStepSize; else state[2];"
            ],
            "inputBoundaryPatches": [
                "plateBlock/noSlipWall"
            ],
            "outputTargetName": "plateBlock"
        }
    ]
}

본 예제에서 평판에 대한 동적 거동은 다음 방정식으로 정의됩니다.

$\dot \omega = [\tau_y -K \cdot(\theta-\theta_0)-2\zeta \omega_N I \omega]/I$
$\omega - \omega_0 = \int_0^t \dot \omega (\acute t) d \acute t$
$\theta - \theta_0 = \int_0^t \omega (\acute t) d \acute t$

위 방정식은 상기 JSON 파일 내에 updateLaw 항목에 있는 3개의 인수에 해당합니다. 각 기호에 대한 정의는 아래와 같습니다.

Symbol	Description
$\theta$	평판의 회전 각도 [ $rad$ ]
$\omega$	평판의 회전 속도 [ $rad/sec$ ]
$\dot \omega$	평판의 회전 가속도 [ $rad/sec^2$ ]
$\tau_y$	평판에 작용하는 공기역학적 모멘트
$K$	구조 지지대와 평판 사이의 스프링 강성(spring stiffness)
$\zeta$	구조적 감쇠비(structural damping ratio)
$\omega_N$	구조적 고유 각 진동수(structural natural angular frequency)
$I$	구조적 관성 모멘트(structual moment of inertia)

inputVars에서 MomentY는 공기역학적 하중에 의해 부과된 회전 모멘트의 y 성분입니다. 모멘트는 "plateBlock/noSlipWall"의 momentCenter를 기준으로 계산됩니다.

제어하려는 회전 영역의 이름은 plateBlock이므로 "outputTargetName"은 plateBlock으로 설정합니다. 첫 번째, 두 번째 및 세 번째 상태 변수는 각각 $\dot \omega$ , $\omega$ 그리고 $\theta$ 에 해당됩니다.

7.10.3. 사용자 정의 힘 및 모멘트 세트

본 예제에서는 조종면(에일러론, 러더)이 있는 간단한 비행기 형상을 사용하여 해석이 수렴함에 따라 조종면의 토크 변화를 추적하기 위해 사용자 정의 힘과 모멘트를 계산하는 방법에 대해 다룹니다.

형상과 관련 격자를 생성하려면, 다음 링크에 포함된 입력 csm, 표면 격자 정의, 체적 격자 정의를 사용합니다. 격자 생성에 대한 자세한 방법은 튜토리얼1.4나 튜토리얼1.4를 참고하세요.

아래 이미지에서 확인할 수 있듯이 예제에서 사용한 비행기는 2개의 에일러론과 2개의 러더가 장착되어 있습니다. 4개의 조종면은 모두 각각의 회전 축을 중심으로 회전합니다. 본 예제에서는 각 조종면의 회전 축을 기준으로 토크를 모니터링 하는 방법을 설명합니다.

[ Note]

본 예제에서 사용한 격자 설정은 최소한의 격자만을 생성하도록 되어있습니다. 이로 인해 해석 케이스 실행 시 Flow360 솔버가 경고할 수 있으나 이를 무시하고 진행 합니다.

Fig. 7.10.4 Isometric view of a simple airplane with 4 control surfaces

체적 격자 준비가 끝났으면, 솔버 매개변수를 정의해야 합니다. 본 예제에서는 userDefinedDynamics 정의에 대해 집중적으로 설명합니다. 예제용 JSON 파일을 다운로드한 뒤, 내용을 확인합니다. 받음각은 10 $\degree$ 로 설정되어 있으며, 비행기의 속도는 50 m/sec로 설정되어 있습니다.

조종면의 회전 중심과 회전 축은 아래 표와 같습니다.

Name	rotation Center	rotation Axis
right Aileron	(5.7542, 7, 0)	(0, 1, 0)
left Aileron	(5.7542, -7, 0)	(0, -1, 0)
right Rudder	(12.01, 0.861, 0.861)	(0, 0.7071, 0.7071)
left Rudder	(12.01, -0.861, 0.861)	(0, -0.7071, 0.7071)

[ Note]

본 예제에서 geometry/MomentConter는 에일러론 회전 축 상에 임의의 X 및 Z 위치에 설정되었습니다. 이 방법은 에일러론 모멘트와 Flow360에서 Y 축에 대해 자동 계산된 모멘트를 비교하여 검증할 수 있습니다(에일러론 힌지 축이 Y 축에 정렬되어 있으므로, Flow360에서 자동 계산한 Y 축 모멘트와 동일한 값이 도출됨).

사용자 정의 변수를 통한 조종면 힌지 모멘트를 계산하면, 해석이 진행되는 동안 올바르게 정렬된 힌지 모멘트 값의 변화와 수렴성을 확인할 수 있습니다. 모니터링된 토크는 SI 단위의 힘(N)과 모멘트(Nm)로 변환할 수 있습니다.

입력 JSON 파일 내의 userDefinedDynamics(UDD) 항목에는 4개의 유사한 내용이 포함되어 있습니다. 각 조종면에는 고유한 userDefinedDynamics 인스턴스가 필요하며, 각 인스턴스에는 위 표에 나열된 조종면 이름과 회전 중심 및 회전 축 값이 각각 입력됩니다.

JSON 파일의 "userDefinedDynamics" 항목을 살펴보면 4개 조종면에 대한 각각의 인스턴스 항목은 6개의 변수를 포함하고 있습니다. 각 변수에 대한 설명은 다음과 같습니다.

Fig. 7.10.5 UDD section of the case configuration json file

7.10.3.1 constants

이름이 지정된 상수(Named constat)는 updateLaw에서 사용됩니다. "newCenterX, newCenterY, newCenterZ" 및 "newAxisX, newAxisY, newAxisZ"는 새로운 모멘트 참조 중심 및 회전 축의 좌표입니다. 이 값들은 각각의 러더와 에일러론의 힌지선을 정의합니다.

 1"constants": {
 2    "density_kgpm3": 1.225,
 3    "c_inf_mps": 340.2,
 4    "l_grid_unit": 1,
 5    "newCenterX": 5.7542,
 6    "newCenterY": 7,
 7    "newCenterZ": 0,
 8    "newAxisX": 0,
 9    "newAxisY": 1,
10    "newAxisZ": 0
11}

7.10.3.2 updateLaw

updateLaw 블록은 저장하고자 하는 출력 값을 정의하는데 사용하는 표현식(expression) 목록입니다. 솔버에서 출력되는 결과를 메트릭(metric) 값으로 변환하기 위하여 차원 값 $\rho_\infty C_\infty^2 L_{gridUnit}^3$ 을 계산합니다.

두 번째 줄의 방정식은 차원 값을 정의합니다. 세 번째에서 다섯 번 째 줄은 위에서 정의한 새로운 모멘트 참조 중심과 회전 축에 따라 새로운 모멘트를 정의합니다. 마지막 줄은 총 힌지 모멘트를 정의합니다. 이 값은 조종면을 회전시키기 위해 조종면의 작동 메커니즘이 극복해야할 토크를 의미합니다.

"updateLaw": [
    "density_kgpm3 * c_inf_mps * c_inf_mps * l_grid_unit * l_grid_unit * l_grid_unit",
    "(rotMomentX - ((newCenterY - momentCenterY) * forceZ - (newCenterZ - momentCenterZ) * forceY)) * state [0];",
    "(rotMomentY + ((newCenterX - momentCenterX) * forceZ - (newCenterZ - momentCenterZ) * forceX)) * state [0];",
    "(rotMomentZ - ((newCenterX - momentCenterX) * forceY - (newCenterY - momentCenterY) * forceX)) * state [0];",
    "state[1] * newAxisX + state[2] * newAxisY + state[3] * newAxisZ "
]

7.10.3.3 inputBoundaryPatches

격자 모델에서 계산을 적용하고자 하는 패치. 예를 들면:

"inputBoundaryPatches": [ "fluid/aileronRight" ]

7.10.3.4 visualizing the results

모든 UDD 출력은 results/udd+DYNAMICSNAME_v2.csv에 저장됩니다. 여기서 DYNAMICSNAME은 "dynamicName" 항목에 입력된 값을 가리킵니다.

[ Note]

클라우드 버전에서 Monitors 탭 항목에서 새로 정의된 모니터 결과를 시각화 할 수 있습니다. 힌지 토크의 차원 값은 state[4]에 저장되므로, Monitors 탭에서 해당 값을 체크하여 pseudoStep에 따른 결과를 확인할 수 있습니다.

Fig. 7.10.6 Dimensional hinge torque convergence history for the right Aileron

새롭게 정의된 모멘트 게산이 올바른지 확인하기 위하여, 경계면에 작용하는 전체 힘을 결과 토크 값과 비교합니다. 위에서 언급하였듯이 geometry/MomenCenter는 에일러론 힌지 선에 배치되어 있습니다. 두 값 모두 X=5.7542에 있으며, 에일러론 힌지 벡터 [(0, 1, 0)]이 Y 축과 정렬되어 있기 때문에 왼쪽과 오른쪽 에일러론의 $CMy$ 크기는 각각 왼쪽과 오른쪽 에일러론의 힌지 토크와 동일해야 합니다. 왼쪽 에일러론의 힌지 벡터는 [(0, -1, 0)]으로 설정되어 있으므로, $CMy$ 와 힌지 토크 사이의 부호가 반대로 나타납니다.

새롭게 정의된 오른쪽 에일러론의 힌지 토크 값과 Flow360 솔버에서 자동으로 산출되는 $CMy$ 값을 차원화한 값을 비교합니다.

right Aileron hinge torque	dimensional $CMy$
52.65 $N \cdot m$	-52.65 $N \cdot m$

[ Note]

클라우드 버전에서 아래 그림에 빨간색 원으로 표시된 DataView 아이콘을 클릭합니다. state[4] 열로 스크롤하여 마지막 행까지 내려간 뒤, 테이블 아래의 6번을 클릭하여(빨간색 표시) 테이블의 마지막 부분으로 이동합니다. 하이라이트된 -52.65 $Nm$ 값이 오른쪽 에일러론 힌지 토크에 대한 SI 단위 값입니다.

Fig. 7.10.7 Dimensional hinge torque table for the right Aileron.

Forces 탭으로 돌아간 뒤, Forces and Moments by surface 항목에서 최종 $CMy$ 값을 다음과 같이 확인할 수 있습니다.

$CMy = -5.73 \times 10^{-4}$

Flow360에서 자체적으로 산출되는 힘과 모멘트는 무차원화 된 값이므로, 다음 과정을 통해 SI 단위의 모멘트로 환산 할 수 있습니다.

$-5.73 \times 10^{-4} \cdot \frac{1}{2} \rho_\infty U_{ref}^2 A_{ref} L_{ref} [1]$
$= -5.73 \times 10^{-4} \cdot \frac{1}{2} \cdot 1.225 \frac{kg}{m^3} \cdot (50 \frac{m}{s})^2 \cdot 60m^2 \cdot 1m$
$=-52.65 N \cdot m$