검색
검색 팝업 닫기

Ex) Article Title, Author, Keywords

Article

Split Viewer

연구논문(Research Paper)

2024; 35(2): 49-60

Published online April 25, 2024 https://doi.org/10.3807/KJOP.2024.35.2.049

Copyright © Optical Society of Korea.

Number of Phase Screens Required for Simulation of a High-energy Laser Beam’s Propagation Experiencing Atmospheric Turbulence and Thermal Blooming

Seokyoung Yoon, Woohyeon Moon, Hoon Kim

대기 난류와 열적 블루밍을 겪는 고출력 레이저 빔의 대기 전파 시뮬레이션에 필요한 위상판 개수 분석

윤석영ㆍ문우현ㆍ김 훈

School of Electrical Engineering, Korea Advanced Institute of Science and Technology, Daejeon 34141, Korea

한국과학기술원 전기및전자공학부 ㉾ 34141 대전광역시 유성구 대학로 291

Correspondence to:hoonkim@kaist.ac.kr, ORCID: 0000-0001-7395-3695
Current affiliation: Republic of Korea Army, Gyeryong 32800, Korea

Received: December 26, 2023; Revised: January 23, 2024; Accepted: January 31, 2024

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

We analyze the number of phase screens required for the simulation of a high-energy laser beam’s propagation over an atmospheric channel. For high-energy lasers exceeding tens of kilowatts (kW) in power, the laser beam is mainly affected by atmospheric turbulence and thermal blooming. When using the split-step method to implement losses due to atmospheric absorption and scattering and distortion of the beam due to turbulence and thermal blooming, the number of phase screens is a critical factor in determining the accuracy and time required for the simulation. By comparing simulation results obtained using a large number of phase screens (e.g., 150 screens) under a wide range of atmospheric turbulence conditions, we provide new guidelines for the number of phase screens required for simulating the beam propagation of a high-power laser below 2.5 × 106 W/m2 (e.g., a 500-kW laser beam having a 50-cm diameter).

Keywords: High-energy laser, Laser weapon, Thermal blooming, Turbulence

OCIS codes: (010.1290) Atmospheric optics; (010.1330) Atmospheric turbulence

고출력 레이저(high-energy laser, HEL) 무기는 빛의 속도로 전파되며, 원하는 목표물에 정확히 조사되고, 미사일이나 포탄에 비교할 때 발사 비용이 저렴한 동시에 탄약의 수에 구애받지 않는다는 다양한 장점이 있다. 이러한 장점으로 인해 고출력 레이저 무기는 적의 전 영역에 걸친 군사적 위협에 대응할 수 있는 미래 핵심 전력으로서 전 세계가 주목하고 있다[1]. 그러나 레이저 빔은 조사 후 대기를 전파하는 중에 대기 매질의 굴절률 변화 및 대기와의 상호작용으로 인하여 빔의 세기와 품질 등이 변화하므로 무기의 성능이 감소할 수 있다. 따라서 레이저 무기의 성능을 극대화하기 위해서는 고출력 레이저 빔의 대기 전파 특성을 정확하게 이해하는 것이 필수이다.

레이저 빔은 대기 전파 중 흡수와 산란을 통해 단위 면적당 강도가 감소한다. 또한 대기에 존재하는 불균일한 공기 입자와 시시각각 변화하는 난류로 인하여 신틸레이션(scintillation), 빔 원더링(beam wandering)과 같은 선형적 효과가 발생하고, 그로 인해 목표물에 도달하기 전 에너지가 분산되거나 목표물 조준 오차가 발생하여 타격의 효율이 낮아진다[2]. 또한 레이저 출력이 킬로와트(kW) 급 이상으로 커지면, 강한 빛으로 인해 대기 물질이 열을 흡수하며 대기의 굴절률을 변화시키는 열적 블루밍(thermal blooming) 현상이 발생한다. 열적 블루밍 현상으로 인해 레이저 빔은 공간적으로 퍼지고, 전파 방향이 왜곡되어 빔의 조사가 부정확해진다[3].

레이저 빔이 대기를 전파하며 겪는 광학 현상을 이해하고, 이러한 현상이 레이저 빔 조사에 미치는 영향을 분석하기 위하여 통상적으로 컴퓨터 시뮬레이션이 활용된다[4,5]. 특히 시뮬레이션을 활용하면 실제 실험에서 조정하기 어려운 빔의 출력 세기, 난류의 세기, 전파 거리 등의 여러 조건을 손쉽게 변화시키면서 무기의 성능을 예측할 수 있다. 또한 시뮬레이션을 반복 실행함으로써 난류와 열적 블루밍의 무작위 특성이 빔에 끼치는 영향을 통계학적으로 분석할 수 있다[6]. 예를 들어 고출력 레이저 빔의 목표물에서의 형태, 강도, 조준점은 같은 난류 세기 조건에서도 난류의 분포 형태, 그리고 시간에 따라서 지속적으로 변한다. 따라서 정확한 성능 분석을 위해서 목표 거리에서 조준점에 모인 에너지의 양(power in the bucket, PIB)의 확률 질량 함수(probability mass function, PMF) 등을 활용하여 성능을 분석해야 한다. 하지만 통계학적 특성을 파악하기 위하여 시뮬레이션을 대량으로 반복 실행할 경우, 전체 시뮬레이션 시간이 반복 횟수에 비례하여 크게 증가한다. 따라서 시뮬레이션의 정확성을 희생하지 않으면서도 시뮬레이션 시간을 줄이는 것이 중요하다.

레이저 빔의 대기 전파는 확률적 파동 방정식(stochastic wave equation)으로 기술되며 이 방정식은 통상 split step 방법을 활용하여 계산된다[4]. Split step 방법은 빔의 전파 경로상에 특정 전파 거리에 해당하는 위상 왜곡 양을 가지는 위상판(phase screen)을 일정한 간격으로 배치하고, 위상판과 위상판 사이는 자유 공간으로 구성한다. 위상판은 대기의 난류나 열적 블루밍 특성에 맞추어 공간적으로 다른 위상을 가진 매우 얇은 판으로, 특정 거리의 대기를 전파할 때 발생하는 위상 변화를 모사하였다. 대기 난류를 전파하는 레이저 빔을 이와 같은 방법으로 계산할 때 고려해야 하는 사항은 이미 여러 문헌에서 보고되었는데, 예를 들면 난류 위상판 생성 간 고려되어야 할 와류(eddy)의 내외부 규모(inner and outer scale)에 따른 위상판의 크기, 에일리어싱(aliasing) 방지를 위한 위상판 사이의 거리 및 한 개의 위상판이 가져야 할 난류의 세기 정도에 대한 기준이 알려져 있다[7,8]. 열적 블루밍의 모사 또한 위상판의 형태로 전파 경로에 난류 위상판과 나란히 배치하는 방법이 연구된 바 있다[9]. 하지만 대기 난류 현상과 열적 블루밍을 함께 고려한 기존의 시뮬레이션 연구들은 두 효과를 정확히 모사하는 것에만 집중하고 있다. 그 결과 난류와 열적 블루밍을 모사하는 위상판의 개수가 많을수록 좋다는 것 이외에 시뮬레이션 시간을 줄이면서 두 현상을 효과적으로 모사하는 데 필요한 위상판 개수에 대한 분석은 이루어지지 않았다[9,10]. 위상판의 개수가 증가함에 비례해 시뮬레이션 내에서 자유공간 전파와 각 효과의 위상판 생성 횟수가 증가하는 만큼, 전체 시뮬레이션에 소요되는 시간은 선형적으로 증가한다.

가장 효과적으로 전체 시뮬레이션에 소요되는 시간을 줄이는 방법은 시뮬레이션 실행에 가장 많은 시간이 필요한 위상판 생성 및 자유 공간 전파의 횟수를 줄이는 것이다. 이는 사용하는 위상판의 개수를 줄임으로써 가능하다. 이때 핵심은 적은 개수의 위상판을 사용하더라도 빔이 전파되는 동안 난류와 열적 블루밍으로 인해 변화하는 특성이 많은 개수의 위상판을 사용하였을 때와 같아야 한다는 것이다. 따라서 본 논문에서는 열적 블루밍과 난류 현상을 모사하면서 동시에 변화되는 빔의 특성 또한 유지되는 최소 위상판의 개수를 제시하고, 시뮬레이션 시간을 크게 줄일 수 있는 방법을 보고한다. 비교를 위하여 충분히 많은 수의 위상판(150개)을 사용하여 난류와 열적 블루밍을 모사하고, 대조군으로 본 논문에서 제안하는 시뮬레이션은 전파 환경 및 레이저 전력에 맞춰 수량을 최소화한 위상판을 사용한다. 목표물에서의 빔의 성능은 조준점의 PIB를 통해 계산하며, 두 효과의 모사 여부를 확인하기 위하여 전파 경로 축상의 신틸레이션 지수의 평균 변화량을 측정한다. 시뮬레이션 비교 결과, 전파 거리가 5 km, 굴절률 구조 상수 함수가 10−14 m−2/3 이하, 왜곡 지수(distortion number)가 220 rad 이하일 때, 기존의 방법은 시뮬레이션 10,000회 반복에 18일이 소요된 반면, 제안된 방법은 4일이 소요되어 시뮬레이션 시간을 약 4.5배 단축할 수 있었다.

2.1. 확률적 파동 방정식과 split step 방법

대기를 통과하는 레이저 빔의 전파는 확률적 파동 방정식을 통해 모사할 수 있다. z축 방향으로 전파하는 전기장 E에 대한 확률적 파동 방정식은 근축 근사를 적용하여 식 (1)과 같이 표현된다[11].

2ikEz+2E+2k2nRE=0

여기에서 ▽2 = 2/∂x2+2/∂y2은 라플라스 연산자, k는 파수, 그리고 n(R)은 공간 위치 R에 대한 굴절률을 의미한다. 상기 방정식은 통상 split step 방법으로 계산되며, 대기 난류에 의한 공간적 굴절률 변화는 위상판으로 표현할 수 있다. 위상판은 특정 거리에 해당하는 굴절률 변화를 공간적으로 표현한 매우 얇은 두께의 판으로서 전파 거리에 걸쳐 일정 거리마다 배치할 수 있다. 위상판과 위상판 사이는 자유 공간이다. 따라서 식 (1)을 공간 주파수 영역과 공간 영역으로 나누어 계산하면 대기를 전파하는 레이저 빔의 계산이 가능하다.

그림 1은 고출력 레이저 무기의 레이저 빔이 대기 채널을 통과하여 목표물로 조사되는 과정과 이를 시뮬레이션으로 구현하는 방법을 묘사한 그림이다. 고출력 레이저 빔은 대기 난류를 통과하여 목표물에 도달하는데, 이때 레이저 빔은 대기를 전파하면서 대기 난류에 의한 여러 효과(신틸레이션, 빔 원더링 등)를 겪을 뿐 아니라 고출력에 의한 열적 블루밍도 발생시킨다. 이를 시뮬레이션으로 모사하기 위하여 대기 난류에 의한 공간적 굴절률 변화, 그리고 열적 블루밍을 위상판으로 표현한다. Split step 방법은 전파 거리에 이러한 위상판을 주기적으로 위치시켜 레이저 빔의 전파를 계산하는 방법이다. 그림 1에 도시한 바와 같이 구간 ①은 난류 위상판 구간으로 난류로 인한 위상 왜곡이 반영되는 구간이며, 뒤이어 열적 블루밍으로 인한 위상 왜곡이 구간 ②에서 반영된다. 위상판과 위상판 사이의 구간 ③은 자유 공간 전파 구간으로, 빛의 회절 효과와 해당 거리만큼의 대기 흡수 및 산란으로 인한 손실이 반영된다.

Figure 1.Simulation of beam propagation over an atmospheric channel for high-energy laser beam.

2.2. 자유 공간 및 위상판을 통과하는 빔의 전파

위상판과 위상판 사이의 자유 공간에는 굴절률의 변화가 존재하지 않는다. 따라서 식 (1)에서 n(R) = 0이 되므로 아래 식만 남게 된다.

2ikEz+2E=0

식 (2)는 공간 주파수 영역으로, Huygens-Fresnel 적분을 통해 빔의 회절 효과를 계산할 수 있다. Huygens-Fresnel 적분은 식 (3)과 같이 고속 푸리에 변환을 통해 계산할 수 있다.

Ex,y,z+Δz=F1FEx,y,zexpiπλzkx2+ky2+ikΔz

여기에서 ∆z는 자유 공간 전파 구간의 거리, FF−1은 푸리에 변환과 푸리에 역변환, λ는 파장, kxky는 각각 xy방향의 공간 주파수를 의미한다.

자유 공간 전파 이후 빔은 위상판에 도달한다. Split step 방법에서 위상판은 z축 방향의 작은 굴절률 변화를 설명하기 위해 아주 얇게 생성되며 따라서 식 (1)의 두번째 항 ▽2E가 생략되어 아래 식만 남는다.

2ikEz+2k2nRE=0

식 (4)를 통해 난류와 열적 블루밍에 의한 빔의 위상 변화가 계산된다. 자유 공간 전파 및 위상판을 통과하는 과정을 전파 경로 시작 지점부터 마지막 지점까지 순차적으로 계산을 수행하면 전체 전파 경로 동안 레이저 빔이 받는 광학 효과들이 반영된다.

2.3. 흡수 및 산란을 통한 대기 손실

레이저 빔은 대기 중의 고유한 단면적을 가지는 여러 기체 분자 및 에어로졸에 의해 광자 일부가 흡수되거나 산란된다. 반응하는 물질의 입자 종류 및 크기에 따라 레일리 산란(Rayleigh scattering), 미에 산란(Mie scattering), 비선택적 산란(non-selective scattering) 등으로 나뉘며 이 현상으로 인해 목표물로 향하는 레이저 빔의 에너지가 감소한다. 흡수되거나 산란되는 에너지의 양은 일반적으로 파장에 따른 흡수 및 산란 계수로 나타내어지며 이로 인한 레이저 빔의 투과도(transmittance) τ를 식 (5)와 같이 계산할 수 있다.

τ(L)=exp0L α(z)+σ(z)dz

여기에서 ασ는 각각 파장에 따른 흡수 및 산란 계수를, L은 전파 거리를 의미한다. 흡수 및 산란 계수 값은 모두 빛의 파장에 의존한다. 이 계수 값은 이론적으로 도출하는 것이 까다로울 뿐만 아니라, 이론에 포함된 파라미터를 정확하게 측정하는 것이 현실적으로 거의 불가능하므로 통상 가시 거리를 활용하여 해당 파장에 대한 흡수 및 산란 계수를 구해서 사용한다. 그 대표적인 예로 Kim 모델이 있다[12].

2.4. 대기 난류 위상판 생성 방법

대기 난류의 특성을 설명하는 난류의 에너지 캐스케이드 이론(energy cascade theory of turbulence)에 의하면 대기의 풍속은 레이놀즈 수가 임계점을 초과할 때까지 증가한다. 이러한 작용으로 외부 규모(L0) 크기의 불안정한 공기 덩어리인 와류가 형성되며, 형성된 와류는 관성에 의해 더 작은 내부 규모(l0)의 와류들로 분해되면서 받아들인 에너지를 손실 없이 전달한다. 이때 L0의 크기보다 작은 와류는 통계적으로 균질하고 등방성인 것으로 가정한다. 와류의 에너지 전달 과정이 L0 크기의 와류에서 l0 크기의 와류까지 진행되면 받아들인 에너지는 열로 소실되며 최종적으로 생성된 와류는 소멸한다. 와류들은 모두 서로 다른 온도와 압력을 가지는 셀(cell)로 구성되어 있으며, 각 셀의 굴절률 또한 모두 다르다. 이러한 셀로 구성된 와류가 분포된 대기의 굴절률은 시공간적으로 계속 변화하며, 이를 광학 난류라고 한다[11].

시뮬레이션에서 구현하는 광학 난류는 통계적으로 균질하고 등방성을 가진다고 가정되었으며, 생성되는 모든 난류의 형태는 무작위 굴절률 분포를 띄어야 한다. 이는 공간 전력 스펙트럼으로 표현될 수 있으며, 광학 난류를 설명하는 데 대표적으로 사용되는 굴절률 공간 전력 스펙트럼은 난류의 내외부 규모를 고려한 모델인 수정된 본 카르만 모델(modified von Karman spectrum)을 따른다고 알려졌다[11].

Φnκ=0.033Cn2exp(κ2/κm2)κ2+κ0211/6,0κ<;κm=5.93/l0;κ0=2π/L0

식 (6)에서 Cn2는 굴절률 구조상수로서 대기 난류의 강도를 나타내는 값이고, κmκ0은 각각 스칼라 공간 파수로 와류의 내외부 규모를 고려한 상수이다. 또한 κ는 광축에 대하여 수직인 평면의 2차원 공간 파수이다. 식 (6)의 2차원 무작위 위상 분포를 가지는 판을 생성함으로써 대기 채널의 난류를 모사한 난류 위상판을 얻을 수 있다. 난류 위상판 생성 방법으로는 고속 푸리에 변환 방법, 공분산 행렬 방법, Zernike 다항식 방법, 저조파 방법 등이 있으나 본 논문에서는 가장 간단하고 빠르게 계산할 수 있는 고속 푸리에 변환 방법을 사용하였다[13]. 고속 푸리에 변환 방법으로 샘플링하여 구한 위상판의 위치 (x, y)에서의 위상은 아래와 같다.

ΦT(x,y)= κx κys κ x, κ yFϕ κ x, κ yexpi κ xx+ κ yyΔκxΔκy

여기서 s(κx, κy)는 평균 0, 분산 1의 허미시안(Hermitian) 복소 가우시안 과정으로, 이 랜덤 함수를 통해 대기의 무작위성이 부여된다. Fϕ(κx, κy)은 굴절률 변동의 2차원 파워 스펙트럼을 의미한다. 위 과정으로 생성된 난류 위상판은 전체 전파 경로를 동일한 거리 구간으로 나눈 후 각 구간의 중심에 배치되며, 난류 위상판을 통과한 전기장 E'(R)은 다음과 같이 계산된다.

E'(x,y,z)=E(x,y,z)expiΦT(x,y)

2.5. 열적 블루밍 위상판 생성 방법

열적 블루밍은 고출력 레이저(HEL) 빔이 대기로 전파되면서 발생하는 열로 인한 비선형 광학 현상이다. 고출력 레이저 빔이 대기를 통과할 때 대기 중의 특정 분자와 미립자들이 레이저 에너지의 일부를 흡수하는데, 이때 대기에 흡수된 에너지는 주변 대기를 가열하여 팽창시키고, 고출력 레이저 빔 전파 경로상에 분산된 열적 렌즈(thermal lens)를 형성한다. 이렇게 형성된 열적 렌즈를 통과한 고출력 레이저 빔은 너비가 넓어지고, 전파 경로가 휘어진다. 이로 인해 최초 지점에서 완전한 가우시안 분포의 형태였던 레이저 빔은 목표 지점에서 초승달 모양의 형태를 가지게 된다.

열적 블루밍 현상은 제1 열역학 법칙의 에너지 보존 법칙을 바탕으로 내부 에너지만을 고려한 총 에너지 방정식을 굴절률 항으로 나타낸 식 (9)의 에너지 균형 방정식 풀이를 통해 계산할 수 있다[14].

nx,y,z,tt+vznx,y,z,t=1n0αzρCPT0Ix,y,z,t

여기서 t는 시간, v는 대기 이류(advection)의 속도 벡터, ρ는 대기의 밀도, CP는 정압 비열, n0는 공기의 굴절률, T0는 대기의 온도, I(R,t)는 시간 t에 따른 공간 좌표 상의 레이저 빔의 강도(intensity)를 의미한다.

식 (9)의 물리적인 의미는 고출력 레이저 빔이 전파하며 대기에 흡수되는 에너지와 횡방향으로 불어오는 대기 이류로 인해 냉각되는 에너지가 균형을 이루는 것이다. 식 (9)에서 레이저 빔에 영향을 주는 대기 이류 v는 레이저 빔이 전파되는 방향에 수직으로 불어온다고 볼 수 있다. 이때 v(z) = [v(z), 0, 0]으로 나타낼 수 있다. 또한 식 (9)로 설명되는 고출력 레이저 빔 주변 대기의 에너지 변화는 일정 시간 경과 시 레이저 빔의 강도 변화가 거의 없는 정상 상태(steady state)에 도달한다. 따라서 식 (9)의 시간 항을 모두 생략할 수 있으며 이를 굴절률 n(R)에 대한 식으로 정리하면 아래와 같이 거리에 따라 변화하는 굴절률의 양을 이산적으로 계산할 수 있다.

nx,y,z=1n0αzρCPT0vzx I ξ,y,zdξ

이렇게 구해진 열적 블루밍으로 인해 변화하는 굴절률의 양을 거리에 대해서 적분하면, 난류 위상판과 같이 일정 거리에 대하여 발생하는 위상 왜곡 양을 담은 위상판으로 나타낼 수 있으며 이는 식 (11)과 같이 나타낸다.

ΦBx,y,L=k0L 1n 0 αzτ(z) ρCP T 0 vzdzx I ξ,y,0dξ

열적 블루밍의 세기를 나타내는 지표로서 통상 Bradley & Herrmann 왜곡 지수(Nd)를 사용하며 식 (12)와 같이 표현한다[15,16].

Nd(L)=42kP0D00L n 01αzτ(z) ρC PT 0v(z) dz

여기서 P0는 빔의 최초 출력 세기를, D0는 빔의 최초 직경을 의미한다. Bradley & Herrmann 왜곡 지수를 사용하면 식 (11)은 아래와 같이 간단히 표현된다.

ΦBx,y,L=Nd(L)D042kP0x I ξ,y,0dξ

난류 위상판과 마찬가지로 위 과정으로 생성된 열적 블루밍 위상판은 동일한 거리로 나누어진 각 구간의 중심에 배치되며, 열적 블루밍 위상판을 통과하는 전기장 E'(R)을 계산하면 식 (14)와 같다.

E'(x,y,z)=E(x,y,z)expiΦB(x,y)

그림 2는 2.4절과 2.5절에서 설명한 위상판 생성 방법으로 생성한 위상판의 예시로, 그림 2(a)는 난류 위상판, 그림 2(b)는 열적 블루밍 위상판이다. 열적 블루밍 위상판 생성 간 대기 이류는 2 m/s로 설정하였다. 난류 위상판은 전체 대기에 서로 다른 크기의 와류가 통계적으로 균질하게 분포되어 있으며, 난류 위상판을 통해 난류로 인한 위상 왜곡량이 레이저 빔에 반영된다. 마찬가지로 열적 블루밍 위상판을 통해 대기 이류 방향으로 축적된 위상 왜곡량 또한 전파되는 레이저 빔에 반영된다. 이때 주목할 점은 식 (7)의 s(κx, κy) 함수로 인하여 난류 위상판에 무작위성이 부여되는 반면, 열적 블루밍 위상판은 레이저 전력에 따라 변화할 뿐 무작위성은 없다는 것이다.

Figure 2.Example of phase screens for (a) atmospheric turbulence and (b) thermal blooming.

Split step 방법을 활용한 레이저 빔 전파 시뮬레이션에서 레이저 빔은 2차원 공간 내 특정 크기의 판 안에 표현되며, 위상판 또한 이와 같은 크기로 생성된다. 이러한 특성으로 인해 시뮬레이션이 타당성을 가지기 위해서는 시뮬레이션에 사용된 주요 파라미터들이 특정 조건을 만족해야 한다. 또한 사용하는 각 효과의 위상판 개수에 따라 빔의 산란 및 왜곡되는 정도가 달라지므로, 적절한 개수의 위상판을 판단하는 것이 중요하다.

3.1. Split step 방법을 활용한 시뮬레이션의 조건

Split step 방법을 활용한 레이저 빔 전파 시뮬레이션은 널리 이용되는 방법이기 때문에 주요 파라미터 설정 기준과 조건은 잘 알려져 있다[7,8]. 우선 생성되는 위상판은 와류의 외부 규모와 내부 규모의 영향을 충분히 반영할 수 있어야 한다. 따라서 위상판 크기 A가 외부 규모보다 충분히 커야 하며, 이를 위해 A > L0를 만족해야 한다. 또한 위상판 격자의 크기 ∆x∆y는 내부 규모보다 충분히 작아야 하므로, 각각 ∆x < l0/2와 ∆y < l0/2를 만족해야 한다. 두번째로, 나이퀴스트 샘플링 정리(Nyquist’s sampling theorem)를 만족하기 위해 위상판 내의 이웃한 격자 간의 위상 차이는 π보다 작아야 한다. 세번째로, 빔의 전파 계산에 사용되는 푸리에 변환에서 신호가 왜곡되는 에일리어싱(aliasing) 현상을 방지하기 위해 위상판 간 거리는 기본적으로 ∆z < kA∆xπ를 만족해야 한다. 네번째로, split step 방법에서는 격자(grid)의 불연속적 특성으로 인하여 가장자리에서 빔이 산란되어 판을 벗어나는 경우 옆부분에 해당 빔의 에너지가 표현되는 가장자리 효과(edge effect)가 발생하는데, 이때 가장자리 효과를 최소화하기 위하여 A > 2θmax ∆z를 만족해야 한다. 여기서 θmax = max[√2/k × ∂f/∂x]로서 빔 전파간 빔의 산란 각도의 최대값이다. f는 난류 위상판의 특정 격자의 위상값을 의미한다.

3.2. 대기 난류를 모사하기 위하여 필요한 위상판 개수

레이저 빔이 대기 난류를 통과하면서 변동하는 강도는 전통적으로 약한 변동 이론과 강한 변동 이론으로 분류된다. 여기서 Rytov 분산 값은 변동의 강약을 구분하는 값이며 식 (15)와 같이 쓸 수 있다[11].

σR2=1.23Cn2k7/6L11/6

빔 전파 시작 지점부터 σR2 < 1을 만족하는 거리까지를 빔 변동 강도가 약한 구간(weak fluctuation)이라고 한다. 대기 난류를 통과하는 레이저 빔의 변동 강도는 신틸레이션 지수로 나타낼 수 있으며, σR2 ≥ 1일 때 신틸레이션 지수는 포화 상태에 이른다. 집속된 가우시안 빔에 대한 신틸레이션 지수의 이론식은 잘 알려져 있으며, 포화 상태에 이르기 전의 값은 식 (16)으로 구할 수 있다[11,17].

σIT20,0,L=4.42σR2Γ5/6 σ pe DB 2+3.86σR2Rei5/6 2F156 ,116 ;176 ;1Θ+iΓ1116Γ5/6,σR2<1

여기서 Г = Г0/(1 + Γ02), Г0 = 8L/kD02이며 가우시안 빔 매개변수를 의미한다. 2F1(a, b; c; z)은 가우스 초기하 함수(Gaussian hypergeometric function)이며 Θ = 1/(1 + Γ02)은 빔 곡률 매개변수를 의미한다. DB는 아무런 왜곡 없이 목표 거리에 도달한 이상적인 가우시안 빔의 회절 한계를 고려한 빔의 직경이다. σpe2는 지터(jitter)로 인한 포인팅 오류 분산값으로 식 (17)과 같다.

σpe2~ λLD02 D0r05/3,D0r0<<1 λLD02 r0D01/3,D0r0<<1r0= 0.423Cn 2 k 2L 3/5

여기서 r0는 프라이드 파라미터(Fried’s parameter)로 대기 난류의 영향 정도를 나타내는 척도이다. 또한 신틸레이션 지수는 아래와 같이 시뮬레이션을 통해 생성된 거리별 전파 경로상 중심에서의 빛의 강도 표본값을 바탕으로 계산할 수 있다.

σI2x,y,L=I2x,y,LIx,y,L21

여기서 < >는 앙상블 평균값을 의미한다. σR2 < 1일 때, 해당 구간에서 신틸레이션 지수의 변화 정도를 제한함으로써 레이저 빔의 대기 난류 전파를 정확하게 모사할 수 있다. 이때 기준은 다음과 같다. 우선 난류 위상판 1개를 통과할 때 증가하는 신틸레이션 지수는 각 위상판이 충분히 작은 위상 변동만 일으켜야 한다는 조건을 만족시키기 위하여 0.1을 초과할 수 없다[8]. 이 기준은 대기 전파 시뮬레이션 연구에서 종합적으로 경험에 의해 얻어진 수치로서 널리 사용되고 있으며 이를 만족하는 위상판 개수(NT1)는 아래와 같이 나타낼 수 있다.

NT1=10LσIT20,0,1.23Cn2k7/66/111.23Cn2k7/66/11

또한 전체 거리를 전파하며 변화하는 신틸레이션 지수가 작더라도, 충분한 개수의 위상판을 사용하기 위한 조건을 만족시키기 위하여 난류 위상판 1개를 통과할 때 변화하는 신틸레이션 지수의 변화량은 전체 거리를 전파하는 동안 변화하는 신틸레이션 지수의 총 변화량의 10% 이하여야만 한다[8,18]. 위상판과 위상판의 거리는 모두 동일하므로, 각 위상판이 변화시키는 위상 변화의 총량 또한 동일하다. 따라서 위 조건을 만족하는 위상판 개수(NT2)는 아래와 같다.

NT2=10

위의 두 조건을 모두 만족하는 경우 전체 경로에 걸쳐 대기 난류를 모사하는 데 필요한 위상판의 최소 개수(NT)를 구할 수 있으며, 아래와 같이 표현할 수 있다.

NT=maxNT1,10

3.3. 열적 블루밍을 모사하기 위하여 필요한 위상판 개수

열적 블루밍의 경우 전파 거리 L에 대하여 고출력 레이저 빔의 중심에 주는 왜곡량의 최댓값은 아래와 같이 표현된다[14].

ΦB_max=12πNdL

또한 총 M개의 열적 블루밍 위상판을 사용하면 m번째 열적 블루밍 위상판이 가지는 왜곡의 양은 식 (23)과 같다.

ΦBx,y,mLM=α(z)k1n0eα(z)+σ(z)CPT0ρv(z)α(z)+σ(z)emL/Mem1L/Mx I ξ,y,0dξ

본 논문에서는 열적 블루밍에 필요한 위상판 개수를 찾기 위해 고출력 레이저 빔의 성능을 측정할 때 주로 사용되는 power in the bucket (PIB) 지표를 활용하였다. PIB는 빔 전파 최초 지점에서 지름이 D0인 원형 버킷 안에 있는 에너지의 총량과, 전파된 후 목표 거리의 중앙에서의 회절 한계를 고려한 지름 DB의 원형 버킷 안에 들어간 에너지의 총량을 측정하여 식 (24)와 같이 계산한다.

PIB= DB/2 DB/2Ix,y,Ldxdy D0/2 D0/2Ix,y,0dxdy

그림 3은 열적 블루밍 세기와 위상판 개수를 변화시켜가며 측정한 PIB 값을 표준화한 결과이다. 그래프에 표시된 흰 점선은 모든 열적 블루밍 세기 범위에서 표준화한 PIB 값이 1로 수렴하는 데 필요한 최소 개수의 위상판 개수를 나타낸 것이다. 그림 3에 따르면 위상판의 개수에 따라 PIB 값의 차이가 큰 것을 볼 수 있는데, 이는 위상판 간 거리가 멀어서 자유 공간 전파 중에 빔이 많이 분산되기 때문이다. 반면 충분한 개수의 위상판이 사용되면 표준화한 PIB 값이 1로 수렴하는데, 이 지점에서 위상판이 가지는 위상 왜곡량의 비율을 식 (22)와 식 (23)을 통해 측정함으로써 위상판 한 개가 가져야 할 위상 왜곡량이 총 왜곡량의 5% 미만임을 확인하였다. 이 조건을 만족하는 위상판 개수(NB1)는 식 (25)와 같이 나타낼 수 있다.

Figure 3.Normalized power in the bucket value as functions of distortion number, Nd, and number of thermal blooming screens NB.

NB1=52πD0P0eαz+σ(z)LΔzD0 I ξ,0,0dξ

또한 나이퀴스트 샘플링 정리(Nyquist sampling theorem)를 만족하기 위해 열적 블루밍 위상판의 위상차는 π를 넘을 수 없다[7]. 이를 만족하는 위상판 개수(NB2)는 아래와 같이 나타낼 수 있다.

NB2=kn01α(z)LCPT0ρv(z)πeα(z)+σ(z)ΔzD0 I ξ,0,0dξ

앞선 두 조건을 모두 만족시키면 전체 경로에 걸쳐 열적 블루밍을 모사하는 데 필요한 위상판의 최소 개수(NB)를 구할 수 있으며, 이는 식 (27)과 같이 표현할 수 있다.

NB=maxNB1,NB2

3.4. 대기 난류와 열적 블루밍이 함께 고려된 시뮬레이션에 필요한 위상판 개수

지금까지 대기 난류와 열적 블루밍이 개별적으로 고려된 시뮬레이션에 필요한 위상판 개수를 구하였다. 그러나 실제 대기 공간에서는 고출력 레이저 빔 전파 시 난류와 열적 블루밍이 동시에 발생한다. 본 논문에서는 계산의 용이성과 시뮬레이션의 단순화를 위하여 그림 1과 같이 난류 위상판과 열적 블루밍 위상판을 전파 경로상에 함께 배치하였다. 또한 이때 사용되는 각 효과의 위상판 개수는 시뮬레이션 조건(난류 세기, 열적 블루밍 세기 등)에 따라 달라지는 위상판 개수 중 가장 큰 값(Nps)을 따르는 것으로 설정하였다. 이를 수식으로 표현하면 식 (28)과 같다.

Nps=maxNT,NB=maxNT1,10,NB1,NB2

또한 난류 위상판과 열적 블루밍 위상판을 통과하는 전기장 E'(R)의 계산 결과는 아래와 같다.

E'(x,y,z)=E(x,y,z)expiΦT(x,y)expiΦB(x,y)

이번 장에서 우리는 150개의 위상판을 사용하는 기존 방법과 본 논문에서 제안하는 방법의 시뮬레이션 정확도 및 소요 시간을 비교하였다. 시뮬레이션에 사용된 주요 매개변수의 값은 표 1에 요약되어 있다. 먼저, 시뮬레이션은 3장에서 제시한 기준을 모두 만족하는 조건에서 수행하였다. 고출력 레이저 빔은 집속된 가우시안 빔을 사용한다고 가정하며 전송 거리는 5 km이다. 난류 위상판은 식 (5)의 수정된 본 카르만 모델을 사용하여 생성하고, 와류의 내외부 규모는 각각 4 mm, 1.7 m로 설정하였다. 위상판의 크기는 4 × 4 m2이고, 이를 4,096 × 4,096개의 격자 포인트 해상도로 구현하였다. 또한 대기 이류의 속도는 한국 45개 지역에서 측정된 연평균 풍속인 2 m/s로 설정하였으며, 난류 세기(굴절률 구조상수) 또한 한국 도시별 굴절률 구조상수의 값을 고려하여 10−16–10−13 m−2/3으로 설정하였다[19].

Table 1 Simulation parameters

ParameterValue
Simulation Grid4,096 × 4,096
Side Length of Phase Screen (m)4
Wavelength (nm)1,070
Initial Beam Diameter (m)0.5
Propagation Distance (km)5
Transmitter Beam Power (kW)1–1,000
Absorption Coefficient (m−1)5 × 10−6
Scattering Coefficient (m−1)5 × 10−5
Refractive Index of Air1.0003
Specific Heat at Constant Pressure (J/kg/K)1,004
Density of Air at Constant Pressure (kg/m3)1.293
Initial Air Temperature (K)300
Refractive Index Structure Constant (m−2/3)10−13–10−16
Advection Wind Speed (Direction) (m/s)2 (Left to Right)
Beam Propagation DirectionHorizontal Direction


두 시뮬레이션의 정확성을 비교하기 위하여 난류와 열적 블루밍 세기별로 총 10,000회의 시뮬레이션을 시행하였다. 이는 무작위 특성의 난류와 열적 블루밍이 빔에 끼치는 영향을 통계학적으로 분석하기 위함이다. 이렇게 획득한 10,000개의 표본에 대해 PIB를 식 (24)를 사용하여 계산하였다.

다음으로 PIB 값 분포를 통해 구해진 PMF의 유사도를 결정계수(R-squared)를 측정하여 비교하였다.

R2=1 i=1Mgihi2 i=1Mgig¯2

여기서 gi는 기존 방법으로 구한 PIB 분포 그래프에서 각 bin의 발생 확률이며 gi의 평균값을, hi는 제안된 방법으로 구한 PIB 분포 그래프의 각 bin의 발생 확률을 의미한다. 결정계수 값이 0.9 이상이면 통상 통계적 유사도가 높은 것으로 볼 수 있다[20,21]. 반면 결정계수 값이 낮으면 오차가 많으므로 두 분포의 유사도가 낮다고 본다.

또한 전파되는 동안 거리별 전파 경로 축상 신틸레이션 지수(on-axis σI2)를 식 (18)을 사용하여 계산하고, 두 시뮬레이션의 on-axis σI2의 평균 오차(η)를 다음과 같이 구하였다.

η= i=1MσIc2(0,0,zi)σIp2(0,0,zi)M

여기서 σIc2σIp2는 각각 기존 방법과 제안된 방법으로 계산된 신틸레이션 지수이다.

그림 4는 모두 열적 블루밍의 강도(Nd)가 50 rad인 조건에서 난류 세기를 달리했을 때 목표 거리에서의 고출력 레이저 빔 형태의 예시이다. 화면 중앙의 흰색 원은 PIB 값을 측정하기 위한 3 cm 지름의 원형 버킷이다. 그림 4(a)는 난류 세기가 약한 대기 채널(Cn2 = 10−16 m−2/3)을 통과했을 때의 모습이다. 난류에 의한 신틸레이션과 빔 원더링 현상이 또렷이 확인되며, 열적 블루밍에 의해 빔의 형태가 대기 이류의 역방향으로 초승달처럼 휜 것을 확인할 수 있다. 그림 4(b)는 난류 세기가 강할 때(Cn2 = 10−13 m−2/3) 대기를 통과한 빔의 형태이다. 그림 4(a)보다 신텔레이션과 빔 원더링 현상이 강하게 발현되었으나 열적 블루밍은 상대적으로 미미한 것을 알 수 있다. 이는 강한 난류로 인해 빔이 많이 퍼져 있기 때문이다.

Figure 4.Exemplary beam intensity profiles after 5-km propagation over an atmospheric channel when the distortion number is 50 rad: (a) weak turbulence (Cn2 = 10−16 m−2/3) and (b) strong turbulence (Cn2 = 10−13 m−2/3). The white circle represents a bucket having a 3-cm diameter.

그림 5는 기존 방법 대비 제안된 방법의 PIB 분포별 PMF 결과로, 많은 횟수의 시뮬레이션을 시행했을 때 두 방법 간의 확률적인 차이를 비교할 수 있다. 이때 PMF의 bin 개수는 결정계수 값 비교를 위하여 난류 및 열적 블루밍 세기와 상관없이 50개로 설정하였다. 그림 5(a)는 난류 세기가 작고(Cn2 = 10−15 m−2/3), 열적 블루밍의 세기는 126 rad으로, 이때 결정계수는 0.9749이다. 그림 5(b)는 난류 세기가 강하고(Cn2 = 10−13 m−2/3), 열적 블루밍의 세기가 84 rad인 경우이며, 이때의 결정계수는 0.9906이다. 따라서 두 경우 모두 기존의 방법과 제안된 방법의 PIB 분포별 PMF 히스토그램의 결과가 잘 일치함을 알 수 있다.

Figure 5.Probability mass functions of power in the bucket values under (a) weak turbulence (Cn2 = 10−15 m−2/3) when the distortion number is 126 rad and (b) strong turbulence (Cn2 = 10−13 m−2/3) when the distortion number is 84 rad. The total number of simulations is 10,000.

표 2는 조건을 다르게 하여 측정한 결정계수와 on-axis σI2 의 평균 오차값의 결과이다. 표 2에서 레이저 출력 밀도는 레이저 출력과 빔 단면적으로 계산된 평균값을 사용하였다. 그림 6(a)표 2의 결정계수 결과를, 그림 6(b)표 2의 on-axis σI2 의 평균 오차값의 결과를 레이저 출력 밀도에 대한 그래프로 나타낸 것이다. 결과를 보면 열적 블루밍 세기(레이저 출력)가 커질수록, 난류 세기는 작아질수록, 결정계수 값이 점차 감소하는 것을 볼 수 있다. 이는 난류 세기가 약할수록 열적 블루밍이 상대적으로 빔의 요동 및 형태 변화에 더 큰 영향을 미치기 때문이다. 열적 블루밍의 강도가 일정량 이상이 되면, 난류와 열적 블루밍의 상호작용에 의한 섭동(thermal Rayleigh scattering)이 추가되어 본 논문에서 제시한 기준보다 더 많은 개수의 위상판이 필요하다[22]. 반면 난류 세기가 강하면 난류가 열적 블루밍보다 상대적으로 빔의 변화에 더 큰 영향을 미치며, 빔의 섭동이 열적 블루밍에 의해 크게 영향을 받지 않는다. on-axis σI2 의 평균 오차값은 열적 블루밍 세기가 210 rad (레이저 출력 밀도 2.5 × 106 W/m2)를 초과할 경우 값이 현저하게 증가함을 알 수 있다. 따라서 기존 방법의 시뮬레이션을 기준으로 할 때, 우리가 제시한 위상판 개수는 모든 난류 세기 강도(Cn2 = 10−16 – 10−13 m−2/3)에서 레이저 출력 밀도가 2.5 × 106 W/m2 이하일 때 사용하기 적합하다.

Figure 6.(a) R-squared and (b) average error of on-axis scintillation indexes for various atmospheric conditions and laser power densities.

Table 2 R-squared and average error of on-axis scintillation indexes (η) for various atmospheric conditions and laser output powers

P0(Nd)1 kW
(0.4 rad)
10 kW
(4 rad)
50 kW
(21 rad)
100 kW
(42 rad)
200 kW
(42 rad)
300 kW
(126 rad)
500 kW
(210 rad)
600 kW
(252 rad)
700 kW
(295 rad)
800 kW
(337 rad)
900 kW
(379 rad)
1 MW
(420 rad)
Cn210−16
(m−2/3)
R20.98590.97750.99980.95550.9428-0.90720.9001---0.8064
η0.00210.00200.00590.00820.0084-0.06500.0830---0.1983
10−15
(m−2/3)
R20.99000.98370.9723-0.97580.97490.95370.94160.9325---
η0.01240.00970.0090-0.00850.00800.09520.14250.1948---
10−14
(m−2/3)
R20.98360.99440.98530.9927--0.96370.9455-0.8923--
η0.05000.07260.08250.0708--0.08820.1130-0.1783--
10−13
(m−2/3)
R20.99250.9962--0.9906-0.98050.9638--0.9617-
η0.06500.0860--0.0689-0.08250.0868--0.0990-
Laser Power Density (kW/m2)5.15.1 × 102.5 × 1025.1 × 1021.0 × 1031.5 × 1032.5 × 1033.1 × 1033.6 × 1034.1 × 1034.6 × 1035.1 × 103


대기 난류와 열적 블루밍이 함께 고려된 레이저 빔 전파 시뮬레이션의 과정은 크게 난류 위상판 및 열적 블루밍 위상판의 생성, 자유 공간 전파, 위상판 통과 및 기타 부수 과정(PIB 값 계산, 난류 세기 및 열적 블루밍 세기 연산 과정 등)으로 구분할 수 있다. 따라서 시뮬레이션 총 진행 횟수(N)에 따른 총 소요 시간(Ttotal)은 아래와 같이 계산할 수 있다.

Ttotal=Nps×Ttb+Ttul+Tprop+Tpass+Λ×N

여기서 TtbTtul은 각각 열적 블루밍 및 난류 위상판을 생성하는 시간, Tprop은 자유 공간 전파 구현시간, Tpass는 각 위상판 통과 연산 시간, Λ는 기타 부수 과정에 소요되는 시간을 의미한다. 대부분의 시간은 TtbTtul에서 소요된다.

식 (32)에서 알 수 있듯이 위상판의 개수가 증가할수록 시뮬레이션 내에서 자유공간 전파 시간과 각 효과의 위상판 생성 횟수도 그와 비례하여 증가하며, 전체 시뮬레이션에 소요되는 시간 또한 선형적으로 증가하게 된다. 이를 바탕으로 기존의 방법과 본 논문에서 제안한 방법을 사용한 두 시뮬레이션의 소요 시간을 비교하였으며 이를 그림 7에 나타내었다. 두 방법 모두 클락 속도 4.7 GHz의 Intel® CoreTM i7-8,700K 프로세서를 갖춘 개인용 컴퓨터에서 실행되었다.

Figure 7.(a) Total simulation time measured by using a personal computer versus the total number of simulations, N, under strong turbulence (Cn2 = 10−13 m−2/3) and weak turbulence (Cn2 = 10−16 m−2/3). (b) Total simulation time versus the refractive index structure constant, Cn2. The total number of simulations is 10,000.

그림 7(a)는 시뮬레이션 총 진행 횟수(N)에 따른 소요 시간을 나타낸다. 기존의 방법은 모든 난류 세기에 대해서 동일한 시간이 소요되었으나, 제안된 방법의 경우 난류 세기에 따라 소요 시간이 다르므로 약한 난류 세기(Cn2 = 10−16 m−2/3)일 때와 강한 난류 세기(Cn2 = 10−13 m−2/3)일 때의 소요 시간을 따로 나타내었다. 시뮬레이션 결과, 제안된 방법을 사용할 시 기존의 방법에 비하여 총 소요 시간이 최대 6배 감소하였다. 예를 들면 시뮬레이션 총 실행 횟수가 100,000회일 때, 제안된 방법은 난류 세기가 강할 경우 약 200일이, 난류 세기가 약할 경우에는 약 54일이 걸리지만, 기존의 방법은 모든 난류 조건에서 330일이 걸린다. 이는 제안된 방법에서 필요 위상판 개수에 따라 위상판 생성 및 자유공간 전파 횟수가 현저하게 줄어들기 때문이다.

그림 7(b)는 난류의 세기인 굴절률 구조상수 값을 달리해 가며 시뮬레이션 10,000회를 시행하는 시간을 나타낸다. 제안된 방법은 난류의 세기가 약할수록 필요한 위상판 개수가 적어지므로 더 적은 시간이 소요되며, 모든 난류 세기를 기준으로 보았을 때에조차 기존의 방법보다 소요 시간이 적은 것을 알 수 있다. 결과적으로 제안된 방법은 우리나라의 도시별 평균 최대 난류 세기(Cn2 = 10−13 m−2/3) 이하의 범위에서 실행 횟수가 많을수록, 그리고 난류 세기가 약할수록 기존 방법보다 시간 면에서 효율적이다.

본 논문에서는 대기 난류와 열적 블루밍을 겪는 고출력 레이저 빔의 대기 전파를 빠르고 정확하게 시뮬레이션 수행하기 위하여 요구되는 위상판의 개수를 제시하였다. Split step 방법을 사용한 빔 전파 시뮬레이션의 시간 대부분은 자유공간 빔 전파와 위상판 생성에 소요되기 때문에, 시뮬레이션에 사용되는 위상판의 총 개수를 줄임으로써 시뮬레이션의 시간을 크게 단축할 수 있다. 먼저 전파 거리, 대기 난류의 세기, 레이저 빔의 출력 전력 등에 따라서 각 위상판이 레이저 신호에 미치는 적절한 위상의 양을 결정함으로써 정확한 시뮬레이션 구현에 요구되는 최소 위상판의 개수를 도출하였다. 본 논문에서는 레이저 출력 밀도 2.5 × 106 W/m2 (빔 직경 50 cm 레이저의 경우 출력 전력 500 kW) 미만의 빔이 대기를 전파할 때 이를 모사하기 위하여 필요한 최소 위상판의 개수를 제시하였으므로, 이 방법은 고출력 레이저 빔의 대기 전파를 통계적으로 분석할 때 유용하게 사용될 것으로 기대된다.

본 연구의 결과 분석 및 생성된 데이터는 모두 본 논문 내 명시되어 있으며 공공의 이용이 가능하다. 데이터에 접근하거나 사용하고자 하는 이는 저자에게 타당한 이유를 밝히고 허가를 득해 사용 가능하다.

  1. P. E. Ross, "Economics drives a ray-gun resurgence: Lasers, cheaper by the shot, should work well against drones and cruise missiles," IEEE Spectrum 60, 40-41 (2023).
    CrossRef
  2. B. Zohuri, Directed Energy Weapons (Springer Charm, USA, 2016), 1st ed..
    CrossRef
  3. P. E. Nielsen, Effects of Directed Energy Weapons (CreateSpace Independent Publishing, USA, 2012), 1st ed..
  4. W. Moon and H. Kim, "Standard deviation of fiber-coupling efficiency for free-space optical communication through atmospheric turbulence," IEEE Photonics J. 15, 7302507 (2023).
    CrossRef
  5. R. Frehlich, "Simulation of laser propagation in a turbulent atmosphere," Appl. Opt. 39, 393-397 (2020).
    Pubmed CrossRef
  6. Z. Chen, D. Zhang, C. Xiao, and M. Qin, "Precision analysis of turbulence phase screens and their influence on the simulation of Gaussian beam propagation in turbulent atmosphere," Appl. Opt. 59, 3726-3735 (2020).
    Pubmed CrossRef
  7. D. L. Knepp, "Multiple phase-screen calculation of the temporal behavior of stochastic waves," Proc. IEEE 71, 722-737 (1983).
    CrossRef
  8. J. M. Martin and S. M. Flatté, "Intensity images and statistics from numerical simulation of wave propagation in 3-D random media," Appl. Opt. 27, 2111-2126 (1988).
    Pubmed CrossRef
  9. M. F. Spencer, "Wave-optics investigation of turbulence thermal blooming interaction: II. Using time-dependent simulations," Opt. Eng. 59, 081805 (2020).
    CrossRef
  10. C. E. Murphy and M. F. Spencer, "Investigation of turbulence thermal blooming interaction using the split-step beam propagation method," Proc. SPIE 10772, 1077208 (2018).
  11. L. C. Andrews and R. L. Phillips, Laser Beam Propagation Through Random Media, 2nd ed. (SPIE Press, USA, 2005).
    CrossRef
  12. I. I. Kim, B. McArthur, and E. J. Korevaar, "Comparison of laser beam propagation at 785 nm and 1550 nm in fog and haze for optical wireless communications," Proc. SPIE 4214, 26-37 (2001).
  13. D. T. Ha, V. V. Mai, and H. Kim, "Comparison of phase-screen-generation methods for simulating the effects of atmospheric turbulence," Korean J. Opt. Photonics 30, 87-93 (2019).
  14. F. G. Gebhardt, "Twenty-five years of thermal blooming: an overview," Proc. SPIE 1221, 2-25 (1990).
    CrossRef
  15. L. C. Bradley and J. Herrmann, "Phase compensation for thermal blooming," Appl. Opt. 13, 331-334 (1974).
    Pubmed CrossRef
  16. P. A. Konyaev and V. P. Lukin, "Thermal distortions of focused laser beams in the atmosphere," Appl. Opt. 24, 415-421 (1985).
    Pubmed CrossRef
  17. L. C. Andrews, R. L. Philips, R. J. Sasiela, and R. Parenti, "Beam wander effects on the scintillation index of a focused beam," Proc. SPIE 5793, 28-37 (2005).
    CrossRef
  18. J. M. Martin and S. M. Flatté, "Simulation of point-source scintillation through three-dimensional random media," J. Opt. Soc. Am. A 7, 838-847 (1990).
    CrossRef
  19. V. V. Mai, D. T. Ha, and H. Kim, "Link availability of terrestrial free-space optical communication systems in Korea," Korean J. Opt. Photonics 29, 77-84 (2018).
  20. E. Ostertagová, "Modelling using polynomial regression," Procedia Engineering 48, 500-506 (2012).
    CrossRef
  21. H. Kotake, Y. Abe, D. R. Kolev, Y. Saito, Y. Takahashi, T. Fuse, Y. Satoh, T. Itahashi, S. Yamakawa, H. Tsuji, and M. Toyoshima, "Experimental analysis of atmospheric channel model with misalignment fading for GEO satellite-to-ground optical link using 'LUCAS' onboard optical data relay satellite," Opt. Express 31, 21351-21366 (2023).
    Pubmed CrossRef
  22. T. J. Karr, J. R. Morris, D. H. Chambers, J. A. Viecelli, and P. G. Cramer, "Perturbation growth by thermal blooming in turbulence," J. Opt. Soc. Am. B 7, 1103-1124 (1990).
    CrossRef

Article

연구논문(Research Paper)

2024; 35(2): 49-60

Published online April 25, 2024 https://doi.org/10.3807/KJOP.2024.35.2.049

Copyright © Optical Society of Korea.

Number of Phase Screens Required for Simulation of a High-energy Laser Beam’s Propagation Experiencing Atmospheric Turbulence and Thermal Blooming

Seokyoung Yoon, Woohyeon Moon, Hoon Kim

School of Electrical Engineering, Korea Advanced Institute of Science and Technology, Daejeon 34141, Korea

Correspondence to:hoonkim@kaist.ac.kr, ORCID: 0000-0001-7395-3695
Current affiliation: Republic of Korea Army, Gyeryong 32800, Korea

Received: December 26, 2023; Revised: January 23, 2024; Accepted: January 31, 2024

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

We analyze the number of phase screens required for the simulation of a high-energy laser beam’s propagation over an atmospheric channel. For high-energy lasers exceeding tens of kilowatts (kW) in power, the laser beam is mainly affected by atmospheric turbulence and thermal blooming. When using the split-step method to implement losses due to atmospheric absorption and scattering and distortion of the beam due to turbulence and thermal blooming, the number of phase screens is a critical factor in determining the accuracy and time required for the simulation. By comparing simulation results obtained using a large number of phase screens (e.g., 150 screens) under a wide range of atmospheric turbulence conditions, we provide new guidelines for the number of phase screens required for simulating the beam propagation of a high-power laser below 2.5 × 106 W/m2 (e.g., a 500-kW laser beam having a 50-cm diameter).

Keywords: High-energy laser, Laser weapon, Thermal blooming, Turbulence

I. 서 론

고출력 레이저(high-energy laser, HEL) 무기는 빛의 속도로 전파되며, 원하는 목표물에 정확히 조사되고, 미사일이나 포탄에 비교할 때 발사 비용이 저렴한 동시에 탄약의 수에 구애받지 않는다는 다양한 장점이 있다. 이러한 장점으로 인해 고출력 레이저 무기는 적의 전 영역에 걸친 군사적 위협에 대응할 수 있는 미래 핵심 전력으로서 전 세계가 주목하고 있다[1]. 그러나 레이저 빔은 조사 후 대기를 전파하는 중에 대기 매질의 굴절률 변화 및 대기와의 상호작용으로 인하여 빔의 세기와 품질 등이 변화하므로 무기의 성능이 감소할 수 있다. 따라서 레이저 무기의 성능을 극대화하기 위해서는 고출력 레이저 빔의 대기 전파 특성을 정확하게 이해하는 것이 필수이다.

레이저 빔은 대기 전파 중 흡수와 산란을 통해 단위 면적당 강도가 감소한다. 또한 대기에 존재하는 불균일한 공기 입자와 시시각각 변화하는 난류로 인하여 신틸레이션(scintillation), 빔 원더링(beam wandering)과 같은 선형적 효과가 발생하고, 그로 인해 목표물에 도달하기 전 에너지가 분산되거나 목표물 조준 오차가 발생하여 타격의 효율이 낮아진다[2]. 또한 레이저 출력이 킬로와트(kW) 급 이상으로 커지면, 강한 빛으로 인해 대기 물질이 열을 흡수하며 대기의 굴절률을 변화시키는 열적 블루밍(thermal blooming) 현상이 발생한다. 열적 블루밍 현상으로 인해 레이저 빔은 공간적으로 퍼지고, 전파 방향이 왜곡되어 빔의 조사가 부정확해진다[3].

레이저 빔이 대기를 전파하며 겪는 광학 현상을 이해하고, 이러한 현상이 레이저 빔 조사에 미치는 영향을 분석하기 위하여 통상적으로 컴퓨터 시뮬레이션이 활용된다[4,5]. 특히 시뮬레이션을 활용하면 실제 실험에서 조정하기 어려운 빔의 출력 세기, 난류의 세기, 전파 거리 등의 여러 조건을 손쉽게 변화시키면서 무기의 성능을 예측할 수 있다. 또한 시뮬레이션을 반복 실행함으로써 난류와 열적 블루밍의 무작위 특성이 빔에 끼치는 영향을 통계학적으로 분석할 수 있다[6]. 예를 들어 고출력 레이저 빔의 목표물에서의 형태, 강도, 조준점은 같은 난류 세기 조건에서도 난류의 분포 형태, 그리고 시간에 따라서 지속적으로 변한다. 따라서 정확한 성능 분석을 위해서 목표 거리에서 조준점에 모인 에너지의 양(power in the bucket, PIB)의 확률 질량 함수(probability mass function, PMF) 등을 활용하여 성능을 분석해야 한다. 하지만 통계학적 특성을 파악하기 위하여 시뮬레이션을 대량으로 반복 실행할 경우, 전체 시뮬레이션 시간이 반복 횟수에 비례하여 크게 증가한다. 따라서 시뮬레이션의 정확성을 희생하지 않으면서도 시뮬레이션 시간을 줄이는 것이 중요하다.

레이저 빔의 대기 전파는 확률적 파동 방정식(stochastic wave equation)으로 기술되며 이 방정식은 통상 split step 방법을 활용하여 계산된다[4]. Split step 방법은 빔의 전파 경로상에 특정 전파 거리에 해당하는 위상 왜곡 양을 가지는 위상판(phase screen)을 일정한 간격으로 배치하고, 위상판과 위상판 사이는 자유 공간으로 구성한다. 위상판은 대기의 난류나 열적 블루밍 특성에 맞추어 공간적으로 다른 위상을 가진 매우 얇은 판으로, 특정 거리의 대기를 전파할 때 발생하는 위상 변화를 모사하였다. 대기 난류를 전파하는 레이저 빔을 이와 같은 방법으로 계산할 때 고려해야 하는 사항은 이미 여러 문헌에서 보고되었는데, 예를 들면 난류 위상판 생성 간 고려되어야 할 와류(eddy)의 내외부 규모(inner and outer scale)에 따른 위상판의 크기, 에일리어싱(aliasing) 방지를 위한 위상판 사이의 거리 및 한 개의 위상판이 가져야 할 난류의 세기 정도에 대한 기준이 알려져 있다[7,8]. 열적 블루밍의 모사 또한 위상판의 형태로 전파 경로에 난류 위상판과 나란히 배치하는 방법이 연구된 바 있다[9]. 하지만 대기 난류 현상과 열적 블루밍을 함께 고려한 기존의 시뮬레이션 연구들은 두 효과를 정확히 모사하는 것에만 집중하고 있다. 그 결과 난류와 열적 블루밍을 모사하는 위상판의 개수가 많을수록 좋다는 것 이외에 시뮬레이션 시간을 줄이면서 두 현상을 효과적으로 모사하는 데 필요한 위상판 개수에 대한 분석은 이루어지지 않았다[9,10]. 위상판의 개수가 증가함에 비례해 시뮬레이션 내에서 자유공간 전파와 각 효과의 위상판 생성 횟수가 증가하는 만큼, 전체 시뮬레이션에 소요되는 시간은 선형적으로 증가한다.

가장 효과적으로 전체 시뮬레이션에 소요되는 시간을 줄이는 방법은 시뮬레이션 실행에 가장 많은 시간이 필요한 위상판 생성 및 자유 공간 전파의 횟수를 줄이는 것이다. 이는 사용하는 위상판의 개수를 줄임으로써 가능하다. 이때 핵심은 적은 개수의 위상판을 사용하더라도 빔이 전파되는 동안 난류와 열적 블루밍으로 인해 변화하는 특성이 많은 개수의 위상판을 사용하였을 때와 같아야 한다는 것이다. 따라서 본 논문에서는 열적 블루밍과 난류 현상을 모사하면서 동시에 변화되는 빔의 특성 또한 유지되는 최소 위상판의 개수를 제시하고, 시뮬레이션 시간을 크게 줄일 수 있는 방법을 보고한다. 비교를 위하여 충분히 많은 수의 위상판(150개)을 사용하여 난류와 열적 블루밍을 모사하고, 대조군으로 본 논문에서 제안하는 시뮬레이션은 전파 환경 및 레이저 전력에 맞춰 수량을 최소화한 위상판을 사용한다. 목표물에서의 빔의 성능은 조준점의 PIB를 통해 계산하며, 두 효과의 모사 여부를 확인하기 위하여 전파 경로 축상의 신틸레이션 지수의 평균 변화량을 측정한다. 시뮬레이션 비교 결과, 전파 거리가 5 km, 굴절률 구조 상수 함수가 10−14 m−2/3 이하, 왜곡 지수(distortion number)가 220 rad 이하일 때, 기존의 방법은 시뮬레이션 10,000회 반복에 18일이 소요된 반면, 제안된 방법은 4일이 소요되어 시뮬레이션 시간을 약 4.5배 단축할 수 있었다.

II. 대기 난류 채널에서의 열적 블루밍이 고려된 빔 전파 이론

2.1. 확률적 파동 방정식과 split step 방법

대기를 통과하는 레이저 빔의 전파는 확률적 파동 방정식을 통해 모사할 수 있다. z축 방향으로 전파하는 전기장 E에 대한 확률적 파동 방정식은 근축 근사를 적용하여 식 (1)과 같이 표현된다[11].

2ikEz+2E+2k2nRE=0

여기에서 ▽2 = 2/∂x2+2/∂y2은 라플라스 연산자, k는 파수, 그리고 n(R)은 공간 위치 R에 대한 굴절률을 의미한다. 상기 방정식은 통상 split step 방법으로 계산되며, 대기 난류에 의한 공간적 굴절률 변화는 위상판으로 표현할 수 있다. 위상판은 특정 거리에 해당하는 굴절률 변화를 공간적으로 표현한 매우 얇은 두께의 판으로서 전파 거리에 걸쳐 일정 거리마다 배치할 수 있다. 위상판과 위상판 사이는 자유 공간이다. 따라서 식 (1)을 공간 주파수 영역과 공간 영역으로 나누어 계산하면 대기를 전파하는 레이저 빔의 계산이 가능하다.

그림 1은 고출력 레이저 무기의 레이저 빔이 대기 채널을 통과하여 목표물로 조사되는 과정과 이를 시뮬레이션으로 구현하는 방법을 묘사한 그림이다. 고출력 레이저 빔은 대기 난류를 통과하여 목표물에 도달하는데, 이때 레이저 빔은 대기를 전파하면서 대기 난류에 의한 여러 효과(신틸레이션, 빔 원더링 등)를 겪을 뿐 아니라 고출력에 의한 열적 블루밍도 발생시킨다. 이를 시뮬레이션으로 모사하기 위하여 대기 난류에 의한 공간적 굴절률 변화, 그리고 열적 블루밍을 위상판으로 표현한다. Split step 방법은 전파 거리에 이러한 위상판을 주기적으로 위치시켜 레이저 빔의 전파를 계산하는 방법이다. 그림 1에 도시한 바와 같이 구간 ①은 난류 위상판 구간으로 난류로 인한 위상 왜곡이 반영되는 구간이며, 뒤이어 열적 블루밍으로 인한 위상 왜곡이 구간 ②에서 반영된다. 위상판과 위상판 사이의 구간 ③은 자유 공간 전파 구간으로, 빛의 회절 효과와 해당 거리만큼의 대기 흡수 및 산란으로 인한 손실이 반영된다.

Figure 1. Simulation of beam propagation over an atmospheric channel for high-energy laser beam.

2.2. 자유 공간 및 위상판을 통과하는 빔의 전파

위상판과 위상판 사이의 자유 공간에는 굴절률의 변화가 존재하지 않는다. 따라서 식 (1)에서 n(R) = 0이 되므로 아래 식만 남게 된다.

2ikEz+2E=0

식 (2)는 공간 주파수 영역으로, Huygens-Fresnel 적분을 통해 빔의 회절 효과를 계산할 수 있다. Huygens-Fresnel 적분은 식 (3)과 같이 고속 푸리에 변환을 통해 계산할 수 있다.

Ex,y,z+Δz=F1FEx,y,zexpiπλzkx2+ky2+ikΔz

여기에서 ∆z는 자유 공간 전파 구간의 거리, FF−1은 푸리에 변환과 푸리에 역변환, λ는 파장, kxky는 각각 xy방향의 공간 주파수를 의미한다.

자유 공간 전파 이후 빔은 위상판에 도달한다. Split step 방법에서 위상판은 z축 방향의 작은 굴절률 변화를 설명하기 위해 아주 얇게 생성되며 따라서 식 (1)의 두번째 항 ▽2E가 생략되어 아래 식만 남는다.

2ikEz+2k2nRE=0

식 (4)를 통해 난류와 열적 블루밍에 의한 빔의 위상 변화가 계산된다. 자유 공간 전파 및 위상판을 통과하는 과정을 전파 경로 시작 지점부터 마지막 지점까지 순차적으로 계산을 수행하면 전체 전파 경로 동안 레이저 빔이 받는 광학 효과들이 반영된다.

2.3. 흡수 및 산란을 통한 대기 손실

레이저 빔은 대기 중의 고유한 단면적을 가지는 여러 기체 분자 및 에어로졸에 의해 광자 일부가 흡수되거나 산란된다. 반응하는 물질의 입자 종류 및 크기에 따라 레일리 산란(Rayleigh scattering), 미에 산란(Mie scattering), 비선택적 산란(non-selective scattering) 등으로 나뉘며 이 현상으로 인해 목표물로 향하는 레이저 빔의 에너지가 감소한다. 흡수되거나 산란되는 에너지의 양은 일반적으로 파장에 따른 흡수 및 산란 계수로 나타내어지며 이로 인한 레이저 빔의 투과도(transmittance) τ를 식 (5)와 같이 계산할 수 있다.

τ(L)=exp0L α(z)+σ(z)dz

여기에서 ασ는 각각 파장에 따른 흡수 및 산란 계수를, L은 전파 거리를 의미한다. 흡수 및 산란 계수 값은 모두 빛의 파장에 의존한다. 이 계수 값은 이론적으로 도출하는 것이 까다로울 뿐만 아니라, 이론에 포함된 파라미터를 정확하게 측정하는 것이 현실적으로 거의 불가능하므로 통상 가시 거리를 활용하여 해당 파장에 대한 흡수 및 산란 계수를 구해서 사용한다. 그 대표적인 예로 Kim 모델이 있다[12].

2.4. 대기 난류 위상판 생성 방법

대기 난류의 특성을 설명하는 난류의 에너지 캐스케이드 이론(energy cascade theory of turbulence)에 의하면 대기의 풍속은 레이놀즈 수가 임계점을 초과할 때까지 증가한다. 이러한 작용으로 외부 규모(L0) 크기의 불안정한 공기 덩어리인 와류가 형성되며, 형성된 와류는 관성에 의해 더 작은 내부 규모(l0)의 와류들로 분해되면서 받아들인 에너지를 손실 없이 전달한다. 이때 L0의 크기보다 작은 와류는 통계적으로 균질하고 등방성인 것으로 가정한다. 와류의 에너지 전달 과정이 L0 크기의 와류에서 l0 크기의 와류까지 진행되면 받아들인 에너지는 열로 소실되며 최종적으로 생성된 와류는 소멸한다. 와류들은 모두 서로 다른 온도와 압력을 가지는 셀(cell)로 구성되어 있으며, 각 셀의 굴절률 또한 모두 다르다. 이러한 셀로 구성된 와류가 분포된 대기의 굴절률은 시공간적으로 계속 변화하며, 이를 광학 난류라고 한다[11].

시뮬레이션에서 구현하는 광학 난류는 통계적으로 균질하고 등방성을 가진다고 가정되었으며, 생성되는 모든 난류의 형태는 무작위 굴절률 분포를 띄어야 한다. 이는 공간 전력 스펙트럼으로 표현될 수 있으며, 광학 난류를 설명하는 데 대표적으로 사용되는 굴절률 공간 전력 스펙트럼은 난류의 내외부 규모를 고려한 모델인 수정된 본 카르만 모델(modified von Karman spectrum)을 따른다고 알려졌다[11].

Φnκ=0.033Cn2exp(κ2/κm2)κ2+κ0211/6,0κ<;κm=5.93/l0;κ0=2π/L0

식 (6)에서 Cn2는 굴절률 구조상수로서 대기 난류의 강도를 나타내는 값이고, κmκ0은 각각 스칼라 공간 파수로 와류의 내외부 규모를 고려한 상수이다. 또한 κ는 광축에 대하여 수직인 평면의 2차원 공간 파수이다. 식 (6)의 2차원 무작위 위상 분포를 가지는 판을 생성함으로써 대기 채널의 난류를 모사한 난류 위상판을 얻을 수 있다. 난류 위상판 생성 방법으로는 고속 푸리에 변환 방법, 공분산 행렬 방법, Zernike 다항식 방법, 저조파 방법 등이 있으나 본 논문에서는 가장 간단하고 빠르게 계산할 수 있는 고속 푸리에 변환 방법을 사용하였다[13]. 고속 푸리에 변환 방법으로 샘플링하여 구한 위상판의 위치 (x, y)에서의 위상은 아래와 같다.

ΦT(x,y)= κx κys κ x, κ yFϕ κ x, κ yexpi κ xx+ κ yyΔκxΔκy

여기서 s(κx, κy)는 평균 0, 분산 1의 허미시안(Hermitian) 복소 가우시안 과정으로, 이 랜덤 함수를 통해 대기의 무작위성이 부여된다. Fϕ(κx, κy)은 굴절률 변동의 2차원 파워 스펙트럼을 의미한다. 위 과정으로 생성된 난류 위상판은 전체 전파 경로를 동일한 거리 구간으로 나눈 후 각 구간의 중심에 배치되며, 난류 위상판을 통과한 전기장 E'(R)은 다음과 같이 계산된다.

E'(x,y,z)=E(x,y,z)expiΦT(x,y)

2.5. 열적 블루밍 위상판 생성 방법

열적 블루밍은 고출력 레이저(HEL) 빔이 대기로 전파되면서 발생하는 열로 인한 비선형 광학 현상이다. 고출력 레이저 빔이 대기를 통과할 때 대기 중의 특정 분자와 미립자들이 레이저 에너지의 일부를 흡수하는데, 이때 대기에 흡수된 에너지는 주변 대기를 가열하여 팽창시키고, 고출력 레이저 빔 전파 경로상에 분산된 열적 렌즈(thermal lens)를 형성한다. 이렇게 형성된 열적 렌즈를 통과한 고출력 레이저 빔은 너비가 넓어지고, 전파 경로가 휘어진다. 이로 인해 최초 지점에서 완전한 가우시안 분포의 형태였던 레이저 빔은 목표 지점에서 초승달 모양의 형태를 가지게 된다.

열적 블루밍 현상은 제1 열역학 법칙의 에너지 보존 법칙을 바탕으로 내부 에너지만을 고려한 총 에너지 방정식을 굴절률 항으로 나타낸 식 (9)의 에너지 균형 방정식 풀이를 통해 계산할 수 있다[14].

nx,y,z,tt+vznx,y,z,t=1n0αzρCPT0Ix,y,z,t

여기서 t는 시간, v는 대기 이류(advection)의 속도 벡터, ρ는 대기의 밀도, CP는 정압 비열, n0는 공기의 굴절률, T0는 대기의 온도, I(R,t)는 시간 t에 따른 공간 좌표 상의 레이저 빔의 강도(intensity)를 의미한다.

식 (9)의 물리적인 의미는 고출력 레이저 빔이 전파하며 대기에 흡수되는 에너지와 횡방향으로 불어오는 대기 이류로 인해 냉각되는 에너지가 균형을 이루는 것이다. 식 (9)에서 레이저 빔에 영향을 주는 대기 이류 v는 레이저 빔이 전파되는 방향에 수직으로 불어온다고 볼 수 있다. 이때 v(z) = [v(z), 0, 0]으로 나타낼 수 있다. 또한 식 (9)로 설명되는 고출력 레이저 빔 주변 대기의 에너지 변화는 일정 시간 경과 시 레이저 빔의 강도 변화가 거의 없는 정상 상태(steady state)에 도달한다. 따라서 식 (9)의 시간 항을 모두 생략할 수 있으며 이를 굴절률 n(R)에 대한 식으로 정리하면 아래와 같이 거리에 따라 변화하는 굴절률의 양을 이산적으로 계산할 수 있다.

nx,y,z=1n0αzρCPT0vzx I ξ,y,zdξ

이렇게 구해진 열적 블루밍으로 인해 변화하는 굴절률의 양을 거리에 대해서 적분하면, 난류 위상판과 같이 일정 거리에 대하여 발생하는 위상 왜곡 양을 담은 위상판으로 나타낼 수 있으며 이는 식 (11)과 같이 나타낸다.

ΦBx,y,L=k0L 1n 0 αzτ(z) ρCP T 0 vzdzx I ξ,y,0dξ

열적 블루밍의 세기를 나타내는 지표로서 통상 Bradley & Herrmann 왜곡 지수(Nd)를 사용하며 식 (12)와 같이 표현한다[15,16].

Nd(L)=42kP0D00L n 01αzτ(z) ρC PT 0v(z) dz

여기서 P0는 빔의 최초 출력 세기를, D0는 빔의 최초 직경을 의미한다. Bradley & Herrmann 왜곡 지수를 사용하면 식 (11)은 아래와 같이 간단히 표현된다.

ΦBx,y,L=Nd(L)D042kP0x I ξ,y,0dξ

난류 위상판과 마찬가지로 위 과정으로 생성된 열적 블루밍 위상판은 동일한 거리로 나누어진 각 구간의 중심에 배치되며, 열적 블루밍 위상판을 통과하는 전기장 E'(R)을 계산하면 식 (14)와 같다.

E'(x,y,z)=E(x,y,z)expiΦB(x,y)

그림 2는 2.4절과 2.5절에서 설명한 위상판 생성 방법으로 생성한 위상판의 예시로, 그림 2(a)는 난류 위상판, 그림 2(b)는 열적 블루밍 위상판이다. 열적 블루밍 위상판 생성 간 대기 이류는 2 m/s로 설정하였다. 난류 위상판은 전체 대기에 서로 다른 크기의 와류가 통계적으로 균질하게 분포되어 있으며, 난류 위상판을 통해 난류로 인한 위상 왜곡량이 레이저 빔에 반영된다. 마찬가지로 열적 블루밍 위상판을 통해 대기 이류 방향으로 축적된 위상 왜곡량 또한 전파되는 레이저 빔에 반영된다. 이때 주목할 점은 식 (7)의 s(κx, κy) 함수로 인하여 난류 위상판에 무작위성이 부여되는 반면, 열적 블루밍 위상판은 레이저 전력에 따라 변화할 뿐 무작위성은 없다는 것이다.

Figure 2. Example of phase screens for (a) atmospheric turbulence and (b) thermal blooming.

III. Split step 방법을 활용한 시뮬레이션 및 위상판 개수 판단의 조건

Split step 방법을 활용한 레이저 빔 전파 시뮬레이션에서 레이저 빔은 2차원 공간 내 특정 크기의 판 안에 표현되며, 위상판 또한 이와 같은 크기로 생성된다. 이러한 특성으로 인해 시뮬레이션이 타당성을 가지기 위해서는 시뮬레이션에 사용된 주요 파라미터들이 특정 조건을 만족해야 한다. 또한 사용하는 각 효과의 위상판 개수에 따라 빔의 산란 및 왜곡되는 정도가 달라지므로, 적절한 개수의 위상판을 판단하는 것이 중요하다.

3.1. Split step 방법을 활용한 시뮬레이션의 조건

Split step 방법을 활용한 레이저 빔 전파 시뮬레이션은 널리 이용되는 방법이기 때문에 주요 파라미터 설정 기준과 조건은 잘 알려져 있다[7,8]. 우선 생성되는 위상판은 와류의 외부 규모와 내부 규모의 영향을 충분히 반영할 수 있어야 한다. 따라서 위상판 크기 A가 외부 규모보다 충분히 커야 하며, 이를 위해 A > L0를 만족해야 한다. 또한 위상판 격자의 크기 ∆x∆y는 내부 규모보다 충분히 작아야 하므로, 각각 ∆x < l0/2와 ∆y < l0/2를 만족해야 한다. 두번째로, 나이퀴스트 샘플링 정리(Nyquist’s sampling theorem)를 만족하기 위해 위상판 내의 이웃한 격자 간의 위상 차이는 π보다 작아야 한다. 세번째로, 빔의 전파 계산에 사용되는 푸리에 변환에서 신호가 왜곡되는 에일리어싱(aliasing) 현상을 방지하기 위해 위상판 간 거리는 기본적으로 ∆z < kA∆xπ를 만족해야 한다. 네번째로, split step 방법에서는 격자(grid)의 불연속적 특성으로 인하여 가장자리에서 빔이 산란되어 판을 벗어나는 경우 옆부분에 해당 빔의 에너지가 표현되는 가장자리 효과(edge effect)가 발생하는데, 이때 가장자리 효과를 최소화하기 위하여 A > 2θmax ∆z를 만족해야 한다. 여기서 θmax = max[√2/k × ∂f/∂x]로서 빔 전파간 빔의 산란 각도의 최대값이다. f는 난류 위상판의 특정 격자의 위상값을 의미한다.

3.2. 대기 난류를 모사하기 위하여 필요한 위상판 개수

레이저 빔이 대기 난류를 통과하면서 변동하는 강도는 전통적으로 약한 변동 이론과 강한 변동 이론으로 분류된다. 여기서 Rytov 분산 값은 변동의 강약을 구분하는 값이며 식 (15)와 같이 쓸 수 있다[11].

σR2=1.23Cn2k7/6L11/6

빔 전파 시작 지점부터 σR2 < 1을 만족하는 거리까지를 빔 변동 강도가 약한 구간(weak fluctuation)이라고 한다. 대기 난류를 통과하는 레이저 빔의 변동 강도는 신틸레이션 지수로 나타낼 수 있으며, σR2 ≥ 1일 때 신틸레이션 지수는 포화 상태에 이른다. 집속된 가우시안 빔에 대한 신틸레이션 지수의 이론식은 잘 알려져 있으며, 포화 상태에 이르기 전의 값은 식 (16)으로 구할 수 있다[11,17].

σIT20,0,L=4.42σR2Γ5/6 σ pe DB 2+3.86σR2Rei5/6 2F156 ,116 ;176 ;1Θ+iΓ1116Γ5/6,σR2<1

여기서 Г = Г0/(1 + Γ02), Г0 = 8L/kD02이며 가우시안 빔 매개변수를 의미한다. 2F1(a, b; c; z)은 가우스 초기하 함수(Gaussian hypergeometric function)이며 Θ = 1/(1 + Γ02)은 빔 곡률 매개변수를 의미한다. DB는 아무런 왜곡 없이 목표 거리에 도달한 이상적인 가우시안 빔의 회절 한계를 고려한 빔의 직경이다. σpe2는 지터(jitter)로 인한 포인팅 오류 분산값으로 식 (17)과 같다.

σpe2~ λLD02 D0r05/3,D0r0<<1 λLD02 r0D01/3,D0r0<<1r0= 0.423Cn 2 k 2L 3/5

여기서 r0는 프라이드 파라미터(Fried’s parameter)로 대기 난류의 영향 정도를 나타내는 척도이다. 또한 신틸레이션 지수는 아래와 같이 시뮬레이션을 통해 생성된 거리별 전파 경로상 중심에서의 빛의 강도 표본값을 바탕으로 계산할 수 있다.

σI2x,y,L=I2x,y,LIx,y,L21

여기서 < >는 앙상블 평균값을 의미한다. σR2 < 1일 때, 해당 구간에서 신틸레이션 지수의 변화 정도를 제한함으로써 레이저 빔의 대기 난류 전파를 정확하게 모사할 수 있다. 이때 기준은 다음과 같다. 우선 난류 위상판 1개를 통과할 때 증가하는 신틸레이션 지수는 각 위상판이 충분히 작은 위상 변동만 일으켜야 한다는 조건을 만족시키기 위하여 0.1을 초과할 수 없다[8]. 이 기준은 대기 전파 시뮬레이션 연구에서 종합적으로 경험에 의해 얻어진 수치로서 널리 사용되고 있으며 이를 만족하는 위상판 개수(NT1)는 아래와 같이 나타낼 수 있다.

NT1=10LσIT20,0,1.23Cn2k7/66/111.23Cn2k7/66/11

또한 전체 거리를 전파하며 변화하는 신틸레이션 지수가 작더라도, 충분한 개수의 위상판을 사용하기 위한 조건을 만족시키기 위하여 난류 위상판 1개를 통과할 때 변화하는 신틸레이션 지수의 변화량은 전체 거리를 전파하는 동안 변화하는 신틸레이션 지수의 총 변화량의 10% 이하여야만 한다[8,18]. 위상판과 위상판의 거리는 모두 동일하므로, 각 위상판이 변화시키는 위상 변화의 총량 또한 동일하다. 따라서 위 조건을 만족하는 위상판 개수(NT2)는 아래와 같다.

NT2=10

위의 두 조건을 모두 만족하는 경우 전체 경로에 걸쳐 대기 난류를 모사하는 데 필요한 위상판의 최소 개수(NT)를 구할 수 있으며, 아래와 같이 표현할 수 있다.

NT=maxNT1,10

3.3. 열적 블루밍을 모사하기 위하여 필요한 위상판 개수

열적 블루밍의 경우 전파 거리 L에 대하여 고출력 레이저 빔의 중심에 주는 왜곡량의 최댓값은 아래와 같이 표현된다[14].

ΦB_max=12πNdL

또한 총 M개의 열적 블루밍 위상판을 사용하면 m번째 열적 블루밍 위상판이 가지는 왜곡의 양은 식 (23)과 같다.

ΦBx,y,mLM=α(z)k1n0eα(z)+σ(z)CPT0ρv(z)α(z)+σ(z)emL/Mem1L/Mx I ξ,y,0dξ

본 논문에서는 열적 블루밍에 필요한 위상판 개수를 찾기 위해 고출력 레이저 빔의 성능을 측정할 때 주로 사용되는 power in the bucket (PIB) 지표를 활용하였다. PIB는 빔 전파 최초 지점에서 지름이 D0인 원형 버킷 안에 있는 에너지의 총량과, 전파된 후 목표 거리의 중앙에서의 회절 한계를 고려한 지름 DB의 원형 버킷 안에 들어간 에너지의 총량을 측정하여 식 (24)와 같이 계산한다.

PIB= DB/2 DB/2Ix,y,Ldxdy D0/2 D0/2Ix,y,0dxdy

그림 3은 열적 블루밍 세기와 위상판 개수를 변화시켜가며 측정한 PIB 값을 표준화한 결과이다. 그래프에 표시된 흰 점선은 모든 열적 블루밍 세기 범위에서 표준화한 PIB 값이 1로 수렴하는 데 필요한 최소 개수의 위상판 개수를 나타낸 것이다. 그림 3에 따르면 위상판의 개수에 따라 PIB 값의 차이가 큰 것을 볼 수 있는데, 이는 위상판 간 거리가 멀어서 자유 공간 전파 중에 빔이 많이 분산되기 때문이다. 반면 충분한 개수의 위상판이 사용되면 표준화한 PIB 값이 1로 수렴하는데, 이 지점에서 위상판이 가지는 위상 왜곡량의 비율을 식 (22)와 식 (23)을 통해 측정함으로써 위상판 한 개가 가져야 할 위상 왜곡량이 총 왜곡량의 5% 미만임을 확인하였다. 이 조건을 만족하는 위상판 개수(NB1)는 식 (25)와 같이 나타낼 수 있다.

Figure 3. Normalized power in the bucket value as functions of distortion number, Nd, and number of thermal blooming screens NB.

NB1=52πD0P0eαz+σ(z)LΔzD0 I ξ,0,0dξ

또한 나이퀴스트 샘플링 정리(Nyquist sampling theorem)를 만족하기 위해 열적 블루밍 위상판의 위상차는 π를 넘을 수 없다[7]. 이를 만족하는 위상판 개수(NB2)는 아래와 같이 나타낼 수 있다.

NB2=kn01α(z)LCPT0ρv(z)πeα(z)+σ(z)ΔzD0 I ξ,0,0dξ

앞선 두 조건을 모두 만족시키면 전체 경로에 걸쳐 열적 블루밍을 모사하는 데 필요한 위상판의 최소 개수(NB)를 구할 수 있으며, 이는 식 (27)과 같이 표현할 수 있다.

NB=maxNB1,NB2

3.4. 대기 난류와 열적 블루밍이 함께 고려된 시뮬레이션에 필요한 위상판 개수

지금까지 대기 난류와 열적 블루밍이 개별적으로 고려된 시뮬레이션에 필요한 위상판 개수를 구하였다. 그러나 실제 대기 공간에서는 고출력 레이저 빔 전파 시 난류와 열적 블루밍이 동시에 발생한다. 본 논문에서는 계산의 용이성과 시뮬레이션의 단순화를 위하여 그림 1과 같이 난류 위상판과 열적 블루밍 위상판을 전파 경로상에 함께 배치하였다. 또한 이때 사용되는 각 효과의 위상판 개수는 시뮬레이션 조건(난류 세기, 열적 블루밍 세기 등)에 따라 달라지는 위상판 개수 중 가장 큰 값(Nps)을 따르는 것으로 설정하였다. 이를 수식으로 표현하면 식 (28)과 같다.

Nps=maxNT,NB=maxNT1,10,NB1,NB2

또한 난류 위상판과 열적 블루밍 위상판을 통과하는 전기장 E'(R)의 계산 결과는 아래와 같다.

E'(x,y,z)=E(x,y,z)expiΦT(x,y)expiΦB(x,y)

IV. 시뮬레이션 결과

이번 장에서 우리는 150개의 위상판을 사용하는 기존 방법과 본 논문에서 제안하는 방법의 시뮬레이션 정확도 및 소요 시간을 비교하였다. 시뮬레이션에 사용된 주요 매개변수의 값은 표 1에 요약되어 있다. 먼저, 시뮬레이션은 3장에서 제시한 기준을 모두 만족하는 조건에서 수행하였다. 고출력 레이저 빔은 집속된 가우시안 빔을 사용한다고 가정하며 전송 거리는 5 km이다. 난류 위상판은 식 (5)의 수정된 본 카르만 모델을 사용하여 생성하고, 와류의 내외부 규모는 각각 4 mm, 1.7 m로 설정하였다. 위상판의 크기는 4 × 4 m2이고, 이를 4,096 × 4,096개의 격자 포인트 해상도로 구현하였다. 또한 대기 이류의 속도는 한국 45개 지역에서 측정된 연평균 풍속인 2 m/s로 설정하였으며, 난류 세기(굴절률 구조상수) 또한 한국 도시별 굴절률 구조상수의 값을 고려하여 10−16–10−13 m−2/3으로 설정하였다[19].

Table 1 . Simulation parameters.

ParameterValue
Simulation Grid4,096 × 4,096
Side Length of Phase Screen (m)4
Wavelength (nm)1,070
Initial Beam Diameter (m)0.5
Propagation Distance (km)5
Transmitter Beam Power (kW)1–1,000
Absorption Coefficient (m−1)5 × 10−6
Scattering Coefficient (m−1)5 × 10−5
Refractive Index of Air1.0003
Specific Heat at Constant Pressure (J/kg/K)1,004
Density of Air at Constant Pressure (kg/m3)1.293
Initial Air Temperature (K)300
Refractive Index Structure Constant (m−2/3)10−13–10−16
Advection Wind Speed (Direction) (m/s)2 (Left to Right)
Beam Propagation DirectionHorizontal Direction


두 시뮬레이션의 정확성을 비교하기 위하여 난류와 열적 블루밍 세기별로 총 10,000회의 시뮬레이션을 시행하였다. 이는 무작위 특성의 난류와 열적 블루밍이 빔에 끼치는 영향을 통계학적으로 분석하기 위함이다. 이렇게 획득한 10,000개의 표본에 대해 PIB를 식 (24)를 사용하여 계산하였다.

다음으로 PIB 값 분포를 통해 구해진 PMF의 유사도를 결정계수(R-squared)를 측정하여 비교하였다.

R2=1 i=1Mgihi2 i=1Mgig¯2

여기서 gi는 기존 방법으로 구한 PIB 분포 그래프에서 각 bin의 발생 확률이며 gi의 평균값을, hi는 제안된 방법으로 구한 PIB 분포 그래프의 각 bin의 발생 확률을 의미한다. 결정계수 값이 0.9 이상이면 통상 통계적 유사도가 높은 것으로 볼 수 있다[20,21]. 반면 결정계수 값이 낮으면 오차가 많으므로 두 분포의 유사도가 낮다고 본다.

또한 전파되는 동안 거리별 전파 경로 축상 신틸레이션 지수(on-axis σI2)를 식 (18)을 사용하여 계산하고, 두 시뮬레이션의 on-axis σI2의 평균 오차(η)를 다음과 같이 구하였다.

η= i=1MσIc2(0,0,zi)σIp2(0,0,zi)M

여기서 σIc2σIp2는 각각 기존 방법과 제안된 방법으로 계산된 신틸레이션 지수이다.

그림 4는 모두 열적 블루밍의 강도(Nd)가 50 rad인 조건에서 난류 세기를 달리했을 때 목표 거리에서의 고출력 레이저 빔 형태의 예시이다. 화면 중앙의 흰색 원은 PIB 값을 측정하기 위한 3 cm 지름의 원형 버킷이다. 그림 4(a)는 난류 세기가 약한 대기 채널(Cn2 = 10−16 m−2/3)을 통과했을 때의 모습이다. 난류에 의한 신틸레이션과 빔 원더링 현상이 또렷이 확인되며, 열적 블루밍에 의해 빔의 형태가 대기 이류의 역방향으로 초승달처럼 휜 것을 확인할 수 있다. 그림 4(b)는 난류 세기가 강할 때(Cn2 = 10−13 m−2/3) 대기를 통과한 빔의 형태이다. 그림 4(a)보다 신텔레이션과 빔 원더링 현상이 강하게 발현되었으나 열적 블루밍은 상대적으로 미미한 것을 알 수 있다. 이는 강한 난류로 인해 빔이 많이 퍼져 있기 때문이다.

Figure 4. Exemplary beam intensity profiles after 5-km propagation over an atmospheric channel when the distortion number is 50 rad: (a) weak turbulence (Cn2 = 10−16 m−2/3) and (b) strong turbulence (Cn2 = 10−13 m−2/3). The white circle represents a bucket having a 3-cm diameter.

그림 5는 기존 방법 대비 제안된 방법의 PIB 분포별 PMF 결과로, 많은 횟수의 시뮬레이션을 시행했을 때 두 방법 간의 확률적인 차이를 비교할 수 있다. 이때 PMF의 bin 개수는 결정계수 값 비교를 위하여 난류 및 열적 블루밍 세기와 상관없이 50개로 설정하였다. 그림 5(a)는 난류 세기가 작고(Cn2 = 10−15 m−2/3), 열적 블루밍의 세기는 126 rad으로, 이때 결정계수는 0.9749이다. 그림 5(b)는 난류 세기가 강하고(Cn2 = 10−13 m−2/3), 열적 블루밍의 세기가 84 rad인 경우이며, 이때의 결정계수는 0.9906이다. 따라서 두 경우 모두 기존의 방법과 제안된 방법의 PIB 분포별 PMF 히스토그램의 결과가 잘 일치함을 알 수 있다.

Figure 5. Probability mass functions of power in the bucket values under (a) weak turbulence (Cn2 = 10−15 m−2/3) when the distortion number is 126 rad and (b) strong turbulence (Cn2 = 10−13 m−2/3) when the distortion number is 84 rad. The total number of simulations is 10,000.

표 2는 조건을 다르게 하여 측정한 결정계수와 on-axis σI2 의 평균 오차값의 결과이다. 표 2에서 레이저 출력 밀도는 레이저 출력과 빔 단면적으로 계산된 평균값을 사용하였다. 그림 6(a)표 2의 결정계수 결과를, 그림 6(b)표 2의 on-axis σI2 의 평균 오차값의 결과를 레이저 출력 밀도에 대한 그래프로 나타낸 것이다. 결과를 보면 열적 블루밍 세기(레이저 출력)가 커질수록, 난류 세기는 작아질수록, 결정계수 값이 점차 감소하는 것을 볼 수 있다. 이는 난류 세기가 약할수록 열적 블루밍이 상대적으로 빔의 요동 및 형태 변화에 더 큰 영향을 미치기 때문이다. 열적 블루밍의 강도가 일정량 이상이 되면, 난류와 열적 블루밍의 상호작용에 의한 섭동(thermal Rayleigh scattering)이 추가되어 본 논문에서 제시한 기준보다 더 많은 개수의 위상판이 필요하다[22]. 반면 난류 세기가 강하면 난류가 열적 블루밍보다 상대적으로 빔의 변화에 더 큰 영향을 미치며, 빔의 섭동이 열적 블루밍에 의해 크게 영향을 받지 않는다. on-axis σI2 의 평균 오차값은 열적 블루밍 세기가 210 rad (레이저 출력 밀도 2.5 × 106 W/m2)를 초과할 경우 값이 현저하게 증가함을 알 수 있다. 따라서 기존 방법의 시뮬레이션을 기준으로 할 때, 우리가 제시한 위상판 개수는 모든 난류 세기 강도(Cn2 = 10−16 – 10−13 m−2/3)에서 레이저 출력 밀도가 2.5 × 106 W/m2 이하일 때 사용하기 적합하다.

Figure 6. (a) R-squared and (b) average error of on-axis scintillation indexes for various atmospheric conditions and laser power densities.

Table 2 . R-squared and average error of on-axis scintillation indexes (η) for various atmospheric conditions and laser output powers.

P0(Nd)1 kW
(0.4 rad)
10 kW
(4 rad)
50 kW
(21 rad)
100 kW
(42 rad)
200 kW
(42 rad)
300 kW
(126 rad)
500 kW
(210 rad)
600 kW
(252 rad)
700 kW
(295 rad)
800 kW
(337 rad)
900 kW
(379 rad)
1 MW
(420 rad)
Cn210−16
(m−2/3)
R20.98590.97750.99980.95550.9428-0.90720.9001---0.8064
η0.00210.00200.00590.00820.0084-0.06500.0830---0.1983
10−15
(m−2/3)
R20.99000.98370.9723-0.97580.97490.95370.94160.9325---
η0.01240.00970.0090-0.00850.00800.09520.14250.1948---
10−14
(m−2/3)
R20.98360.99440.98530.9927--0.96370.9455-0.8923--
η0.05000.07260.08250.0708--0.08820.1130-0.1783--
10−13
(m−2/3)
R20.99250.9962--0.9906-0.98050.9638--0.9617-
η0.06500.0860--0.0689-0.08250.0868--0.0990-
Laser Power Density (kW/m2)5.15.1 × 102.5 × 1025.1 × 1021.0 × 1031.5 × 1032.5 × 1033.1 × 1033.6 × 1034.1 × 1034.6 × 1035.1 × 103


대기 난류와 열적 블루밍이 함께 고려된 레이저 빔 전파 시뮬레이션의 과정은 크게 난류 위상판 및 열적 블루밍 위상판의 생성, 자유 공간 전파, 위상판 통과 및 기타 부수 과정(PIB 값 계산, 난류 세기 및 열적 블루밍 세기 연산 과정 등)으로 구분할 수 있다. 따라서 시뮬레이션 총 진행 횟수(N)에 따른 총 소요 시간(Ttotal)은 아래와 같이 계산할 수 있다.

Ttotal=Nps×Ttb+Ttul+Tprop+Tpass+Λ×N

여기서 TtbTtul은 각각 열적 블루밍 및 난류 위상판을 생성하는 시간, Tprop은 자유 공간 전파 구현시간, Tpass는 각 위상판 통과 연산 시간, Λ는 기타 부수 과정에 소요되는 시간을 의미한다. 대부분의 시간은 TtbTtul에서 소요된다.

식 (32)에서 알 수 있듯이 위상판의 개수가 증가할수록 시뮬레이션 내에서 자유공간 전파 시간과 각 효과의 위상판 생성 횟수도 그와 비례하여 증가하며, 전체 시뮬레이션에 소요되는 시간 또한 선형적으로 증가하게 된다. 이를 바탕으로 기존의 방법과 본 논문에서 제안한 방법을 사용한 두 시뮬레이션의 소요 시간을 비교하였으며 이를 그림 7에 나타내었다. 두 방법 모두 클락 속도 4.7 GHz의 Intel® CoreTM i7-8,700K 프로세서를 갖춘 개인용 컴퓨터에서 실행되었다.

Figure 7. (a) Total simulation time measured by using a personal computer versus the total number of simulations, N, under strong turbulence (Cn2 = 10−13 m−2/3) and weak turbulence (Cn2 = 10−16 m−2/3). (b) Total simulation time versus the refractive index structure constant, Cn2. The total number of simulations is 10,000.

그림 7(a)는 시뮬레이션 총 진행 횟수(N)에 따른 소요 시간을 나타낸다. 기존의 방법은 모든 난류 세기에 대해서 동일한 시간이 소요되었으나, 제안된 방법의 경우 난류 세기에 따라 소요 시간이 다르므로 약한 난류 세기(Cn2 = 10−16 m−2/3)일 때와 강한 난류 세기(Cn2 = 10−13 m−2/3)일 때의 소요 시간을 따로 나타내었다. 시뮬레이션 결과, 제안된 방법을 사용할 시 기존의 방법에 비하여 총 소요 시간이 최대 6배 감소하였다. 예를 들면 시뮬레이션 총 실행 횟수가 100,000회일 때, 제안된 방법은 난류 세기가 강할 경우 약 200일이, 난류 세기가 약할 경우에는 약 54일이 걸리지만, 기존의 방법은 모든 난류 조건에서 330일이 걸린다. 이는 제안된 방법에서 필요 위상판 개수에 따라 위상판 생성 및 자유공간 전파 횟수가 현저하게 줄어들기 때문이다.

그림 7(b)는 난류의 세기인 굴절률 구조상수 값을 달리해 가며 시뮬레이션 10,000회를 시행하는 시간을 나타낸다. 제안된 방법은 난류의 세기가 약할수록 필요한 위상판 개수가 적어지므로 더 적은 시간이 소요되며, 모든 난류 세기를 기준으로 보았을 때에조차 기존의 방법보다 소요 시간이 적은 것을 알 수 있다. 결과적으로 제안된 방법은 우리나라의 도시별 평균 최대 난류 세기(Cn2 = 10−13 m−2/3) 이하의 범위에서 실행 횟수가 많을수록, 그리고 난류 세기가 약할수록 기존 방법보다 시간 면에서 효율적이다.

V. 결 론

본 논문에서는 대기 난류와 열적 블루밍을 겪는 고출력 레이저 빔의 대기 전파를 빠르고 정확하게 시뮬레이션 수행하기 위하여 요구되는 위상판의 개수를 제시하였다. Split step 방법을 사용한 빔 전파 시뮬레이션의 시간 대부분은 자유공간 빔 전파와 위상판 생성에 소요되기 때문에, 시뮬레이션에 사용되는 위상판의 총 개수를 줄임으로써 시뮬레이션의 시간을 크게 단축할 수 있다. 먼저 전파 거리, 대기 난류의 세기, 레이저 빔의 출력 전력 등에 따라서 각 위상판이 레이저 신호에 미치는 적절한 위상의 양을 결정함으로써 정확한 시뮬레이션 구현에 요구되는 최소 위상판의 개수를 도출하였다. 본 논문에서는 레이저 출력 밀도 2.5 × 106 W/m2 (빔 직경 50 cm 레이저의 경우 출력 전력 500 kW) 미만의 빔이 대기를 전파할 때 이를 모사하기 위하여 필요한 최소 위상판의 개수를 제시하였으므로, 이 방법은 고출력 레이저 빔의 대기 전파를 통계적으로 분석할 때 유용하게 사용될 것으로 기대된다.

재정지원

BK21 Four 사업.

이해상충

저자는 본 논문과 관련된 어떠한 이해충돌 사항도 없었음을 밝힌다.

데이터 가용성

본 연구의 결과 분석 및 생성된 데이터는 모두 본 논문 내 명시되어 있으며 공공의 이용이 가능하다. 데이터에 접근하거나 사용하고자 하는 이는 저자에게 타당한 이유를 밝히고 허가를 득해 사용 가능하다.

Fig 1.

Figure 1.Simulation of beam propagation over an atmospheric channel for high-energy laser beam.
Korean Journal of Optics and Photonics 2024; 35: 49-60https://doi.org/10.3807/KJOP.2024.35.2.049

Fig 2.

Figure 2.Example of phase screens for (a) atmospheric turbulence and (b) thermal blooming.
Korean Journal of Optics and Photonics 2024; 35: 49-60https://doi.org/10.3807/KJOP.2024.35.2.049

Fig 3.

Figure 3.Normalized power in the bucket value as functions of distortion number, Nd, and number of thermal blooming screens NB.
Korean Journal of Optics and Photonics 2024; 35: 49-60https://doi.org/10.3807/KJOP.2024.35.2.049

Fig 4.

Figure 4.Exemplary beam intensity profiles after 5-km propagation over an atmospheric channel when the distortion number is 50 rad: (a) weak turbulence (Cn2 = 10−16 m−2/3) and (b) strong turbulence (Cn2 = 10−13 m−2/3). The white circle represents a bucket having a 3-cm diameter.
Korean Journal of Optics and Photonics 2024; 35: 49-60https://doi.org/10.3807/KJOP.2024.35.2.049

Fig 5.

Figure 5.Probability mass functions of power in the bucket values under (a) weak turbulence (Cn2 = 10−15 m−2/3) when the distortion number is 126 rad and (b) strong turbulence (Cn2 = 10−13 m−2/3) when the distortion number is 84 rad. The total number of simulations is 10,000.
Korean Journal of Optics and Photonics 2024; 35: 49-60https://doi.org/10.3807/KJOP.2024.35.2.049

Fig 6.

Figure 6.(a) R-squared and (b) average error of on-axis scintillation indexes for various atmospheric conditions and laser power densities.
Korean Journal of Optics and Photonics 2024; 35: 49-60https://doi.org/10.3807/KJOP.2024.35.2.049

Fig 7.

Figure 7.(a) Total simulation time measured by using a personal computer versus the total number of simulations, N, under strong turbulence (Cn2 = 10−13 m−2/3) and weak turbulence (Cn2 = 10−16 m−2/3). (b) Total simulation time versus the refractive index structure constant, Cn2. The total number of simulations is 10,000.
Korean Journal of Optics and Photonics 2024; 35: 49-60https://doi.org/10.3807/KJOP.2024.35.2.049

Table 1 Simulation parameters

ParameterValue
Simulation Grid4,096 × 4,096
Side Length of Phase Screen (m)4
Wavelength (nm)1,070
Initial Beam Diameter (m)0.5
Propagation Distance (km)5
Transmitter Beam Power (kW)1–1,000
Absorption Coefficient (m−1)5 × 10−6
Scattering Coefficient (m−1)5 × 10−5
Refractive Index of Air1.0003
Specific Heat at Constant Pressure (J/kg/K)1,004
Density of Air at Constant Pressure (kg/m3)1.293
Initial Air Temperature (K)300
Refractive Index Structure Constant (m−2/3)10−13–10−16
Advection Wind Speed (Direction) (m/s)2 (Left to Right)
Beam Propagation DirectionHorizontal Direction

Table 2 R-squared and average error of on-axis scintillation indexes (η) for various atmospheric conditions and laser output powers

P0(Nd)1 kW
(0.4 rad)
10 kW
(4 rad)
50 kW
(21 rad)
100 kW
(42 rad)
200 kW
(42 rad)
300 kW
(126 rad)
500 kW
(210 rad)
600 kW
(252 rad)
700 kW
(295 rad)
800 kW
(337 rad)
900 kW
(379 rad)
1 MW
(420 rad)
Cn210−16
(m−2/3)
R20.98590.97750.99980.95550.9428-0.90720.9001---0.8064
η0.00210.00200.00590.00820.0084-0.06500.0830---0.1983
10−15
(m−2/3)
R20.99000.98370.9723-0.97580.97490.95370.94160.9325---
η0.01240.00970.0090-0.00850.00800.09520.14250.1948---
10−14
(m−2/3)
R20.98360.99440.98530.9927--0.96370.9455-0.8923--
η0.05000.07260.08250.0708--0.08820.1130-0.1783--
10−13
(m−2/3)
R20.99250.9962--0.9906-0.98050.9638--0.9617-
η0.06500.0860--0.0689-0.08250.0868--0.0990-
Laser Power Density (kW/m2)5.15.1 × 102.5 × 1025.1 × 1021.0 × 1031.5 × 1032.5 × 1033.1 × 1033.6 × 1034.1 × 1034.6 × 1035.1 × 103

References

  1. P. E. Ross, "Economics drives a ray-gun resurgence: Lasers, cheaper by the shot, should work well against drones and cruise missiles," IEEE Spectrum 60, 40-41 (2023).
    CrossRef
  2. B. Zohuri, Directed Energy Weapons (Springer Charm, USA, 2016), 1st ed..
    CrossRef
  3. P. E. Nielsen, Effects of Directed Energy Weapons (CreateSpace Independent Publishing, USA, 2012), 1st ed..
  4. W. Moon and H. Kim, "Standard deviation of fiber-coupling efficiency for free-space optical communication through atmospheric turbulence," IEEE Photonics J. 15, 7302507 (2023).
    CrossRef
  5. R. Frehlich, "Simulation of laser propagation in a turbulent atmosphere," Appl. Opt. 39, 393-397 (2020).
    Pubmed CrossRef
  6. Z. Chen, D. Zhang, C. Xiao, and M. Qin, "Precision analysis of turbulence phase screens and their influence on the simulation of Gaussian beam propagation in turbulent atmosphere," Appl. Opt. 59, 3726-3735 (2020).
    Pubmed CrossRef
  7. D. L. Knepp, "Multiple phase-screen calculation of the temporal behavior of stochastic waves," Proc. IEEE 71, 722-737 (1983).
    CrossRef
  8. J. M. Martin and S. M. Flatté, "Intensity images and statistics from numerical simulation of wave propagation in 3-D random media," Appl. Opt. 27, 2111-2126 (1988).
    Pubmed CrossRef
  9. M. F. Spencer, "Wave-optics investigation of turbulence thermal blooming interaction: II. Using time-dependent simulations," Opt. Eng. 59, 081805 (2020).
    CrossRef
  10. C. E. Murphy and M. F. Spencer, "Investigation of turbulence thermal blooming interaction using the split-step beam propagation method," Proc. SPIE 10772, 1077208 (2018).
  11. L. C. Andrews and R. L. Phillips, Laser Beam Propagation Through Random Media, 2nd ed. (SPIE Press, USA, 2005).
    CrossRef
  12. I. I. Kim, B. McArthur, and E. J. Korevaar, "Comparison of laser beam propagation at 785 nm and 1550 nm in fog and haze for optical wireless communications," Proc. SPIE 4214, 26-37 (2001).
  13. D. T. Ha, V. V. Mai, and H. Kim, "Comparison of phase-screen-generation methods for simulating the effects of atmospheric turbulence," Korean J. Opt. Photonics 30, 87-93 (2019).
  14. F. G. Gebhardt, "Twenty-five years of thermal blooming: an overview," Proc. SPIE 1221, 2-25 (1990).
    CrossRef
  15. L. C. Bradley and J. Herrmann, "Phase compensation for thermal blooming," Appl. Opt. 13, 331-334 (1974).
    Pubmed CrossRef
  16. P. A. Konyaev and V. P. Lukin, "Thermal distortions of focused laser beams in the atmosphere," Appl. Opt. 24, 415-421 (1985).
    Pubmed CrossRef
  17. L. C. Andrews, R. L. Philips, R. J. Sasiela, and R. Parenti, "Beam wander effects on the scintillation index of a focused beam," Proc. SPIE 5793, 28-37 (2005).
    CrossRef
  18. J. M. Martin and S. M. Flatté, "Simulation of point-source scintillation through three-dimensional random media," J. Opt. Soc. Am. A 7, 838-847 (1990).
    CrossRef
  19. V. V. Mai, D. T. Ha, and H. Kim, "Link availability of terrestrial free-space optical communication systems in Korea," Korean J. Opt. Photonics 29, 77-84 (2018).
  20. E. Ostertagová, "Modelling using polynomial regression," Procedia Engineering 48, 500-506 (2012).
    CrossRef
  21. H. Kotake, Y. Abe, D. R. Kolev, Y. Saito, Y. Takahashi, T. Fuse, Y. Satoh, T. Itahashi, S. Yamakawa, H. Tsuji, and M. Toyoshima, "Experimental analysis of atmospheric channel model with misalignment fading for GEO satellite-to-ground optical link using 'LUCAS' onboard optical data relay satellite," Opt. Express 31, 21351-21366 (2023).
    Pubmed CrossRef
  22. T. J. Karr, J. R. Morris, D. H. Chambers, J. A. Viecelli, and P. G. Cramer, "Perturbation growth by thermal blooming in turbulence," J. Opt. Soc. Am. B 7, 1103-1124 (1990).
    CrossRef

저널정보

Optical Society of Korea

April 2024
Vol.35 No.2

pISSN 1225-6285
eISSN 2287-321X

Title: Korean Journal of Optics and Photonics
Abbreviation: Korean J. Opt. Photon.

Share this article on :

  • line