C/C++ Runge-Kutta 방법으로 알아보는 심해 압력

여기서는 바다속의 수압이 깊이에 따라 어떻게 달라지는지 수치해석으로 알아봅시다. 정역학적 평형 (hydrostatic equilibrium)과 물의 상태 방정식 (equation of state)을 결부지어서 미분방정식을 세우고, 이를 Runge-Kutta 방법을 통해 풀어보겠습니다.

 

반응형

 

유체의 정역학적 평형에 대해서 간단히 말하자면, 위치에 따라 달라지는 압력으로 인해 생기는 힘이 외력과 상쇄가 되었을 때 알짜힘이 0이 되어 유체의 운동상태가 일정하게 유지되는 것을 의미하는데요. 이러한 정역학적 평형 상태를 정의하기 위한 미분방정식을 세우는 방법은 다음 포스팅에 더 자세하게 소개되어 있습니다.

 

물리학 상식 : 유체의 정역학적 평형

 

물질의 특징을 이해하기 위해 필요한 요소 중 하나인 상태방정식은 밀도와 압력 간의 상관관계를 나타낸 것입니다. 이들을 수학적으로 정의하는 방법에 대한 더 자세한 내용은 다음 포스팅에 소개되어 있습니다.

 

물리학 상식 : 압력, 밀도와 상태 방정식

 

정역학적 평형상태와 물의 상태 방정식을 결합하면, 수압을 깊이에 대한 함수로 구하는 데 필요한 미분방정식이 등장하게 됩니다. 해수면에서의 압력이 1기압이라는 초기조건과 더불어서 미분방정식을 풀면 되는 것입니다.

 

 

본격적인 계산에 들어가기에 앞서 미리 짚고 넘어갈 점이 있다면, 이 문제는 근사적으로 풀 수 있는 방법이 있다는 것입니다. 단위 질량당 물의 부피가 항상 일정하다고 가정하고 미분방정식을 푸는 것인데요. 그러면 수심이 10m 깊어질때마다 수압이 1기압만큼 증가한다는 결론을 얻을 수 있습니다. 지구상에서 가장 깊은 곳이라고 알려진 마리아나 해구의 수심이 10km 가량인데, 여기서는 수압이 1000기압에 달하게 되겠죠. 인간의 입장에서 보면 극한의 환경이고, 심해 탐사가 어려운 이유 중 하나이기도 합니다.

 

결론적인 얘기를 먼저 하자면, 물이 일정한 밀도를 가진다고 가정하는 것은 매우 타당한 근사법입니다. 상태 방정식과 정역학적 평형을 통한 정공법으로 계산한 결과와 비교해보면 그 차이는 5%보다 작습니다. 이런 얘기를 학술 저널같은데 들이밀면 바로 리젝을 먹겠지만, 여기는 제 개인 블로그니까 얘기를 계속하기로 하겠습니다.

 

C++ 소스 코드입니다.

 

ODESolveRK.h [다운로드]

 

ODESolveRK.cpp [다운로드]

 

초기조건이 주어진 연립 상미분 방정식을 4차 Runge-Kutta 방법으로 풀기 위한 ODESolveRK 클래스를 정의한 헤더 파일과 소스 파일입니다. 함수 포인터를 통해서 미분방정식의 형태를 정의하고, init 함수를 통해 초기화할 수 있습니다. 그 다음 next 함수를 반복적으로 호출하면서, 정의역으로 들어가는 변수의 값이 달라질 때 함수들의 값이 어떻게 변하는지를 알아내는 구조입니다.

 

여기에 소개된 ODESolveRK 클래스를 이용해서 미분방정식을 풀기 위한 더 자세한 방법은 다음 포스팅에 진자운동의 예시를 통해 소개되어 있습니다.

 

C/C++ Runge-Kutta 방법으로 알아보는 진자운동

 

test1_press_sea_RK4.cpp [다운로드]

 

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<vector>
#include"ODESolveRK.h"

/* length (z) is multiplied by g/(p_0 V_0),
 * while pressure (p) is divided by 1 atm.
 * By doing that,
 * we can deal with dimensionless quantities. */

const double z_ini = 0.;  // initial depth
const double z_fin = -1500.;  // final depth
const double delta_z = -0.02;  // size of timestep in RK
const double p_ini = 1.;  // initial pressure at the sea level

const double coeff_a = 0.315 / log(10.);
const double coeff_b = 2750.;

double vol_specific(double p);

// derivative of pressure
double dp_dz(double z, double *y);

int main(int argc, char *argv[]) {
    /* array of function pointer for derivative
     * of pressure with respect to depth */
    func_oderk_deriv *ptr_in_func;
    ptr_in_func =
        (func_oderk_deriv *)malloc(sizeof(func_oderk_deriv) * 1);
    ptr_in_func[0] = &dp_dz;

    /* arrays to store the history
     * of pressure */
    int nbin_z =
        (int)floor(fabs(z_fin - z_ini) / fabs(delta_z) + 1e-8);
    double *z_history;
    std::vector<double> *p_history;
    z_history = new double [nbin_z + 1];
    p_history = new std::vector<double> [nbin_z + 1];

    z_history[0] = z_ini;
    p_history[0].clear();
    p_history[0].push_back(p_ini);

    ODESolveRK ode_solver;
    // initialize the ODE Runge-Kutta solver
    bool initialized =
        ode_solver.init(1, z_history[0], p_history[0], ptr_in_func);

    for (int iz = 1; iz <= nbin_z; iz++) {
        bool go_next = ode_solver.next_RK4(delta_z);
        if (!go_next) {
            break;
        }

        ode_solver.get_system_current(z_history[iz], p_history[iz]);
        fprintf(stderr, "  %e  %e\n",
            z_history[iz], p_history[iz].at(0));
    }

    FILE *fout;
    char filename_out[200];
    strcpy(filename_out, "pressure_sea.dat");
    fout = fopen(filename_out, "w");
    fprintf(fout, "# coeff_a/V_0 = %e\n", coeff_a);
    fprintf(fout, "# coeff_b/p_0 = %e\n", coeff_b);
    fprintf(fout,
        "# [(g z)/(p_0 V_0)]  [p/p_0 numerical]  [p/p_0 approx]\n");

    for (int iz = 0; iz <= nbin_z; iz++) {
        double p_approx = 1. - z_history[iz];
        fprintf(fout, "  %e  %e  %e\n",
            z_history[iz], p_history[iz].at(0), p_approx);
    }

    fclose(fout);

    delete [] z_history;
    delete [] p_history;
    free(ptr_in_func);

    return 0;
}

double vol_specific(double p) {
    return 1. - coeff_a * log((p + coeff_b) / (1. + coeff_b));
}

double dp_dz(double z, double *p) {
    return -1. / vol_specific(p[0]);
}

 

main 함수가 포함된 소스파일입니다. 여기서 ODESolveRK 클래스 객체를 선언하고, 정역학적 평형상태를 나타내는 미분방정식의 형태를 알려주게 됩니다. 프로그램을 컴파일하고 실행시켜 보면 깊이에 따른 압력을 텍스트 파일의 형태로 얻을 수 있습니다.

 

• 컴파일
  g++ ODESolveRK.cpp -c
  g++ test1_press_sea_RK4.cpp -c
• 링크
  g++ test1_press_sea_RK4.o ODESolveRK.o -lm -o [실행파일 이름]

 

 

결과를 살펴보면, 일정한 밀도를 가정한 근사법에 비해서 차이가 매우 작습니다. 물은 압축률이 낮아서 결과적으로 특종이 아닌 낙종이 되어버렸습니다만, 다른 종류의 액체에 대해서도 비슷한 방법으로 계산이 가능합니다. 압력에 따라 밀도가 크게 달라지는 액체라면 다른 결과가 나올 수도 있겠죠.