편미분방정식 - 2차원 열 방정식 유한차분
이 장에서 배우는 것
앞 장에서 시간에 따라 변하는 상태를 상미분방정식으로 계산했다. 금속판의 온도는 시간뿐 아니라 위치에 따라서도 달라진다. 판의 가운데를 가열하면 열은 주변으로 퍼지고, 가장자리가 차갑게 유지되면 일부 열은 판 밖으로 빠져나간다. 이 현상을 계산하려면 시간 간격과 함께 공간 간격을 정해야 한다.
이 장에서는 편미분방정식(partial differential equation)을 격자 위의 계산으로 바꾸는 과정을 다룬다. 작은 판의 온도를 명시적으로 갱신하는 프로그램을 만들고, 계산이 안정되기 위한 조건을 확인한다. 암시적 방법은 계산 구조와 선택 기준까지만 살펴본다.
- 2차원 열 방정식의 공간 미분을 유한차분으로 나타낸다.
- 격자점의 인덱스와 실제 좌표를 연결하고 경계 조건을 적용한다.
- 명시적 갱신식의 계수와 시간 간격의 안정 조건을 계산한다.
- 이전 온도와 다음 온도를 구분하여 판의 온도 분포를 구한다.
- 명시적 방법과 암시적 방법의 계산 비용 및 제약을 비교한다.
문제 상황
센서가 달린 얇은 금속판의 가운데를 짧게 가열했다고 하자. 가열이 끝난 순간을 계산의 시작으로 잡는다. 이후 내부 발열은 없고, 판의 네 가장자리는 냉각 장치에 의해 20도로 유지된다. 알고 싶은 것은 가운데 온도가 얼마나 내려가는지, 주변이 어떤 순서로 따뜻해지는지다.
실제 판에서는 두께 방향의 온도 차이, 표면에서 공기로 빠지는 열, 재료 특성의 온도 의존성도 고려할 수 있다. 여기서는 두께 방향으로 온도가 같고, 판 내부의 열확산율이 일정하며, 열은 판의 두 방향으로만 이동한다고 가정한다. 이런 가정은 계산 결과가 무엇을 뜻하는지 정하는 모델의 일부다.
예제의 판은 가로와 세로가 각각 0.04 m다. 두 방향 모두 네 구간으로 나누므로 격자 간격은 0.01 m이고, 경계를 포함한 격자점은 5×5개다. 초기 온도는 모두 20도이며 가운데 격자점만 120도로 둔다. 한 점의 높은 온도는 작은 가열 영역을 거칠게 표현한 초기 조건이다. 실제 가열 영역을 더 정확하게 표현하려면 공간 격자를 더 촘촘하게 해야 한다.
출력은 두 번의 시간 갱신으로 제한한다. 작은 격자와 짧은 계산은 가운데에서 주변으로 온도가 전달되는 과정을 손으로 검산하기에 알맞다. 긴 시간의 분포를 얻는 일은 같은 갱신식을 더 많이 반복하는 작업이다.
격자와 경계 조건으로 문제를 정한다
온도를 T(x, y, t), 열확산율을 α라고 쓰면 내부 발열이 없는 2차원 열 방정식은 다음과 같다. 열확산율의 단위는 m²/s이며, 온도 차이가 공간적으로 고르게 퍼지는 속도를 결정한다.
∂T/∂t = α (∂²T/∂x² + ∂²T/∂y²)온도를 섭씨로 표현해도 이 식을 사용할 수 있다. 여기서는 온도 자체가 아니라 온도 차이와 미분이 계산에 들어가기 때문이다. 다만 다른 물리 법칙과 결합하여 절대 온도가 필요한 항을 추가한다면 온도 단위를 다시 검토해야 한다.
구간 수와 격자점 수를 구분한다
x 방향의 구간 수를 nx, y 방향의 구간 수를 ny라고 하자. 격자점은 각각 nx+1개와 ny+1개다. 배열을 temp(0:nx, 0:ny)로 선언하면 좌표의 시작과 인덱스의 시작을 함께 0으로 둘 수 있다.
dx = length_x / real(nx, real64)
dy = length_y / real(ny, real64)
x_i = i * dx
y_j = j * dytemp(i,j)는 좌표 (i dx, j dy)의 온도다. 내부점은 i=1부터 nx-1까지, j=1부터 ny-1까지다. i=0, i=nx, j=0, j=ny에 놓인 점들은 경계점이다. nx가 4라는 것은 점이 네 개라는 뜻이 아니라 구간이 네 개라는 뜻이다.
첫 번째 배열 인덱스를 x 방향으로 정하는 것은 프로그램의 약속이다. 배열 출력에서는 y가 큰 행부터 내려오도록 작성하여 위쪽이 판의 위쪽처럼 보이게 한다. 계산할 때 사용하는 인덱스 순서와 화면에 표시하는 행 순서는 구분해서 생각해야 한다.
초기 조건과 경계 조건은 역할이 다르다
초기 조건은 시작 시각의 전체 온도 분포를 정한다. 경계 조건은 계산하는 동안 가장자리에서 어떤 물리 상황을 유지할지 정한다. 이 예제에서는 가운데만 높은 초기 분포를 주고, 모든 시간에 가장자리를 20도로 고정한다.
| 종류 | 정하는 양 | 물리 상황 |
|---|---|---|
| 디리클레 조건 | 경계의 온도 | 가장자리를 일정 온도로 유지한다 |
| 노이만 조건 | 경계의 법선 방향 온도 미분 | 열유속을 지정하며, 미분이 0이면 단열을 나타낸다 |
| 로빈 조건 | 온도와 법선 방향 미분의 관계 | 주변 공기와의 열교환을 모델링한다 |
완성 코드에서는 디리클레 조건만 구현한다. 경계 온도를 고정하면 냉각 장치가 열을 받거나 공급할 수 있으므로 판 내부의 열량이 일정할 이유가 없다. 특히 내부가 경계보다 뜨거운 현재 문제에서는 시간이 지나면서 열이 경계를 통해 빠져나간다.
노이만 조건을 사용하려면 경계 근처에서 미분 조건을 만족하도록 별도의 식을 만들어야 한다. 내부점 갱신식을 가장자리까지 그대로 확장하는 것으로 경계 조건이 구현되지는 않는다. 어떤 조건을 선택했는지 먼저 정하고, 그 조건에 맞는 경계 계산을 작성해야 한다.
명시적 방법과 안정 조건
공간 미분을 이웃 온도의 차이로 바꾼다
유한차분법(finite difference method)은 미분을 유한한 간격에서 얻은 값의 차이로 근사한다. x 방향의 두 번째 미분은 현재 점과 좌우 이웃을 사용하고, y 방향의 두 번째 미분은 현재 점과 아래위 이웃을 사용한다.
∂²T/∂x² ≈ (T(i+1,j) - 2T(i,j) + T(i-1,j)) / dx²
∂²T/∂y² ≈ (T(i,j+1) - 2T(i,j) + T(i,j-1)) / dy²시간 미분에는 앞 장의 오일러 방법처럼 현재 시각에서 다음 시각으로 나아가는 차분을 적용한다. 시간 간격을 dt로 두면 다음 온도는 현재 온도와 현재 시각의 공간 차분만으로 계산할 수 있다. 이처럼 다음 값을 이미 알려진 값들로 구하는 방법이 명시적 방법이다.
rx = α dt / dx²
ry = α dt / dy²
T_new(i,j) = T_old(i,j)
+ rx * (T_old(i+1,j) - 2T_old(i,j) + T_old(i-1,j))
+ ry * (T_old(i,j+1) - 2T_old(i,j) + T_old(i,j-1))rx와 ry는 단위가 없는 수다. 같은 물성과 공간 격자에서 dt가 커지면 두 계수도 커진다. 같은 dt에서 공간 간격을 절반으로 줄이면 해당 방향의 계수는 네 배가 된다. 공간 해상도를 높이는 선택은 시간 간격의 제한에도 영향을 준다.
충분히 매끄러운 해에 대해 이 방법은 시간 방향으로 1차, 각 공간 방향으로 2차 정확도를 가진다. 이는 오차의 크기가 대체로 dt, dx², dy²의 크기에 따라 줄어드는 근사 구조라는 뜻이다. 불연속에 가까운 초기 분포나 성긴 격자에서는 이 차수만으로 실제 오차를 판단하기 어렵다.
계수의 합으로 안정 조건을 확인한다
갱신식을 다시 묶으면 현재 점의 계수는 1-2rx-2ry이고, 좌우 이웃의 계수는 각각 rx, 아래위 이웃의 계수는 각각 ry다. 열확산율과 시간 간격이 양수일 때 모든 계수가 음수가 아니려면 다음 조건이 필요하다.
rx + ry ≤ 1/2
dt ≤ 1 / (2α(1/dx² + 1/dy²))
dx = dy = h인 경우:
dt ≤ h² / (4α)이 조건은 여기서 사용하는 일정한 열확산율과 직사각형 균일 격자의 명시적 차분에 대한 안정 조건이다. 계수가 음수가 아니고 합이 1이면 새 온도는 다섯 온도의 가중 평균이다. 따라서 이전 분포의 최솟값과 최댓값을 벗어나지 않는다. 현재 예제처럼 일정한 경계 온도까지 이전 배열에 포함하면 온도가 초기 범위 밖으로 튀는 현상을 막을 수 있다.
안정 조건을 만족한다는 사실은 결과의 오차가 충분히 작다는 뜻은 아니다. 안정성은 오차가 계산 과정에서 부적절하게 증폭되는 문제와 관련되고, 정확도는 선택한 격자와 시간 간격이 원하는 해를 얼마나 잘 근사하는지와 관련된다. 안정 조건을 통과한 뒤에도 간격을 줄인 계산과 비교해야 한다.
예제에서는 α=0.0001 m²/s, dx=dy=0.01 m, dt=0.1 s다. 따라서 rx=ry=0.1이고 합은 0.2다. 허용되는 최대 시간 간격은 0.25 s이므로 선택한 시간 간격은 조건을 만족한다. 반면 dt=0.3 s로 바꾸면 합이 0.6이 되어 이 방법의 안정 조건을 벗어난다.
이전 배열과 다음 배열을 분리한다
명시적 갱신식에 들어가는 다섯 온도는 모두 같은 시각의 값이어야 한다. temp를 읽으면서 temp 자체에 새 값을 저장하면, 반복문 뒤쪽의 점은 이미 갱신된 이웃을 읽게 된다. 그러면 작성한 수식과 다른 계산이 되고 방문 순서에 따라 결과가 달라질 수 있다.
이 문제를 피하기 위해 temp에는 이전 온도를 보관하고 next_temp에는 다음 온도를 저장한다. 모든 내부점 계산이 끝난 뒤 temp=next_temp를 실행한다. 두 배열을 사용하는 구조는 각 시간 단계가 어디서 시작하고 끝나는지 드러내므로 계산을 검토하기도 쉽다.
암시적 방법은 무엇을 바꾸는가
명시적 방법은 간단하지만 작은 공간 간격에서 시간 간격이 크게 제한된다. 암시적 방법은 공간 미분을 다음 시각의 온도로 계산한다. 뒤로 오일러 방법을 적용하면 다음과 같은 식을 얻는다.
(1 + 2rx + 2ry) T_new(i,j)
- rx T_new(i-1,j) - rx T_new(i+1,j)
- ry T_new(i,j-1) - ry T_new(i,j+1)
= T_old(i,j)이번에는 이웃의 다음 온도도 미지수다. 한 점씩 식에 대입하는 것만으로 전체 해를 구할 수 없으며, 내부 격자점에 대한 식을 함께 풀어야 한다. 고정된 경계 온도는 알려진 값이므로 해당 항을 우변으로 옮긴다. 앞서 다룬 선형 연립방정식이 공간 전체의 온도 갱신에 나타나는 것이다.
2차원에서는 격자점 하나가 최대 네 이웃과 연결된다. 미지수를 한 줄로 번호 매기면 행렬에 여러 대각선이 생기므로 전체 문제를 하나의 삼중대각 연립방정식으로 바로 취급할 수는 없다. 격자가 커지면 연결된 항만 저장하거나 반복적으로 해를 구하는 방식이 계산 비용 면에서 중요해진다.
| 비교 항목 | 명시적 방법 | 뒤로 오일러 암시적 방법 |
|---|---|---|
| 사용하는 이웃 온도 | 이전 시각의 알려진 값 | 다음 시각의 미지수 |
| 한 단계의 핵심 계산 | 격자점마다 갱신식 평가 | 연립방정식 풀이 |
| 선형 확산 문제의 안정성 | 시간 간격 제한이 있다 | 명시적 방법과 같은 제한은 없다 |
| 시간 간격 선택 | 안정성과 정확도를 확인한다 | 정확도와 풀이 오차를 확인한다 |
이 선형 확산 문제의 뒤로 오일러 방법은 큰 시간 간격에서도 안정적이지만, 빠른 초기 변화를 거친 시간 간격으로 정확하게 표현할 수 있다는 뜻은 아니다. 연립방정식을 근사적으로 푼다면 그 풀이 오차도 관리해야 한다. 이 장의 완성 코드는 계산 흐름을 직접 확인하기 쉬운 명시적 방법으로 작성한다.
완성 코드
모든 코드는 main.f90 한 파일에 넣는다. 입력은 소스의 상수로 정하며 표준 입력은 읽지 않는다. 각 단계에서 다음 배열 전체를 경계 온도로 채우고 내부점만 갱신하여 네 가장자리의 온도를 유지한다.
main.f90
program plate_heat
use iso_fortran_env, only : real64
implicit none
integer, parameter :: nx = 4, ny = 4
integer, parameter :: nsteps = 2
real(real64), parameter :: length_x = 0.04_real64
real(real64), parameter :: length_y = 0.04_real64
real(real64), parameter :: alpha = 0.0001_real64
real(real64), parameter :: dt = 0.1_real64
real(real64), parameter :: boundary_temp = 20.0_real64
real(real64), parameter :: hot_temp = 120.0_real64
real(real64) :: temp(0:nx, 0:ny)
real(real64) :: next_temp(0:nx, 0:ny)
real(real64) :: dx, dy, rx, ry, time, change
integer :: i, j, step
if (nx < 2 .or. ny < 2) then
error stop 'Each direction needs at least two intervals.'
end if
if (length_x <= 0.0_real64 .or. length_y <= 0.0_real64) then
error stop 'Plate lengths must be positive.'
end if
if (alpha <= 0.0_real64 .or. dt <= 0.0_real64) then
error stop 'Alpha and dt must be positive.'
end if
dx = length_x / real(nx, real64)
dy = length_y / real(ny, real64)
rx = alpha * dt / dx**2
ry = alpha * dt / dy**2
if (rx + ry > 0.5_real64) then
error stop 'Explicit stability condition is violated.'
end if
temp = boundary_temp
temp(nx/2, ny/2) = hot_temp
write (*, '("rx=",F6.3," ry=",F6.3)') rx, ry
write (*, '("initial center=",F8.2)') temp(nx/2, ny/2)
do step = 1, nsteps
next_temp = boundary_temp
do j = 1, ny - 1
do i = 1, nx - 1
next_temp(i,j) = temp(i,j) &
+ rx * (temp(i+1,j) - 2.0_real64 * temp(i,j) &
+ temp(i-1,j)) &
+ ry * (temp(i,j+1) - 2.0_real64 * temp(i,j) &
+ temp(i,j-1))
end do
end do
change = maxval(abs(next_temp - temp))
temp = next_temp
time = real(step, real64) * dt
write (*, '("step=",I2," time=",F6.2," center=",F8.2, &
&" change=",F8.2)') &
step, time, temp(nx/2, ny/2), change
end do
write (*, '(A)') 'final temperature (top to bottom):'
do j = ny, 0, -1
write (*, '(*(F8.2))') (temp(i,j), i = 0, nx)
end do
end program plate_heat줄별 해설
use iso_fortran_env, only : real64는 실수 종류를 일정하게 선택한다. 길이, 열확산율, 시간 간격, 온도에 모두 같은 종류를 사용하고 실수 리터럴에도 _real64를 붙인다. 식을 계산한 뒤 저장할 때만 정밀도를 높이는 방식으로는 중간 계산의 정밀도를 확보할 수 없다.
nx와 ny는 구간 수이고 nsteps는 시간 갱신 횟수다. 배열의 경계를 0:nx와 0:ny로 지정했으므로 실제 배열 크기는 5×5다. 이 작은 예제에서는 크기가 컴파일 시점에 정해지는 배열로 충분하다.
세 묶음의 검사문은 내부점이 존재하는지, 판 길이가 양수인지, 열확산율과 시간 간격이 양수인지 확인한다. 이 검사는 나눗셈 전에 실행한다. 길이가 0인 상태에서 먼저 격자 간격을 계산하고 나중에 검사하면 계산 과정에 이미 잘못된 값이 들어갈 수 있다.
real(nx, real64)와 real(ny, real64)는 구간 수를 선택한 실수 종류로 명시적으로 변환한다. 이어지는 두 줄은 rx와 ry를 계산한다. 안정 조건 검사는 두 방향의 계수를 합쳐서 수행한다. 한 방향의 계수만 확인해서는 2차원 확산의 제한을 판단할 수 없다.
temp = boundary_temp는 경계와 내부를 모두 20도로 초기화한다. 다음 줄은 가운데 점의 온도만 120도로 바꾼다. 현재는 구간 수가 짝수여서 nx/2와 ny/2가 기하학적 가운데를 가리킨다. 홀수 구간으로 바꾸면 중앙 좌표가 격자점 사이에 있으므로 초기 가열 영역의 정의도 함께 검토해야 한다.
시간 반복문의 첫 줄 next_temp = boundary_temp는 새 시각의 경계를 준비한다. 이후 내부점만 계산하므로 가장자리에 저장한 20도는 바뀌지 않는다. 내부 전체를 새로 덮어쓰므로 이전 단계의 내부 값이 next_temp에 남아 영향을 주는 일도 없다.
바깥 반복문은 j, 안쪽 반복문은 i를 변화시킨다. Fortran 배열에서는 첫 번째 인덱스가 연속으로 저장되므로 이 순서는 연속한 원소를 따라 접근한다. 정확한 결과를 결정하는 핵심은 반복 순서보다 읽는 배열과 쓰는 배열을 구분하는 데 있다.
갱신식의 오른쪽에는 temp만 있고 왼쪽에는 next_temp만 있다. 줄 끝의 &는 자유 형식 소스의 이어쓰기다. 두 방향의 두 번째 차분을 별도로 읽을 수 있도록 식을 나눠 적었다.
change는 전체 격자에서 한 단계 동안 변한 온도의 절댓값 중 최댓값이다. 경계는 변하지 않으므로 경계의 기여는 0이다. 이 값은 해의 정확한 오차가 아니라 연속한 두 수치 상태의 차이다. 정상 상태에 가까워지는 정도를 살펴보는 보조 지표로 사용할 수 있다.
temp = next_temp는 모든 내부점 계산이 끝난 뒤 수행한다. 시간은 dt를 반복해서 더하지 않고 단계 번호와 dt의 곱으로 계산한다. 이 방식도 실수 표현 오차는 있지만, 반복 덧셈에서 생기는 누적 효과를 줄이고 단계와 출력 시각의 관계를 명확하게 한다.
마지막 반복문은 j를 ny부터 0까지 감소시킨다. 한 행에서는 i가 증가하므로 왼쪽에서 오른쪽으로 출력한다. *(F8.2)는 필요한 만큼 폭 8, 소수 둘째 자리 형식을 반복한다. 배열 크기를 바꾸어도 같은 형식으로 행을 출력할 수 있다.
실행 결과
macOS 또는 Linux에서 GNU Fortran 16으로 다음과 같이 컴파일하고 실행한다. 추가한 실행 시 검사 옵션은 배열 경계와 같은 오류를 확인하는 데 사용한다.
gfortran -std=f2018 -Wall -fcheck=all main.f90 -o plate_heat
./plate_heat예상 출력은 다음과 같다. 시간의 단위는 초이고 온도와 change의 단위는 섭씨도다.
rx= 0.100 ry= 0.100
initial center= 120.00
step= 1 time= 0.10 center= 80.00 change= 40.00
step= 2 time= 0.20 center= 60.00 change= 20.00
final temperature (top to bottom):
20.00 20.00 20.00 20.00 20.00
20.00 22.00 32.00 22.00 20.00
20.00 32.00 60.00 32.00 20.00
20.00 22.00 32.00 22.00 20.00
20.00 20.00 20.00 20.00 20.00첫 단계에서 가운데 온도는 120+0.1(20-240+20)+0.1(20-240+20)=80도다. 가운데와 직접 이웃한 네 내부점은 각각 30도가 된다. 대각선의 내부점은 아직 20도다. 이 갱신식은 한 단계에서 직접 연결된 이웃의 정보만 사용한다.
두 번째 단계에서 가운데는 80+0.1(30-160+30)+0.1(30-160+30)=60도가 된다. 직접 이웃은 32도, 대각선 내부점은 22도가 된다. 출력의 좌우 대칭과 위아래 대칭은 초기 조건, 경계 조건, 두 방향의 간격이 대칭이라는 사실과 일치한다.
20도를 뺀 내부 온도를 합하면 초기에는 100, 첫 단계에도 100, 두 번째 단계에는 96이다. 첫 단계 이후 경계 옆 내부점이 따뜻해져서 다음 단계에 열이 경계로 빠져나가기 때문이다. 이 합은 현재 같은 면적 간격을 가진 격자에서 열의 감소를 살펴보는 간단한 지표이며, 열량 단위를 얻으려면 밀도, 비열, 두께와 면적 요소를 함께 고려해야 한다.
격자점 사이의 실제 온도는 이 출력만으로 직접 알 수 없다. 또한 성긴 격자의 한 점에 준 고온은 초기 가열 영역을 세밀하게 표현하지 못한다. 결과를 설계 판단에 사용하려면 같은 물리적 초기 분포를 유지하면서 공간과 시간 간격을 줄여 관심 온도가 얼마나 달라지는지 확인해야 한다.
실무에서 자주 틀리는 것
이전 온도를 계산 도중 덮어쓴다
다음 코드는 오른쪽에서 읽는 temp를 왼쪽에서 즉시 바꾼다. 다음 점을 계산할 때 새 값과 이전 값이 섞인다.
! Wrong: overwrite values still needed by neighboring points.
temp(i,j) = temp(i,j) &
+ rx * (temp(i+1,j) - 2.0_real64 * temp(i,j) + temp(i-1,j)) &
+ ry * (temp(i,j+1) - 2.0_real64 * temp(i,j) + temp(i,j-1))고친 코드는 다음 배열에 저장한다. 배열 복사는 내부 반복문 두 개가 모두 끝난 뒤에 둔다.
next_temp(i,j) = temp(i,j) &
+ rx * (temp(i+1,j) - 2.0_real64 * temp(i,j) + temp(i-1,j)) &
+ ry * (temp(i,j+1) - 2.0_real64 * temp(i,j) + temp(i,j-1))
! After both spatial loops:
temp = next_temp두 방향의 안정 조건을 따로 검사한다
각 계수가 0.5 이하라고 해서 합도 0.5 이하인 것은 아니다. rx=ry=0.3이면 아래의 잘못된 검사는 통과하지만 합은 0.6이다.
! Wrong for this two-dimensional scheme.
if (rx > 0.5_real64 .or. ry > 0.5_real64) then
error stop 'Unstable.'
end if양수 조건을 확인한 뒤 두 계수의 합을 검사한다.
if (rx + ry > 0.5_real64) then
error stop 'Explicit stability condition is violated.'
end if경계점에 내부점 수식을 적용한다
아래 반복문은 i=0에서 i-1을, i=nx에서 i+1을 참조한다. y 방향에도 같은 문제가 있다. 경계점은 내부점과 역할이 다르다.
! Wrong: neighbor references go outside the array.
do j = 0, ny
do i = 0, nx
next_temp(i,j) = temp(i,j) &
+ rx * (temp(i+1,j) - 2.0_real64 * temp(i,j) + temp(i-1,j)) &
+ ry * (temp(i,j+1) - 2.0_real64 * temp(i,j) + temp(i,j-1))
end do
end do경계를 먼저 지정하고 이웃이 모두 존재하는 내부점만 계산한다.
next_temp = boundary_temp
do j = 1, ny - 1
do i = 1, nx - 1
next_temp(i,j) = temp(i,j) &
+ rx * (temp(i+1,j) - 2.0_real64 * temp(i,j) + temp(i-1,j)) &
+ ry * (temp(i,j+1) - 2.0_real64 * temp(i,j) + temp(i,j-1))
end do
end do격자를 늘리고 시간 간격을 그대로 둔다
판 길이를 유지하면서 nx와 ny를 8로 바꾸면 두 간격은 절반이 된다. 기존 dt=0.1 s를 유지하면 rx=ry=0.4가 되어 안정 조건을 벗어난다.
! Wrong combination for the original lengths and alpha.
integer, parameter :: nx = 8, ny = 8
real(real64), parameter :: dt = 0.1_real64기존 계수를 유지하려면 시간 간격을 4분의 1로 줄인다. 같은 종료 시각까지 계산하려면 단계 수는 네 배로 늘린다.
integer, parameter :: nx = 8, ny = 8
integer, parameter :: nsteps = 8
real(real64), parameter :: dt = 0.025_real64이 변경만으로 두 계산의 초기 조건이 같은 물리 분포가 되지는 않는다. 한 격자점만 가열하면 격자 간격에 따라 가열 영역의 크기가 달라진다. 격자 수렴을 조사할 때는 좌표에 대한 동일한 초기 온도 함수를 각 격자에서 평가해야 한다.
한눈에 보기
| 항목 | 식 또는 코드 | 확인할 의미 |
|---|---|---|
| 격자점 수 | (nx+1)×(ny+1) | 구간 수에 경계점을 포함한다 |
| 공간 간격 | dx=Lx/nx, dy=Ly/ny | 실제 길이와 인덱스를 연결한다 |
| 확산 계수 | rx=αdt/dx², ry=αdt/dy² | 시간 및 공간 간격의 영향을 나타낸다 |
| 안정 조건 | rx+ry≤0.5 | 현재 명시적 차분에 적용한다 |
| 내부 범위 | 1:nx-1, 1:ny-1 | 네 이웃이 배열 안에 존재한다 |
| 상태 갱신 | temp=next_temp | 전체 내부 계산 이후 실행한다 |
| 고정 온도 경계 | next_temp=boundary_temp | 내부를 덮어쓰기 전에 준비한다 |
| 정확도 확인 | 공간·시간 간격을 줄여 비교한다 | 동일한 물리 문제를 유지한다 |
계산의 핵심은 격자, 경계, 시간 간격, 상태 배열의 역할을 함께 맞추는 데 있다. 다음 장에서는 현재의 격자점 계산을 바탕으로 여러 점의 계산을 동시에 수행하는 방법을 살펴본다. 이를 위해서도 같은 시각의 배열을 읽고 별도의 배열에 쓰는 구조가 중요하다.
연습 문제
- 완성 코드에서 dt를 0.2 s로 바꾸고 nsteps를 1로 바꾼다. rx와 ry, 가운데 온도, 직접 이웃 네 점의 온도를 계산하라. 안정 조건을 만족하는지도 설명하라.
- 가로 길이 0.04 m, 세로 길이 0.02 m, nx=ny=4, α=0.0001 m²/s인 판을 생각하라. dx와 dy를 구하고 명시적 방법이 허용하는 최대 dt를 계산하라.
- 초기 온도를 경계를 포함하여 모두 20도로 바꾼다. 여러 단계를 계산해도 온도가 유지되는 이유를 식으로 설명하고, 초기화 코드에서 바꿀 부분을 제시하라.
- 원래 완성 코드에서 세 번째 시간 단계의 가운데 온도와 내부 온도의 최댓값을 구하라. 가운데 온도가 내려가는 동안 일부 주변점의 온도는 올라갈 수 있는 이유를 설명하라.
정답과 해설
-
rx=ry=0.2이므로 합은 0.4다. 안정 조건을 만족한다. 가운데 온도는 120+0.2(20-240+20)+0.2(20-240+20)=40도다. 직접 이웃의 초기 온도는 20도이며, 네 이웃 중 가운데만 120도이므로 새 온도는 20+0.2×100=40도다. 대각선 내부점과 경계는 20도로 유지된다.
원래 코드의 두 단계와 이 계산은 모두 0.2 s에 도달하지만 온도 분포는 다르다. 시간 간격에 따른 근사 오차가 다르기 때문이다. 최종 시각이 같다는 사실만으로 결과가 같아지는 것은 아니다.
-
dx=0.01 m, dy=0.005 m다. 따라서 rx+ry=αdt(10000+40000)=5dt이고, 5dt≤0.5에서 dt≤0.1 s를 얻는다. 세로 간격이 더 작아 y 방향 계수가 시간 간격을 더 강하게 제한한다. 동일한 간격을 가정한 h²/(4α) 식을 그대로 사용하면 이 직사각형 격자의 조건을 잘못 계산할 수 있다.
-
가운데를 가열하는 줄을 삭제하여 다음 초기화만 남긴다.
temp = boundary_temp모든 점이 20도이므로 두 방향의 차분은 각각 20-2×20+20=0이다. 내부 온도는 변하지 않고, 경계도 같은 20도로 유지된다. 출력의 center는 모든 단계에서 20.00, change는 0.00이 된다. 이런 일정한 분포 검사는 갱신식과 경계 처리가 일관되는지 확인하는 간단한 방법이다.
-
두 번째 단계의 가운데는 60도이고 직접 이웃은 모두 32도다. 세 번째 가운데 온도는 60+0.1(32-120+32)+0.1(32-120+32)=48.8도다.
직접 이웃 하나에는 가운데 60도, 경계 20도, 대각선 내부점 두 개의 22도가 연결된다. 따라서 새 온도는 0.6×32+0.1×(60+20+22+22)=31.6도다. 대각선 내부점은 0.6×22+0.1×(32+32+20+20)=23.6도다. 내부 최댓값은 가운데의 48.8도다.
대각선 내부점은 22도에서 23.6도로 올라간다. 열이 고온부에서 저온부로 전달되므로 일부 점이 따뜻해지는 동안 전체 최댓값은 내려갈 수 있다. 모든 점의 온도가 매 단계 내려가야 한다는 판단은 확산의 공간적 전달을 놓친 것이다.