Dispersion Volume 기반 탄성 메타물질의 3차원 등주파수면 추출과 유효 탄성 텐서 역추정

2.1 Bloch–Floquet 유한요소 고유치 문제

주기 구조를 이루는 대표 단위셀 Ω와 격자벡터 a₁, a₂, a₃를 정의하고, 선형 탄성체에 대한 주파수 영역 운동방정식을 기반으로 Bloch–Floquet 경계조건을 적용하였다. 주기 경계에 놓인 대응 점들에 대해 변위장은 복소 위상 인자 e^ik·a_m를 만족하도록 연결되며, 이로부터 Bloch 파수 k에 종속되는 고유치 문제가 구성된다.

여기서, [K(k)]은 Bloch–Floquet 위상 조건이 반영되어 k에 의존하는 강성행렬이고, [M]은 단위셀의 밀도 분포와 형상으로부터 정의되는 질량행렬로서 k에 의존하지 않는다. ω와 φ는 각각 각진동수와 모드벡터를 나타내며, m = 1, 2, 3은 식 (2)에 사용된 격자벡터 a_m의 번호를 의미한다. 이 연구에서는 각 Bloch 파수 벡터 k에 대해 저차 모드부터 관심 주파수 대역을 포함하는 충분한 수의 고유진동수를 계산하였다.

2.2 3차원 k-space 샘플링 및 dispersion volume 구성

기존의 band diagram은 일반적으로 irreducible Brillouin zone 내의 고대칭 경로를 따라 고유진동수를 도시한다. 반면 이 연구에서는 제1 Brillouin zone 또는 그 대칭축약 영역 내에 3차원 구조화된 k-space 격자를 정의하고, 각 샘플링 점에서 Bloch-Floquet 유한요소 고유치 해석을 수행하였다. 3차원 파수 격자의 샘플링 점은 식 (4)와 같이 정의된다.

각 샘플링 점 k_ijl에서 계산된 n번째 branch의 고유진동수를 branch-resolved scalar field로 저장하였으며, 이를 식 (5)와 같이 dispersion volume D_n으로 정의한다.

따라서, D_n(k_x, k_y, k_z)는 n번째 branch에 대한 3차원 고유진동수 분포를 나타내는 branch-resolved scalar field이다. 기존의 고대칭 경로 band diagram은 D_n을 특정 1차원 경로를 따라 절단하거나 보간한 단면으로 해석할 수 있으며, 등주파수면은 D_n(k) = ω₀를 만족하는 2차원 등가면으로 정의된다. 즉, dispersion volume은 단순한 3차원 데이터 저장 배열이 아니라, band diagram 재구성, IFS 추출, 방향성 전파 특성 분석, 그리고 Christoffel 방정식 기반 유효 탄성 텐서 역추정에 공통으로 사용되는 reciprocal-space 데이터 구조이다. k-space sampling 간격은 각 방향의 최대 및 최소 파수 범위와 샘플링 개수로부터 정의하였다. a = {x, y, x}에 대해 각 방향의 격자 간격은 식 (6)과 같다.

k-space sampling 간격이 IFS 재구성 결과에 미치는 영향을 확인하기 위해 서로 다른 격자 밀도 N_x = N_y = N_z = N_k에 대해 grid-sensitivity test를 수행하였다. 각 격자에서 동일한 목표 주파수와 branch에 대해 IFS를 재구성하고, 가장 조밀한 fine grid 결과를 기준으로 방향별 파수 반경의 상대 차이를 계산하였다. m번째 branch와 목표 주파수 f₀에 대한 sampling sensitivity error는 식 (7)과 같이 정의하였다.

여기서, $$ k_{Nk}^{\left(m\right)}$$는 $$ N_k^3$$ 격자에서 재구성한 IFS의 방향별 파수 반경이고, $$ k_{\mathrm{r}\mathrm{e}\mathrm{f} }^{\left(m\right)}$$는 가장 조밀한 기준 격자에서 얻은 값이다. Table 1은 f₀ = 0.800 kHz와 f₀ = 13.000 kHz에서 수행한 grid-sensitivity test 결과를 나타낸다. Branch별 IFS error는 fine grid(11 × 11 × 11)에서 재구성한 IFS를 기준으로 방향별 파수 반경의 평균 상대 차이로 계산하였다.

Table 1

Grid-sensitivity test for dispersion-volume-based IFS reconstruction

2.3 Dispersion Volume 기반 3차원 등주파수면 추출

목표 주파수 ω₀가 주어지면, n번째 branch의 IFS는 2.2절에서 정의한 dispersion volume D_n에서 등주파수 조건을 만족하는 파수 벡터의 집합으로 정의된다. 이 연구에서 사용한 k-space 계산 영역을 Ω_k라 하면, n번째 branch의 IFS는 식 (8)과 같이 표현된다.

여기서, Fig. 1과 같이 Ω_k는 제1 Brillouin zone 또는 그 대칭축약 영역을 의미한다. 이 연구에서는 각 IFS 점에서 새로운 Bloch-Floquet 고유치 문제를 반복적으로 푸는 것이 아니라, 사전에 계산된 3차원 dispersion volume D_n(k_x, k_y, k_z)를 사용하였다. 구체적으로, 구조화된 k-space 격자에서 정의된 D_n에 대해 trilinear interpolation을 적용하고, D_n(k)- ω₀ = 0을 만족하는 zero-level surface를 marching-cubes 계열의 등가면 추출 알고리즘을 통해 재구성하였다. 따라서 이 연구에서 IFS는 직접 고유치 재계산 결과가 아니라, 전체 3차원 k-space에서 미리 계산된 FEM 분산 정보를 이용한 interpolation-based iso surface reconstruction 결과이다.

Fig. 1

Unit-cell geometry, Bloch–Floquet boundary condition, and reciprocal-space setting

동일한 목표 주파수에서 서로 분리된 복수의 등주파수면이 존재할 수 있으므로, 추출된 등가면의 연결 성분은 branch별로 독립적으로 식별하였다. 또한 방향성 전파 특성을 해석하기 위해 Fig. 2와 같이 dispersion volume의 기울기로부터 군속도 방향을 계산하였다. 즉, D_n이 n번째 branch의 고유진동수장인 경우 군속도는 식 (9)와 같이 근사된다.

Fig. 2

Construction of high-symmetry-path dispersion curves from the sampled three-dimensional k-space data

Fig. 3은 첫 세 acoustic branches에 대해 구성된 dispersion volume D_n(k_x, k_y, k_z)의 직교 단면과 대표 등가면을 보여준다. 이 그림은 구조화된 Bloch-wave sampling grid에서 계산된 branch-resolved frequency field가 이후 Fig. 4의 IFS 추출과 Fig. 5의 Christoffel 기반 비교에 사용되는 기초 데이터임을 나타낸다.

Fig. 3

Interpolated dispersion volumes of the first three acoustic branches in three-dimensional k-space. Each panel shows orthogonal slices of the frequency field stored on the structured Bloch-wave sampling grid, together with a semi-transparent representative shell that illustrates the volumetric basis used for subsequent isofrequency-surface extraction

Fig. 4

Frequency-dependent evolution of dispersion-volume-based isofrequency surfaces. Representative IFSs reconstructed from the dispersion volume are shown for selected branches at increasing target frequencies. The surfaces remain smooth and nearly ellipsoidal in the low-frequency regime, whereas higher frequencies induce direction-dependent curvature variation and geometric distortion

Fig. 5

Low-k validation and finite-frequency comparison of FEM-reconstructed and Christoffel-based IFSs. The top row compares near-Γ FEM slope points (black dots) with Christoffel predictions (translucent surfaces) at 0.800 kHz, while the middle row compares first-octant IFSs at 13.000 kHz. The bottom row shows the frequency-dependent mean radial mismatch, principal-direction deviation, and normalized shape mismatch. Columns and curve colors denote acoustic branches 1~3.

2.4 Christoffel 방정식 기반 유효 탄성 텐서 역추정

장파장 영역에서는 주기 구조를 방향의존적 유효 탄성체로 근사할 수 있다. 이 연구에서는 2.2절에서 정의한 dispersion volume D_n 중 Γ점 근방의 세 acoustic branches에서 추출한 다방향 위상속도 데이터를 이용하여 유효 탄성 텐서 C^eff를 역추정하였다. 피팅에 사용한 q번째 파수 벡터 k_q에 대해 전파방향 단위벡터와 FEM 기반 위상속도는 식 (10)과 같이 정의된다.

전파방향 nq에 대한 Christoffel 행렬은 유효 탄성 텐서 C^eff를 이용하여 식 (11)과 같이 구성하였다.

여기서, λ_m은 Christoffel 행렬 Γ_ij의 m번째 고유값이고, ρ_eff는 단위셀의 체적 평균 밀도이다. 유효 탄성 텐서의 독립 성분은 parameter vector θ로 매개화하였으며, 선택한 material symmetry와 positive definiteness 조건을 만족하도록 제한하였다. 역추정은 FEM 기반 위상속도와 Christoffel 기반 위상속도 사이의 상대오차를 최소화하는 문제로 식 (12)와 같이 정식화하였다.

여기서, Q는 low-k fitting에 사용한 파수 벡터 집합, W_qm은 branch별 가중치, P(θ)는 positive definiteness 및 선택한 material symmetry constraint 위반에 대한 penalty 항이다. Branch matching은 Γ점 근방 acoustic branches의 속도 순서와 모드 연속성을 기준으로 수행하였다. 이 연구에서는 MATLAB의 비선형 최소제곱 최적화(lsqnonlin)를 사용하여 목적함수를 최소화하였다. 초기값은 별도의 초기 텐서가 지정되지 않은 경우 등방성 초기 탄성 텐서(E₀ = 1GPa, v₀ = 0.3)로 설정하였으며, 실제 최적화는 해당 초기 텐서의 Cholesky 하삼각 인자를 매개변수화한 벡터를 이용하여 수행하였다. 역추정된 C^eff로부터 계산한 Christoffel 기반 IFS는 이후 FEM 기반 IFS와 비교하여 local/quasi-static effective tensor model의 적용 가능성을 평가하였다.

2.5 FEM 기반 IFS와 Christoffel 기반 IFS의 비교 지표

유효물질 모델과 FEM 결과의 차이를 정량적으로 확인하기 위하여 동일 주파수에서의 평균 반경 오차, 주축 방향의 일치성, 그리고 IFS 형상 왜곡 정도를 비교하였다. 특히 저주파 영역에서는 두 IFS의 전체 형상과 장축/단축 비가 유사한지 확인하였고, 주파수가 증가함에 따라 오차가 어떻게 누적되는지 분석하였다.

동일 주파수 f에서 FEM기반 IFS와 Christoffel 기반 IFS의 차이는 세 가지 지표로 평가하였다. 첫째, 평균 반경 mismatch는 방향 n에 따른 파수 반경 차이의 평균 상대오차로 식 (13)과 같이 정의하였다.

둘째, 주축 방향 편차는 FEM 기반 IFS와 Christoffel 기반 IFS의 주성분 방향 벡터 e_Chr, e_FEM 사이의 각도로 식 (14)와 같이 계산하였다.

셋째, normalized shape mismatch는 평균 반경으로 정규화한 IFS 반경 분포의 차이로 식 (15)와 같이 정의하였다.

여기서, $$ \tilde{k}$$는 각 IFS의 평균 반경으로 정규화한 방향별 파수 반경이다.

3.1 기존 대칭경로 분산곡선과 3차원 Dispersion Volume의 비교

고대칭 경로를 따라 계산된 분산곡선은 밴드 구조의 전반적 추세와 밴드갭 위치를 파악하는 데 유용하다. 그러나 이러한 분산곡선은 선택된 1차원 경로상의 고유진동수만을 제공하므로, 동일 주파수에서 허용되는 파수 벡터가 전체 3차원 k-space에서 어떤 곡면을 이루는지는 충분히 보여주지 못한다. 따라서 대칭경로 band diagram만으로는 IFS의 전체 폐곡면 또는 개곡면 구조, 방향별 곡률 변화, 국소 평탄화 및 비구면성을 판단하기 어렵다.

이 연구에서 구성한 dispersion volume은 이러한 한계를 보완한다. 동일한 branch-resolved frequency field D_n(k_x, k_y, k_z)로부터 목표 주파수에서의 등주파수면을 재구성함으로써, 특정 branch의 IFS가 구면에 가까운지, 타원체형인지, 또는 주파수 증가에 따라 방향별 곡률 변화와 국소적 왜곡을 보이는지를 확인할 수 있다. Fig. 4는 dispersion volume에서 추출한 대표 IFS의 주파수별 변화를 보여주며, 저주파 영역에서는 비교적 매끄럽고 타원체에 가까운 형상이 나타나는 반면, 주파수가 증가함에 따라 방향별 곡률 변화와 표면 왜곡이 증가함을 보여준다. 따라서 제안된 방법은 기존 band diagram을 대체하기보다, 대칭경로 분산곡선으로는 확인하기 어려운 3차원 IFS 형상과 방향성 전파 특성을 보완적으로 제공하는 framework로 이해할 수 있다.

3.2 주파수별 3차원 등주파수면 변화

Fig. 4는 dispersion volume에서 재구성한 대표 IFS의 주파수별 변화를 보여준다. 저주파 영역에서는 IFS가 Γ점 주변에서 비교적 매끄럽고 폐곡면 형태를 유지하며, 그 형상은 대체로 완만한 타원체에 가깝다. 이는 장파장 근사 하에서 등가 연속체적 거동이 지배적임을 의미한다. 반면 주파수가 증가함에 따라 IFS의 곡률은 방향별로 다르게 변화하고, 일부 영역에서는 표면의 국소적 평탄화 또는 왜곡이 나타난다.

이러한 주파수별 IFS 변화는 대칭경로 band diagram만으로는 확인하기 어려운 3차원 방향성 정보를 제공한다. 특히 IFS의 비구면성, 방향별 곡률 변화, 국소 평탄화는 군속도 방향의 편향 및 에너지 유동의 비등방성과 관련될 수 있다. 따라서 Fig. 4의 결과는 제안한 dispersion-volume 기반 IFS 재구성이 탄성 메타물질의 방향성 파동 전파 특성을 평가하는 데 유용함을 보여준다.

3.3 저주파 유효 탄성 텐서의 검증

Fig. 5의 상단 행은 0.800 kHz에서 세 acoustic branches에 대한 low-k self-consistency validation 결과를 보여준다. 2.4절에서 역추정한 유효 탄성 텐서 C^eff를 이용하여 Christoffel 방정식 기반 IFS를 계산하고, 이를 dispersion volume에서 재구성한 FEM 기반 near-Γ IFS와 비교하였다. 세 branches 모두에서 FEM으로부터 추출한 near-Γ slope points와 Christoffel 기반 IFS가 주축 방향, 장축/단축 비, 그리고 전체 곡률 경향에서 잘 일치하였다.

이 결과는 역추정된 C^eff가 장파장 영역에서 주기 구조의 방향의존적 위상속도를 적절히 재현함을 보여준다. 다만 이는 유효 탄성 텐서가 모든 유한 주파수 영역에서 동일하게 적용된다는 의미가 아니라, Γ점 근방의 low-k acoustic branch data에 대한 local/quasi-static effective tensor model의 자기일관성 검증으로 해석된다. 따라서 이후 3.4절에서는 동일한 C^eff로부터 계산한 Christoffel 기반 IFS가 주파수 증가에 따라 FEM 기반 IFS를 어느 정도까지 재현할 수 있는지를 비교한다.

3.4 주파수 증가에 따른 local effective tensor model의 한계

Fig. 5의 중간 행은 13.000 kHz에서 dispersion volume으로부터 재구성한 FEM 기반 IFS와 2.4절에서 역추정한 C^eff로부터 계산한 Christoffel 기반 IFS를 비교한 결과를 보여준다. 저주파 Γ점 근방에서는 두 결과가 잘 일치하지만, 주파수가 증가함에 따라 평균 반경, 주축 방향, 그리고 표면 형상에서 branch별 차이가 나타나며, 전반적으로 local effective tensor model의 재현성이 저하되는 경향을 보인다. 이 연구에서는 2.5절에서 정의한 평균 반경 mismatch $$ {ϵ}_r^{\left(m\right)}\left(f\right)$$, 주축 방향 편차 Δθ^(m)(f), normalized shape mismatch $$ {ϵ}_s^{\left(m\right)}\left(f\right)$$를 이용하여 FEM 기반 IFS와 Christoffel 기반 IFS 사이의 차이를 정량화하였다.

Fig. 5의 하단 행은 평균 반경 mismatch, 주축 방향 편차, normalized shape mismatch의 주파수 의존성을 보여준다. 이러한 지표는 저주파에서 식별한 local/quasi-static C^eff가 유한 주파수 IFS를 어느 정도까지 재현할 수 있는지를 평가하는 기준으로 사용하였다.

따라서, FEM 기반 IFS와 Christoffel 기반 IFS 사이의 finite-frequency mismatch는 동적 유효물성의 직접적 증거로 단정하기보다, Γ점 근방에서 식별한 local/quasi-static effective tensor model의 적용 범위가 주파수 증가에 따라 제한됨을 보여주는 결과로 해석하는 것이 타당하다. 이러한 차이는 유효물성의 주파수 의존성을 배제하지 않지만, spatial dispersion, nonlocal effects, finite-wavelength effect 및 branch coupling의 영향과도 관련될 수 있다.

Table 2에서 볼 수 있듯이, 이 연구에서 채택한 $$ {ϵ}_r^{\left(m\right)}$$ < 5 %, Δθ^(m) < 3°, $$ {ϵ}_s^{\left(m\right)}$$ < 5 % 기준을 동시에 적용할 경우, acoustic branches 1, 2, 3의 유효 주파수 범위는 각각 0.800 kHz ~ 2.600 kHz, 0.800 kHz ~ 2.400 kHz, 0.800 kHz ~ 5.400 kHz로 나타났다. 세 branch 모두에서 limiting metric은 평균 반경 mismatch $$ {ϵ}_r^{\left(m\right)}$$로 확인되었으며, 이는 주파수 증가에 따라 IFS의 평균 반경 차이가 local effective tensor model의 적용 범위를 우선적으로 제한함을 의미한다.

Table 2

Branch-wise validity range of the local effective tensor model based on the adopted mismatch criteria

3.5 유효 탄성 텐서 역추정 결과와 대칭성 비교

2.4절의 Christoffel 방정식 기반 역추정 절차를 이용하여 Γ점 근방의 저주파 acoustic branches로부터 simulation frame에서의 유효 탄성 텐서 C^eff를 식별하였다. 이 절에서의 검증은 외부 정적 균질화 reference tensor 또는 문헌 reference value와의 비교가 아니라, 역추정된 C^eff가 fitting에 사용된 raw low-k FEM velocity data를 얼마나 일관되게 재현하는지를 평가하는 자기일관성 검증(self-consistency validation)이다. 따라서 아래의 오차는 역추정된 C^eff가 Γ점 근방의 FEM 기반 저주파 분산 응답을 얼마나 잘 재현하는지를 나타내는 지표로 해석된다.

이 연구에서는 적합 코드와 동일하게 Voigt 순서를 [11 22 33 23 13 12]로 사용하였다. 식별된 full 6 × 6 tensor는 raw low-k velocity data에 대해 median 상대오차 0.54 %와 RMS 상대오차 0.95 %를 보였으며, 이는 장파장 영역의 방향의존적 위상속도를 작은 오차로 재현함을 의미한다. 계산된 유효 탄성 텐서는 식 (16)과 같다.

Fig. 5의 상단 행은 0.800 kHz에서 세 acoustic branches에 대한 low-k self-consistency validation 결과를 보여주며, 역추정된 C^eff로부터 계산한 Christoffel 기반 IFS가 FEM으로부터 얻은 near-Γ slope points를 잘 재현함을 확인할 수 있다. 반면 Fig. 5의 하단 행은 동일한 C^eff를 유한 주파수 영역에 적용할 때 branch별 mismatch가 증가함을 보여주며, 이는 3.4절에서 논의한 local/quasi-static effective tensor model의 적용 한계와 연결된다.

Reference tensor가 없는 경우 material symmetry 선택은 외부 reference value와의 차이에 의존하지 않고, 단위셀의 geometric symmetry, 독립 탄성상수의 수, positive definiteness 만족 여부, 그리고 동일한 low-k FEM velocity data에 대한 재현 오차를 함께 고려하여 수행하였다. 이를 위해 full tensor, 최적 회전을 허용한 orthotropic projection, tetragonal projection 및 cubic reduction을 비교하였다. 그 결과 orthotropic projection은 median 상대오차 0.53 %, RMS 상대오차 0.96 %로 full tensor와 거의 동일한 정확도를 유지하였다. 반면 tetragonal projection은 median 상대오차 1.74 %, RMS 상대오차 2.10 %로 오차가 증가하였으며, cubic reduction은 단위셀의 기하학적 대칭성과 fitting error 경향을 고려할 때 이 구조의 유효 탄성거동을 설명하는 모델로 채택하지 않았다. 따라서 이 구조의 저주파 유효 탄성거동은 회전된 재료 좌표계에서 orthotropic class로 해석하는 것이 가장 타당하며, 본문에서는 simulation frame의 full tensor를 기본 보고값으로 제시하였다.