2.1 열과 유동 모델

레이저 조사로 인한 비금속 재료의 손상을 예측하기 위해 상용 프로그램(COMSOL Multiphysics 6.2)을 사용하여 열-유체 연계 해석을 수행하였다. 비정상 축 대칭 2차원 지배 방정식인 에너지(식 (1)), 질량(식 (2)), 운동량(식 (3)) 보존 방정식은 다음과 같다.

식 (1)에서 S_laser는 레이저 열원이며, Q_vap는 기화 잠열에 의한 손실 열원이다. 발생된 유체의 특성을 고려하여 유체 영역은 비압축성 층류 뉴턴 유체로 가정하였다. 질량 보존 방정식은 기화로 인한 반동 압력을 구현하기 위해 식 (2)와 같이 소스 항을 추가하였다. 질량 보존 방정식의 소스 항에서 m˙은 기화 온도보다 높은 경우에만 값이 존재하며, δ(φ)는 재료와 기체의 경계면에서만 값을 가진다. 경계면에서 기화가 발생할 때, 질량 보존 방정식은 소스 항이 부여되고 반동 압력이 생성된다. 운동량 보존 방정식에는 손상 과정에서 발생되는 물리적 현상을 모사하도록 중력, 부력, Darcy damping force, 표면장력을 고려하였다. ρg→는 중력 항이며, 식 (5)를 통해 재료와 기체를 구분하여 액체 영역에 적용하였다. 부력은 Boussinesq 근사를 사용하였으며, 이는 액체 영역의 밀도 변화가 온도 변화에만 의존한다고 가정하며 액체 영역에 적용된다. Ku→ 는 Darcy damping force로, 용융 및 응고 현상을 구현하기 위해 사용되며, 이를 적용하여 고체 영역의 속도를 0에 가깝게 설정하였다. Shah and Volkov의 연구에 따르면, 캐비티 깊이와 스폿 중심에서의 표면 온도에 미치는 마랑고니 응력은 거의 없으며, 용융 배출은 반동 효과에 의해 지배된다고 알려져 있다^[7-9]. 따라서, 액체 표면에서 마랑고니 응력은 무시하였으며, 표면장력 계수는 일정한 상수로 가정하였다. 또한, 반동 효과와 관련된 증기 속도는 Courtois et al.이 제안한 식 (2)에 따라 계산되지만, 과소평가된 것으로 보고되었다^[5]. 이로 인해 반동 압력도 과소 예측되며, 결과적으로 침투 깊이의 감소를 초래한다. 실제 현상을 정확히 모사하기 위해서는 액체와 기체 간의 경계면에서 올바른 힘 균형을 유지하는 것이 필요하며, 적절한 표면장력 계수를 적용하는 것이 요구된다. 이를 위해 Courtois et al.은 다양한 출력의 정적 레이저를 이용하여 적절한 표면장력 계수를 도출하였다^[5]. 비금속 재료에 대한 표면장력 계수는 Courtois et al.^[5]에 의해 제안된 접근법을 사용하여 결정하였다. Level Set Method는 고정 격자를 사용하며, 기체에서 레벨셋 변수가 0, 재료에서는 1의 값을 갖는다. 재료와 기체의 경계면에서는 레벨셋 변수가 0.5로 간주된다. 식 (4)는 경계면에서 레이저 조사로 인한 기화 현상을 정확하게 해석하기 위해 수정된 레벨셋 기법을 적용하였다. γ_ls는 전체 계산 도메인에서 예상되는 유체의 최대 속도 크기로 정의된다. 만약 γ_ls가 너무 크면 경계면이 부정확하게 해석되는 수치적 문제가 발생할 수 있으며, γ_ls가 너무 작으면 경계면 두께가 유지되지 않고 경계면이 진동할 수 있다. 따라서, 기화가 발생할 경우 γ_ls는 기체의 최고 속도인 50 m/s로 설정하였으며, 기화가 발생하지 않는 경우에는 용융점 기준으로 고체에서는 0.001 m/s, 액체에서는 1 m/s로 적절하게 설정하였다. ε_ls는 일반적으로 격자 크기의 절반으로 설정된다. 만약 ε_ls가 너무 크면 재료와 기체의 경계면을 정확하게 추적하지 못하고, ε_ls가 너무 작으면 조밀한 격자를 사용해야 하므로 계산 비용이 증가한다. 따라서, 격자 의존성 검사를 수행하여 수렴된 격자 크기를 선정한 후, 이 크기의 절반인 0.015 mm로 설정하였다. 식 (5)와 식 (6)은 온도 조건과 레벨셋 변수에 따라 전체 계산 도메인에서 각 상에 맞는 밀도와 점성계수를 정의하며, ρ_material 및 μ_material은 용융점을 기준으로 온도에 따라 고체 물성에서 액체의 물성으로 변화한다. 포화 증기압과 질량 유량은 식 (7)과 식 (8)로 표현된다. β_r은 재확산 계수로, 계산의 편의를 위해 1차 근사를 적용하여 0으로 가정하였다. 에너지 보존 방정식에 식 (9)를 추가하여 기화 잠열에 의한 에너지 손실을 고려하였다.

2.2 용융과 응고 모델

전체 계산 도메인에서는 질량, 운동량, 에너지 보존 방정식을 사용하여 온도, 속도, 압력 값을 계산한다. 고체는 속도를 가질 수 있으므로, 힘이 작용하면 쉽게 움직여 실험 결과와 차이가 발생할 수 있고, 응고 현상을 구현할 수 있는 모델이 필요하다. 따라서, 고체의 속도를 0에 가깝게 만드는 두 가지 방법인 Darcy damping force^[10]와 Enhanced viscosity^[11]를 적용하였다.

Darcy damping force의 계수는 식 (10)으로 표현된다. f₁의 값에 따라, 액체에서는 K가 0, 고체에서는 K가 -C/b가 된다. 이 비율은 Courtois et al.^[5]에서 제시된 값을 참고하여 모든 비금속 재료에 대해 동일한 값을 적용하였다. 식 (11)은 용융 분율을 나타내며, 고상선 온도보다 낮을 때는 0, 액상선 온도보다 높을 때는 1이다. 또한, 고상선 온도와 액상선 온도 사이의 구간에서는 용융 분율이 선형적으로 변화하도록 적용하였다. 식 (12)는 온도 변화에 따른 점성계수를 나타내며, 고체의 점성계수는 고체의 속도를 억제하기 위해 고점도의 유체로 가정하였다. H(T_melting)는 연속 2차 도함수를 가진 헤비사이드 함수이며, 고상선 온도와 액상선 온도 사이의 구간에서 점성계수가 연속적으로 변화하도록 적용하였다. 식 (13)은 등가 비열법을 적용하여 고상선 온도와 액상선 온도 사이의 구간에서 재료의 용융 잠열을 고려하였다.

2.3 레이저 열원 모델

레이저 열원은 재료의 표면에 조사된다. 레이저 빔의 공간 분포는 가우시안 분포를 적용하였으며, 레이저 열원은 식 (14)와 같다.

cos(θ)는 레이저 빔의 입사 각도에 따른 영향을 고려하기 위해 사용하였다. δ(φ)는 재료와 기체의 경계면에 에너지를 적용하기 위한 델타 함수(Delta function)이다. A는 재료별 흡수 계수를 나타낸다. 강철의 흡수 계수는 약 30 %이며, 키홀이 깊게 형성되면 다중 반사로 인해 흡수 계수가 증가한다. Courtois et al.은 기존의 연구를 바탕으로 다중 반사를 고려한 전체 흡수 계수를 70 %로 설정하였으며, DP600 스틸의 흡수 계수(약 30 %)에 계수 2.5를 곱하여 보정하였다^[5]. 게르마늄의 흡수 계수는 Desnijder and Vandewal에 의해 제공된 값을 사용하였다^[12]. 이는 식 (15)로 표현되며, 1.07 μm 파장에서 게르마늄의 광 침투 깊이는 매우 작기 때문에, 입사된 광은 표면에서 모두 흡수된다고 가정하였다. 또한, Desnijder and Vandewal의 연구에서는 레이저 조사로 인해 키홀이 형성되지 않으므로 흡수 계수를 보정하지 않았다. 실리콘의 흡수 계수는 Wang et al.의 연구를 참고하여 0.67로 설정하였다^[13]. 1.064 μm 파장에서 실리콘의 광 침투 깊이는 매우 작기 때문에, 입사된 광은 표면에서 모두 흡수된다고 가정하였다. Wang et al.의 실험에서는 레이저 출력 밀도가 636,500 W/cm²일 때, 2 ms 동안의 레이저 조사로 약 1.4 mm 깊이의 크레이터가 형성되었다^[13]. Courtois et al.^[5]에 의해 제안된 다중 반사를 고려한 흡수 계수 보정은 흡수 계수가 1을 초과하게 되어 물리적으로 적합하지 않다. 현재까지 실리콘에 대해 다중 반사를 고려한 전체 흡수 계수가 보고되지 않았으므로, 흡수 계수는 다중 반사를 고려하지 않은 값인 0.67^[13]과 최대 흡수 계수인 1로 설정하여 해석을 수행하였다.

2.4 계산 도메인 및 경계 조건

Fig. 1은 재료별 계산 도메인을 나타낸다. 게르마늄은 Desnijder and Vandewal의 계산 도메인을 참고하여 반경 40 mm, 두께 4 mm인 실린더로 설정하였다^[12]. 실리콘은 Wang et al.의 연구에서 사용된 시편을 참고하여 반경 15 mm, 두께 4 mm인 실린더로 설정하였다^[13]. 본 해석에 적용된 시편의 반경은 레이저 빔 반경에 비해 충분히 크므로, 추가적인 에너지의 손실은 고려하지 않았다.

Fig. 1.

Computational domain of germanium and silicon

Fig. 2는 경계 조건을 나타낸다. 레이저 열원은 중심 축을 기준으로 조사된다. 전체 도메인의 초기 온도는 300 K, 초기 압력은 1 bar, 초기 속도는 0으로 설정된다. 레이저 빔의 크기에 비해 시편의 크기가 충분히 크기 때문에 시편의 외곽에는 단열 조건이 설정되며, 외곽에서 속도는 0으로 고정된다. 또한, 레이저 조사로 인해 생성된 고온의 증기가 도메인 밖으로 빠져나갈 때 도메인 내부의 온도에 큰 변화가 없을 것으로 예상되므로 외부 경계에는 단열 조건과 함께 Open flow 및 1 bar 조건이 설정된다. 초기 레벨셋 변수는 기체에서 0, 재료에서 1의 값을 가지며, 초기 경계면은 레벨셋 변수가 0.5로 간주된다.

Fig. 2.

Schematic of boundary conditions

2.5 격자 의존성 검사

격자 의존성 검사는 Sukumar and Kar^[14]의 연구를 참고하여 수행하였다. 격자 의존성 검사에서 사용된 계산 조건은 Wang et al.^[13]의 실험을 참고하여 설정되었다. 이 실험에서 레이저 출력 밀도는 636,500 W/cm², 레이저 빔 직경은 0.6 mm이다. 따라서, 레이저 출력은 1800 W로 계산된다. 레이저 조사 시간이 2 ms, 재료는 실리콘이다. 계산 도메인은 Fig. 1(b)와 같으며, Table 1은 계산에 사용된 격자와 경계면 두께를 요약하였다. 0.08 mm, 0.06 mm, 0.04 mm, 0.03 mm 크기의 균일한 자유 삼각형 격자가 사용되며, r 방향과 z 방향의 격자 요소 수는 각각 188×69, 250×92, 375×138, 500×184로 설정된다. 경계면 두께는 격자 크기의 절반으로 설정되어 각각 0.04 mm, 0.03 mm, 0.02 mm, 0.015 mm가 된다. Fig. 3은 계산에 사용된 격자 분포를 나타낸다. 실리콘의 세 가지 상의 열물리적 특성은 Table 2에 나열되어 있다. 고체와 액체의 열전도율, 밀도, 비열은 주로 Wang et al.^[16]이 제공한 값을 사용하였다. 고상선 온도와 액상선 온도 사이의 간격은 5 K 를 넘지 않은 것으로 알려져 있으므로 Tao et al.^[15]에 의해 제시된 용융점(1683 K)을 기준으로 4 K로 가정하였다. 따라서, 고상선 온도와 액상선 온도는 각각 1681 K, 1685 K로 계산된다. 실리콘의 표면장력 계수는 Wang et al.^[13]에 의해 제공된 0.85 N/m 값을 기준으로 Courtois et al.^[5]에 의해 제안된 접근법을 사용하여 0.085 N/m로 계산된다.

Table 1.

Mesh and Interface thickness

Fig. 3.

Mesh configuration about case

Table 2.

Thermophysical properties of the three phases of silicon

손상의 지표인 크레이터 깊이(Crater depth)를 기준으로 격자 요소 수에 따른 결과를 분석하였다. Fig. 4는 다양한 격자 요소 수에 대해 레이저 조사 시간에 따른 크레이터 깊이를 나타낸다. 크레이터 깊이는 초기 경계면에서 액체와 기체의 경계면까지의 거리를 나타낸다. 0.04 mm(Grid 375×138)이거나 이보다 작은 격자 크기에서는 결과에 변화가 없다. 따라서, 이후 해석에서는 0.03 mm(Grid 500×184) 격자 크기가 사용된다.

Fig. 4.

Crater depth versus laser irradiation time for different grid element numbers

3.1 게르마늄

격자 의존성 검사를 통해 선정된 격자 크기를 사용하여 해석을 수행하였다. Desnijder and Vandewal^[12]이 수행한 시뮬레이션 조건에서 레이저 출력은 1000 W, 레이저 조사 시간은 3 s, 레이저 빔 직경은 2.4 cm이다. 계산 도메인은 Fig. 1(a)와 같으며, 열물리적 특성은 Table 3에 나열된 값과 같다. 고체의 열전도율, 밀도, 비열은 Desnijder and Vandewal^[12]이 제공한 값과 동일하다. 고상선 온도와 액상선 온도 사이의 간격은 5 K를 넘지 않은 것으로 알려져 있으므로 Desnijder and Vandewal^[12]에 제시된 용융점(1213 K)을 기준으로 4 K로 가정하였다. 따라서, 고상선 온도와 액상선 온도는 각각 1211 K, 1215 K로 계산된다. 게르마늄의 표면장력 계수는 Aqra and Ayyad^[18]에 제시된 0.583 N/m 값을 기준으로 Courtois et al.^[5]에 의해 제안된 접근법을 사용하여 0.0583 N/m로 계산된다.

Table 3.

Thermophysical properties of the three phases of germanium

Fig. 5는 레이저 조사 시간에 따른 온도 분포를 나타낸다. 온도 분포는 가우시안 분포와 유사하게 중심부에서는 온도가 더 높고, 주변으로 갈수록 감소하는 경향을 보였다. 2 s에서 최대 표면 온도가 1000 K에 도달하였으며, 3 s가 경과한 후 레이저 조사가 중단되면서 시편은 냉각되었다. 8 s에는 표면 온도가 400 K 에 이르렀다. Fig. 6은 게르마늄에 대한 레이저 조사 시간에 따른 시편 중심부의 표면 온도를 나타낸다. 본 수치 모델로 해석한 온도는 Desnijder and Vandewal^[12]의 시뮬레이션 결과보다 대체로 높게 나타났으며, 레이저를 2.5 s 동안 조사했을 때 최대 450 K의 차이를 보였다. Fig. 7은 다양한 레이저 조사로 인해 발생한 게르마늄의 손상을 보여준다. 빨간 원으로 표시된 손상 영역은 약 300 W/cm의 선형 출력 밀도와 2 s의 조사 시간으로 발생한 것으로 보고되었다^[21]. 선형 출력 밀도는 피크 레이저 출력을 빔 직경으로 나눈 값이며, 가우시안 프로파일의 경우 선형 출력 밀도에 2배를 곱하여 계산하였다^[21]. 이를 역으로 계산하면 피크 레이저 출력은 약 360 W에 해당하며, 이는 현재 시뮬레이션한 레이저 출력보다 낮은 수준이다. 해석에서는 레이저가 재료 표면에 집중되어 급격한 온도 상승이 발생하지만, 열응력 현상은 고려되지 않으므로 균열이 발생하지 않는다. 그러나 실험에서는 레이저가 표면에 집중되어 급격한 온도 상승이 일어나고, 이로 인해 표면과 내부 사이에 온도 차이가 발생하여 열팽창과 응력이 발생하게 된다. 이 응력이 일정 한계를 초과하면 균열이 발생할 수 있다. 또한, 재료 표면의 중심부에서 용융점 이하의 온도에도 불구하고, 균열 부위의 열전도율과 흡수 계수의 변화로 인해 열이 집중되어 해당 부위에서 국소적인 용융이 일어나 구멍이 형성될 수 있다. 전반적으로, 본 수치 모델은 Desnijder 와 Vandewal^[12]의 시뮬레이션 결과에 비해 더 높은 표면 온도를 예측하며, 이는 열응력 증가로 이어질 수 있다. 이러한 결과는 Fig. 7에서 관찰된 게르마늄 렌즈의 구멍 및 균열을 본 수치 모델이 더 정확하게 모사할 수 있음을 보여준다.

Fig. 5.

Temperature distribution for germanium

Fig. 6.

Comparison with simulation results

Fig. 7.

Photograph of the germanium used in the experiment^[21]

3.2 실리콘

격자 의존성 검사를 통해 선정된 격자 크기를 사용하여 해석을 수행하였으며, Wang et al.^[13]의 실험 결과와 비교하였다. Fig. 8은 레이저 조사 시간에 따른 실리콘의 상 분포를 나타낸다. 상 분포는 액체의 부피 분율을 나타내며, 물질이 고체 상태일 때는 0, 액체 상태일 때는 1이다. Fig. 8에서 보이는 바와 같이, 액체 영역은 반동 압력에 의해 수직 및 수평 방향으로 변형되었다. Fig. 9는 각각 왼쪽과 오른쪽에서 압력 및 온도 분포를 나타낸다. 압력 분포에 따르면, 재료의 기화로 인해 크레이터 중심부에서 압력이 약 1.05 bar 로 증가하였다. 온도 분포에서는 크레이터 중심부의 최고 온도가 기화 온도인 3514 K에 도달하는 것을 확인할 수 있다. Fig. 10은 레이저 조사 시간에 따른 실리콘의 크레이터 깊이를 나타낸다. 해석 결과에 따르면, 크레이터의 깊이는 시간이 지남에 따라 로그 함수적으로 증가하였다. 이는 발생한 반동 압력에 의해 액체가 분출되고, 외부로 분출되지 못한 액체가 중력에 의해 다시 내부로 되돌아오면서 크레이터 깊이가 점차 감소하는 현상 때문으로 추정된다. 해석 결과, 실험 결과와 비교하여 크레이터 깊이의 오차율은 흡수 계수가 0.67일 때 약 34 %, 1일 때 약 9 %로 확인되었다. 이러한 오차는 두 가지 주요 원인으로 분석된다. 첫째, 균열로 인해 실리콘 표면의 흡수 계수가 증가할 수 있다. 균열이 발생한 부위에서는 레이저 빔이 가두어지면서 흡수 계수가 증가한다. 또한, 레이저 출력 밀도가 충분히 커서 깊은 크레이터를 생성할 경우, 흡수 계수가 추가로 증가할 수 있다. 둘째, 열응력으로 인한 구조적 변화를 고려하지 않았기 때문에 발생한 것으로 분석된다. 레이저 조사에 의한 급격한 온도 변화로 열응력이 발생하고, 균열이 생기면 일부 재료가 외부로 튕겨 나가며, 레이저가 내부로 더 쉽게 침투할 수 있다. 이로 인해 크레이터의 깊이가 더욱 깊어질 수 있다. 현재까지 비금속 재료를 대상으로 다중 반사에 의한 흡수 계수를 실험적으로 측정한 사례는 없지만, 향후 관련 실험이 수행된다면 본 수치 모델에 적절한 흡수 계수 값을 적용할 수 있을 것으로 판단된다.

Fig. 8.

Phase distribution for silicon

Fig. 9.

Pressure and temperature distribution for silicon

Fig. 10.

Comparison of crater depth for silicon