Original Article

Journal of Computational Fluids Engineering. 30 September 2026. 163-173
https://doi.org/10.6112/kscfe.2026.31.3.163

ABSTRACT


MAIN

  • 1. 서 론

  • 2. 수치기법

  •   2.1 Level Contour Reconstruction Method(LCRM)

  •   2.2 지배 방정식

  •   2.3 수치해석 구성

  • 3. 결 과

  •   3.1 수렴성 검증

  •   3.2 Reynolds 수 및 점도비에 따른 최대 퍼짐 거동

  •   3.3 점도비 𝛽 변화에 따른 최대 퍼짐 거동

  •   3.4 Weissenberg 수에 따른 최대 퍼짐 거동

  • 4. 결 론

1. 서 론

잉크젯 프린팅, 스프레이 코팅, 식품 공정 및 생체의학 응용 등 다양한 산업 분야에서 액적이 고체 표면에 충돌하는 현상은 공정 성능을 결정하는 핵심 요소로 작용한다. 액적 충돌 과정에서는 관성력, 점성력, 표면장력 및 고체 표면의 젖음성이 복합적으로 상호작용하여 퍼짐(spreading), 수축(recoiling), 진동(oscillation), 비산(splashing) 등의 다양한 동역학적 거동이 나타나며, 그 결과 액적의 형상 변화와 최종 부착 상태가 결정된다[1]. 특히, 충돌 직후의 관성 지배적 퍼짐 단계를 지나 도달하게 되는 최대 퍼짐 직경은 코팅 균일성, 액적 부착성, 인쇄 해상도 및 공정 효율을 결정하는 핵심 설계 인자로 알려져 있다[1]. 따라서 액적의 최대 퍼짐 거동을 정확히 예측하고 제어하는 것은 다수의 열유체 및 제조 공정에서 중요한 과제이다.

지금까지의 액적 충돌 연구는 주로 뉴턴 유체를 대상으로 수행되어 왔다. 뉴턴 유체의 경우 점도가 일정하고 전단응력과 변형률 속도가 선형 관계를 가지므로, 충돌 후 퍼짐 거동은 관성력, 점성력, 표면장력의 상대적 크기와 표면 젖음성을 바탕으로 비교적 잘 정량화될 수 있다. Clanet 등[2]은 저점도 액적이 초소수성 표면에 충돌할 때 최대 퍼짐 직경이 Dmax/D0~We1/4형태의 단순한 스케일링을 따른다는 사실을 실험적으로 보였으며, 이후 Reynolds 수, Weber 수, Ohnesorge 수와 같은 무차원수 및 접촉각 조건을 이용한 다양한 최대 퍼짐 예측 모델이 실험적·이론적·수치적으로 제안되어 왔으며[1,2], 이러한 연구들은 에너지 보존, 점성 소산을 기반으로 액적 충돌 현상을 이해하는 기준 틀을 제공하였다.

그러나 실제 산업, 식품 및 생물학적 응용에서는 많은 경우 비뉴턴 유체를 다루게 된다[1]. 고분자 용액, 생체유체 등의 유체는 내부 미세구조의 영향으로 비뉴턴 특성을 나타내며, 이 경우 점도와 응력 응답이 비선형적으로 변화한다. 특히 점탄성 유체는 변형에 저항하며 에너지를 소산하는 점성 특성과 변형 후 원래 상태로 회복하려는 탄성 특성을 동시에 갖기 때문에, 충돌 과정에서 미세구조에 기인한 탄성 응력이 추가로 발생하여 유체의 움직임을 파악하기에 어려움이 존재한다. 따라서 뉴턴 유체 기반의 기존 예측 모델로는 점탄성 유체의 시간 의존적 응력 효과를 충분히 반영하기 어렵다.

이러한 한계를 극복하기 위해 최근에는 뉴턴 유체를 대상으로 한 기존 연구를 넘어, 비뉴턴 유체, 특히 점탄성 액적의 충돌 거동을 이해하기 위한 실험적, 수치해석적 연구가 활발히 수행되고 있다. 점탄성 액적 충돌에 대한 초기 연구는 주로 실험을 통해 수행되었다[3,4]. Bergeron 등[3]은 미량의 고분자 첨가만으로 액적 부착 거동이 크게 달라질 수 있음을 보고하였고, Bartolo 등[4]은 소수성 표면에 충돌하는 비뉴턴 액적에서 수축 속도 감소와 반발 억제를 관찰하고, 접촉선 부근의 비뉴턴 법선응력을 고려한 모델로 수축 속도 감소를 정량적으로 설명하였다. 이러한 실험 연구들은 점탄성 응력이 액적 충돌 거동에 중요한 역할을 한다는 사실을 보여주었으며, 비뉴턴 액적 충돌 연구의 필요성을 제시하였다. 그러나 실험적 접근만으로는 액적 충돌 과정에서 점도, 점탄성 응력 등 개별 물성 인자가 충돌 거동에 미치는 영향을 독립적으로 분리하여 평가하는 데 다소 어려움이 존재하였다.

최근 수치해석 기법의 발전과 함께 점탄성 유체의 자유표면 및 다상 유동을 해석하기 위한 다양한 수치 모델이 개발되어 왔다. 초기 연구들은 주로 점탄성 구성방정식을 계면 추적 기법과 결합하여 복잡한 계면 유동을 안정적으로 해석할 수 있는 수치 방법론을 구축하고, 대표적인 문제를 통해 그 적용성을 검증하는 데 초점을 두었다. Pillapakkam과 Singh[5]은 Level-Set 기법과 Oldroyd-B 모델을 결합한 점탄성 2상 유동 해석 기법을 제시하였으며, Chinyoka 등[6]은 VOF 기법을 이용하여 전단 유동에서 점탄성 액적의 변형을 해석하였다. 이후에도 다양한 점탄성 구성 모델과 계면 추적 기법이 개발되었으며, jet buckling 및 기타 점탄성 자유표면 유동과 같은 다양한 2유체 문제를 대상으로 수치 모델의 정확성과 적용 가능성을 검증하는 연구가 수행되어 왔다[7,8,9,10]. 따라서 기존 연구의 상당수는 특정 유동 조건에서 점탄성 효과에 의해 나타나는 거동을 해석함으로써 개발된 수치 모델의 성능과 물리적 특성을 평가하는 방향으로 발전해 왔다.

이러한 수치 모델의 적용 문제 중 하나로 점탄성 액적의 고체 표면 충돌도 연구되어 왔다[7,9,10,11]. Xu 등[9]은 Smoothed Particle Hydrodynamics(SPH) 기법을 이용하여 3차원 점탄성 자유표면 유동 및 액적 충돌을 해석하였으며, Viezel 등[10]은 Oldroyd-B 유체에서 점도비의 변화가 액적 충돌 거동에 미치는 영향을 분석하였다. 또한 Figueiredo 등[11]은 Giesekus 및 eXtended Pom-Pom(XPP) 모델을 이용하여 주요 유변학적 매개변수 변화에 따른 액적 변형과 퍼짐 특성을 조사하였다. 이러한 연구들을 통해 점탄성 액적 충돌을 해석하기 위한 수치 모델의 적용 가능성이 확인되었으며, 점도비와 응력 완화 특성 등 점탄성 관련 매개변수가 충돌 이후의 액적 변형에 미치는 영향이 제시되었다.

그러나 기존의 점탄성 액적 충돌 연구 역시 수치 모델의 개발 및 적용을 중심으로 특정 점탄성 매개변수의 변화에 따른 액적 거동을 분석하는 데 주로 초점이 맞추어져 있다. 이에 따라 점탄성 효과에 의해 발생하는 변화를 뉴턴 유체의 액적 충돌에서 잘 알려진 거동과 직접 연결하고, 서로 다른 특성의 두 유체 사이의 차이점에 관해서 연속적으로 비교 분석한 연구는 상대적으로 제한적이다. 따라서 본 연구에서는 고체 표면에 수직으로 충돌하는 Oldroyd-B 액적을 대상으로, 점탄성 유체만의 독립적인 충돌 특성을 분석하기보다는 기존에 확립된 뉴턴 유체의 액적 충돌 특성과의 연결성을 중심으로 최대 퍼짐 거동을 분석하였다. 이를 통해 점탄성 효과가 추가됨에 따라 기존 뉴턴 액적의 충돌 및 퍼짐 특성이 어떠한 방식으로 변화하는지를 체계적으로 파악하고자 하였다.

2. 수치기법

2.1 Level Contour Reconstruction Method(LCRM)

액적 충돌 시 경계면의 정확한 추적을 위해 기존 연구에서는 주로 Front Capturing 기법이나 Front Tracking 기법이 사용되어 왔다[12]. Front Capturing 기법(예: Volume of Fluid 또는 Level Set)은 고정된 직교격자 상에서 적절한 스칼라 변수(예: Volume of Fluid(VOF)의 color function이나 level set의 거리함수)를 이용하여 충돌 경계면을 표현하기 때문에, 구현이 비교적 단순하며 계면 변형을 다루기 쉽다는 장점이 있다. 그러나 각 기법은 고유한 수치적 문제(예: VOF 기반 기법에서의 매끄러운 계면 표현, Level Set 기반 기법에서의 거리함수 재초기화)를 해결하기 위해 별도의 보조 기법이 필요하다. Front Tracking 기법은 계면을 표현하기 위해 추가적인 이동격자를 도입하고 이를 라그랑지안 방식으로 추적한다. 라그랑지안 계면 요소 간의 논리적 연결성을 유지해야 하는 구현 부담이 존재하지만, Front Tracking 기법은 계면 현상을 매우 정확하게 표현할 수 있는 우수한 능력을 보여 왔다.

본 연구에서는 두 기법(Level Set과 Front Tracking)의 장점을 결합한 하이브리드 기법인 Level Contour Reconstruction Method(LCRM)[12]을 이용하여 액적 계면의 이동을 추적하였다. LCRM은 정확하게 계면을 추적하며, 주어진 계면 정보로부터 거리함수를 효율적으로 계산할 수 있는 장점이 있다. 이를 통해 Front Tracking의 높은 계면 추적 정확도와 Level Set의 안정적인 계면 재구성 능력을 함께 활용할 수 있다. 특히 계면 변형이 크게 발생하는 액적 충돌 문제에서, 라그랑지안 이동격자의 과도한 변형을 방지하고 계면 추적을 가능하게 한다.

Fig. 1은 LCRM의 기본 개념을 설명한다. 그림 좌측과 같이 상 경계면은 기본적으로 라그랑지안 이동격자로 표현된다. 거리함수는 상 경계면 위에서 정확히 0이 되도록 정의되므로, 각 격자의 셀 위에서 거리함수 ϕ = 0이 되는 두 지점을 연결함으로써 라그랑지안 이동격자를 생성할 수 있다(그림 상단의 확대 그림 참고). 이러한 위치는 적절한 보간법[12]을 이용해 손쉽게 계산할 수 있다. 상 계면의 재구성은 계산 격자 사이의 각 셀 면에서 수행되기 때문에, 다시 그려진 선 요소들은 셀 면 위에 위치한 교점을 공유함으로써 암시적으로 연결된다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310310/images/jkscfe_2026_313_163_F1.jpg
Fig. 1.

Basic concept of level contour reconstruction method(LCRM)

2.2 지배 방정식

본 연구에서 사용된 비압축성 유체에 대한 질량 및 운동량 보존을 포함한 Navier-Stokes 방정식은 다음과 같다. 본 연구에서는 액체와 기체를 하나의 계산 영역에서 해석하는 single-field formulation을 사용하였다.

(1)
∇∙u=0
(2)
ρ∂u∂t+u∙∇u=-∇p+∇∙τ+ρg+Fs

여기서 u는 속도, p 는 압력, 𝜌는 밀도, g는 중력가속도, Fs는 계면에 가해지는 표면장력을 나타낸다. 표면장력에 대한 보다 더 자세한 내용은 기존 문헌[12]에서 추가로 확인할 수 있다. 𝜏는 압력항을 제외한 전체 응력을 나타낸다. 전체 응력 𝜏는 뉴턴 점성 응력과 점탄성 응력으로 분해되며, 다음과 같이 표현된다.

(3)
τ=2ηsD+τp

여기서 ηs는 점탄성 액체의 뉴턴 점성 응력에 대응하는 점도이고, D=12∇u+∇uT는 변형률 속도 텐서이며, τp는 점탄성 응력이다. τp는 유동에 따라 시간적으로 변화하기 때문에 별도의 유변학적 구성방정식이 존재하며 아래와 같이 정의된다.

(4)
∂τp∂t+(u·∇)τp-(∇u)T·τp+τp·∇u=1λMτp

여기서 𝜆는 점탄성 유체에 변형이 가해진 후 발생한 점탄성 응력 τp이 시간에 따라 감소하여 평형 상태로 완화되는 특성 시간(relaxation time)을 의미한다. 또한, Mτp는 점탄성 구성 모델에 따라 다르게 정의되고, 대표적인 점탄성 유체 모델 Oldroyd-B에 대한 식은 다음과 같다.

(5)
Mτp=2ηpD-τp

여기서 ηp는 점탄성 응력에 대응하는 점탄성 기여 점도이다. 본 연구에서는 점탄성 유체 중에서도 비교적 간단한 Oldroyd-B 모델을 사용하였다[13]. 다른 모델에 관한 더 자세한 내용은 기존 문헌을 참고할 수 있다[8].

본 연구에서와 같이 액체가 고체 벽면과 만날 경우, 상이 만나는 경계인 접촉선(contact line)의 움직임의 정확한 표현이 중요하다. 상 경계면이 고체와 접촉하는 경우 접촉선은 고체 표면에서의 점착 조건으로 인하여 고체에 고정된다. 고체 주위의 상 경계면은 유동으로 인해 계속 이동하게 되고, 이는 2유체 유동에서 접촉 지점에 무한한 전단 응력을 발생시키면서 불안정한 계면의 움직임을 만들어 낸다. 본 연구에서는 이러한 계면의 불안정성이 발생하는 것을 방지하기 위하여 다음과 같이 접촉점에서 계면을 미끄러지게 하는 Navier-slip 경계 모델[14]을 사용하였다.

(6)
UCL=λs∂u∂nwall 

여기서 ∂u/∂n은 벽면에서의 속도구배를 나타낸다. λs는 미끄럼 길이로 본 연구에서는 격자 크기를 사용하였다.

추가적으로 접촉 지점의 움직임을 정확하게 설명하기 위해 동적 접촉각이 모델링 되어야 한다. 동적 접촉각은 접촉각 히스테리시스를 고려하기 위해, 각 접촉선에 대해 확장된 선 요소를 이용해 모델링된다. LCRM 기법에서는 계면 요소를 고체 벽 내부로 확장하고, 이 확장 요소가 동적 접촉각을 갖도록 함으로써 지배 방정식에 표면장력을 적절히 부여한다. 확장 각도 𝜃는 접촉각 히스테리시스를 고려하여 다음과 같이 결정된다.

(7)
θ=θaθ=θrθr<θ<θa if θ>θa if θ<θr otherwise. 

여기서 𝜃는 접촉각을 의미하고 하첨자 a와 r은 각각 전진각 및 후진각을 의미한다. 계산 시 접촉각이 전진각보다 큰 경우는 전진각으로 후진각보다 작은 경우는 후진각으로 고정되고, 전진각과 후진각 사이의 범위에서는 접촉선의 위치가 유지되도록 모델링하였다. 본 연구에서는 전진 접촉각 θa=119°, 후진 접촉각 θr=74°로 설정하였다. 접촉각 변화에 따른 영향을 확인하기 위해 후진각/전진각을 각각 44°/89°, 74°/119°, 104°/149°로 변화시켜 계산한 결과, 접촉각 히스테리시스에 대한 최대 퍼짐 직경의 영향은 무시할 수 있을 만큼 작은 것으로 확인하였다. 따라서 본 연구 조건에서는 접촉각 변화가 최대 퍼짐 직경에 미치는 영향이 무시할 수 있는 수준이라 판단하고, 이후 해석에서 접촉각 조건을 유지한 상태로 점도 및 점탄성 물성 변화가 최대 퍼짐 거동에 미치는 영향을 분석하였다. 접촉각에 대한 보다 더 자세한 내용은 참고문헌[14]에서 확인할 수 있다.

지배방정식(2)는 Chorin의 Projection Method[15]를 활용하여 수치적으로 계산되었다. 1차 오일러 기법을 활용하여 시간 차분화를 계산하였고, 공간 차분화를 위해 정렬 엇갈림 격자를 활용하였다. 식(4)의 상류대류미분(Upper-convected derivative)을 포함한 보다 자세한 수치기법 내용은 기존 문헌[8]에서 추가로 확인할 수 있다.

2.3 수치해석 구성

본 연구에서의 기하학적 형상 및 경계 조건을 포함한 수치해석 구성은 Fig. 2에서 확인할 수 있다. 2차원 축대칭 좌표계를 활용했으며, 고체 표면에는 점착 조건을 적용했다. 우측 및 상단 경계는 Open condition을 적용하였으며, 계산 영역을 충분히 크게 설정하여 경계 조건이 액적 퍼짐 거동에 미치는 영향을 최소화하였다. 초기 직경 D0 를 갖는 구형 액적이 초기 속도 U로 고체 표면을 향해 낙하하도록 하였다. 액적 초기 위치는 고체 표면으로부터 H만큼 떨어진 위치에 배치하였으며, 이에 따라 액적 초기 중심은 고체 표면으로부터 H+D0/2 높이에 위치한다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310310/images/jkscfe_2026_313_163_F2.jpg
Fig. 2.

Computational domain and boundary conditions setup

본 연구에서 사용된 주요 물성, 해석 변수 및 무차원수를 Table 1과 2에 정리하였다. 다만 Table 1에 제시한 ηs, ηp, η0, 𝜆의 값은 3.1절의 수렴성 검증 및 선행 연구와의 비교에 사용한 기준값이며, 이후 해석에서는 주요 변수로서 각 조건에 따라 변화시켜 적용하였다.

Table 1.

Reference values of physical properties and computational parameters

Parameter Reference value Unit Description
ρl 1000 [kg/m3] Liquid density
ρg 1 [kg/m3] Gas density
ηs 0.4 [Pa•s] Solvent viscosity of the viscoelastic droplet
ηp 3.6 [Pa•s] Polymeric viscosity
η0 4.0 [Pa•s] Total viscosity of the viscoelastic droplet(η0=ηs+ηp)
μg 0.00001 [Pa•s] Gas viscosity
𝜎 0.001 [N/m] Surface tension coefficient
g 9.81 [m/s2] Gravitational acceleration
𝜆 0.02 [s] Relaxation time
U 1.0 [m/s] Initial impact velocity
D0 0.02 [m] Initial droplet diameter
Table 2.

Definitions and physical interpretations of dimensionless numbers

Dimensionless number Definition Physical interpretation
Re ρlUD0/η0 Ratio of inertial forces to viscous forces
Wi λU/D0 Ratio of the relaxation time to the characteristic flow time
We ρlU2D0/σ Ratio of inertial forces to surface tension forces
𝛽 ηs/η0 Ratio of solvent viscosity to total viscosity

3. 결 과

3.1 수렴성 검증

본 연구에서 사용한 수치해석 기법의 정확성과 격자 의존성을 확인하기 위해, 이전 연구[8]에서 제시된 점탄성 액적 충돌 결과와 비교하여 검증을 진행하였다. 이전 연구[8]를 참고하여 초기 직경 D0=0.02 m, 초기 속도를 U=1 m/s, 초기 낙하 높이를 H=0.03 m로 설정하였다. 이 값을 기준으로 관련 무차원수는 Re=5, Wi=1, We=20000, 𝛽=0.1로 정리된다. 참고문헌을 고려한 본 연구에서 사용한 초기 액적 직경이 미소 액적 해석 조건에 비해 비교적 큰 대표길이에 해당하므로 Weber 수가 크게 나타난다. 본 계산은 시간에 따른 액적의 무차원 폭 변화를 기준으로 나타내었으며, 시간에 따라 변화하는 액적의 폭(D)은 초기 액적 직경 D0 로 나누어 D/D0 로 무차원화하였다. 또한 시간은 t*=tU/D0로 무차원화하였다.

Fig. 3은 점탄성 유체에 대한 이전 수치해석 연구 결과(검정색 원형 마커)[8], 뉴턴 유체에 대한 이전 수치해석 연구 결과(검정색 직사각형 마커)[9], 그리고 점탄성 유체에 대한 각 격자별 수치해석 결과(검정색 점선, 빨간색 실선, 초록색 실선) 및 뉴턴 유체에 해당하는 수치해석 결과(파란색 실선)를 함께 비교하여 보였다. 충돌 직후 뚜렷한 수축 없이 완만한 퍼짐 거동을 보이는 뉴턴 유체와 달리 점탄성 유체의 경우 점탄성 응력의 영향으로 수축 및 재팽창 거동을 하게 된다. 뉴턴 유체와 점탄성 유체 모두 이전 연구[8,9]와 동일한 조건에서 해석을 수행하였고, 이전 연구의 경향을 잘 따라가는 것을 확인하였다. 또한 Fig. 3의 격자 해상도에 따른 차이를 살펴보면 100×100 격자에서는 최대 퍼짐 이후의 수축 구간에서 다소 차이가 나타났지만, 200×200과 400×400 격자에서는 충돌 초기의 급격한 퍼짐 구간과 최대 퍼짐 부근의 거동이 거의 유사하게 나타났으며 격자에 대한 수렴성을 확인할 수 있다. 본 연구에서는 수치적 정확성을 최대한 확보하기 위해 400×400 격자를 사용하여 이후 해석을 진행하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310310/images/jkscfe_2026_313_163_F3.jpg
Fig. 3.

Grid convergence and validation for viscoelastic droplet impact compared with Newtonian case

3.2 Reynolds 수 및 점도비에 따른 최대 퍼짐 거동

액적 충돌에서 최대 퍼짐은 충돌 관성에 의해 액적이 퍼지려는 효과와 점성 및 탄성 응력이 그 퍼짐을 억제하거나 조절하는 효과 사이의 균형으로 결정된다. 따라서 점탄성 액적의 최대 퍼짐 거동을 해석하기 위해서는 먼저 퍼짐 과정의 가장 기본적인 인자인 관성력과 점성력의 상대적 크기를 분석할 필요가 있다. Reynolds 수는 이러한 균형을 나타내는 대표적인 무차원수이다. 뉴턴 유체 액적의 충돌에서는 점성 지배 영역에서의 무차원 최대 퍼짐 직경이 Dmax/D0=CRe1/5 에 비례하는 것으로 알려져 있다[2]. 그러나 점탄성 액적의 경우 전체 점도뿐만 아니라 점탄성 응력이 퍼짐 거동에 추가적으로 영향을 미친다. 따라서 기존 뉴턴 유체의 관계가 그대로 성립하는지, 액체의 점도비 𝛽가 최대 퍼짐 직경에 어떠한 영향을 미치는지에 대한 정량적 검토가 필요하다. 이에 본 절에서는 Reynolds 수와 𝛽를 동시에 변화시키며 점탄성 액적의 최대 퍼짐 거동을 분석하였다.

Fig. 4는 Reynolds 수 변화에 따른 최대 퍼짐 직경 변화를 𝛽=1.0(초록색 네모), 𝛽=0.6(파란색 세모), 𝛽=0.2(빨간색 다이아몬드)일 때 나타낸 결과이다. 해석 조건은 3.1절과 동일하게 D0=0.02 m, U=1 m/s, H=0.03 m에서 진행하였다. 𝛽=1.0은 ηp=0이므로 점탄성 응력의 효과가 사라진다(뉴턴 유체에 해당됨). 반면 𝛽가 감소할수록 전체 점도에서 ηp가 차지하는 비율이 증가한다. 𝛽=1.0 조건에서는 기존 뉴턴 유체에 대한 관계식(검정색 점선)과 일치하는 경향을 보였으며, 𝛽가 감소함에 따라 동일한 Reynolds 수에서 최대 퍼짐 직경이 전체적으로 증가하였다. 다만 Reynolds 수에 대한 증가 기울기는 세 조건 모두에서 뉴턴 유체의 관계식 경향을 그대로 따른다. 즉, 𝛽의 변화는 관계식에서의 비례 상수만 변화시키고, Reynolds 수 변화에 따른 증가 경향성은 그대로 유지되는 형태로 나타났다. 이는 점탄성 액적의 최대 퍼짐 거동이 기본적으로 관성-점성 균형에 의해 지배되어 뉴턴 유체와 동일한 Reynolds 수 의존성을 따르되, 점도 성분비에 따라 최대 퍼짐 직경의 크기가 달라질 수 있음을 보여준다. 따라서 점탄성 액적의 최대 퍼짐 직경은 Reynolds 수만으로 일반화하기 어렵고, 𝛽를 함께 고려해야 함을 확인하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310310/images/jkscfe_2026_313_163_F4.jpg
Fig. 4.

Effect of Reynolds number and viscosity ratio 𝛽 on the maximum spreading diameter

3.3 점도비 𝛽 변화에 따른 최대 퍼짐 거동

앞 절에서는 𝛽와 Reynolds 수가 액적 최대 퍼짐 직경에 동시에 영향을 미치는 것을 확인하였다. 그러나 이 경우 전체 점도의 변화와 점도 성분비의 변화가 동시에 포함되므로, 점탄성 유체를 구성하는 각 점도 성분의 영향을 명확히 구분하기 어렵다. 이러한 분리 해석은 점탄성 액적의 최대 퍼짐 변화가 단순한 전체 점도 변화에 의한 것인지, 또는 ηp와 ηs의 상대적 기여에 의한 것인지를 구분하기 위해 필요하다. 따라서 본 절에서는 Oldroyd-B 유체의 전체 점도η0=ηs+ηp를 4로 일정하게 유지한 상태에서 점도비 𝛽를 변화시켜, 전체 점도의 변화 없이 ηp와 ηs의 상대적 비율이 최대 퍼짐 거동에 미치는 영향을 분석하였다.

Fig. 5는 전체 점도를 η0=ηs+ηp=4로 일정하게 유지한 상태에서 점도비 𝛽의 변화에 따른 Oldroyd-B 유체의 최대 퍼짐 직경을 나타낸다. 𝛽가 증가할수록 전체 점도에서 뉴턴 점도 성분 ηs의 비중은 증가하고 점탄성 기여 점도 ηp의 비중은 감소한다. 계산 결과, 𝛽가 증가함에 따라 최대 퍼짐 직경은 전반적으로 감소하였다. 또한 𝛽=1에서는 ηp=0이 되어 점탄성 응력이 사라지므로 Oldroyd-B 모델은 뉴턴 유체의 조건으로 수렴한다.

특히 𝛽 증가에 따른 최대 퍼짐 직경의 변화율은 전 구간에서 일정하지 않았으며, 𝛽≈0.5 부근을 기준으로 경향이 달라지는 특징이 나타났다. 𝛽<0.5에서는 𝛽 증가에 따라 최대 퍼짐 직경이 비교적 크게 감소한 반면, 𝛽>0.5에서는 그 변화가 상대적으로 완만하게 나타났다. 𝛽=0.5는 ηs=ηp=2로 두 점도 성분의 크기가 동일해지는 조건으로, 본 결과는 전체 점도가 동일하더라도 뉴턴 점도 성분과 점탄성 기여 점도의 상대적 구성에 따라 최대 퍼짐 거동의 변화 정도가 달라질 수 있음을 보여준다. 이러한 변화는 점도 성분의 상대적 구성과 관련된 것으로 보이며, 구체적인 물리적 원인은 향후 점탄성 응력 및 에너지 거동에 대한 추가 분석을 통해 규명하고자 한다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310310/images/jkscfe_2026_313_163_F5.jpg
Fig. 5.

Effect of solvent viscosity ratio 𝛽 on maximum spreading diameter at fixed total viscosity

결과적으로 Oldroyd-B 유체의 전체 점도를 일정하게 유지하더라도 ηs와 ηp의 상대적 구성 변화에 따라 최대 퍼짐 직경이 달라지는 것을 확인하였다. 이는 점탄성 액적의 최대 퍼짐 거동을 전체 점도만으로 설명하기 어렵고, 전체 점도를 구성하는 뉴턴 점도 성분과 점탄성 기여 점도의 상대적 비율과 함께 고려할 필요가 있음을 보여준다.

3.4 Weissenberg 수에 따른 최대 퍼짐 거동

앞 절에서는 전체 점도를 고정한 상태에서 𝛽를 변화시켜, 점탄성 기여 점도의 상대적 비율이 최대 퍼짐 거동에 미치는 영향을 분석하였다. 그러나 점탄성 액적의 퍼짐 거동은 점탄성 응력의 상대적 크기뿐만 아니라, 충돌 과정에서 형성된 응력이 유동 중 얼마나 빠르게 완화되는지에도 영향을 받을 수 있다. Oldroyd-B 유체에서 이러한 시간 의존적 특성은 완화 시간(𝜆)으로 나타나며, 이는 Weissenberg 수(Wi)를 통해 유동 시간 척도와 비교할 수 있다.

Fig. 6은 전체 점도 η0=4로 고정한 상태에서 Weissenberg 수의 변화에 따른 최대 퍼짐 직경 변화(Dmax/D0)를 나타낸다. Weissenberg 수의 변화는 다른 충돌 조건(초기 속도, 초기 직경)은 유지한 상태로 완화 시간(𝜆)만을 변화시켜 나타내었다. 계산 결과, Weissenberg 수가 증가함에 따라 최대 퍼짐 직경이 증가하는 경향이 나타났다. 특히 𝛽가 작은 조건일수록 Weissenberg 수의 증가에 따른 최대 퍼짐 직경의 증가가 상대적으로 크게 나타났으며, 𝛽가 1에 가까워질수록 그 변화는 감소하였다. 이러한 결과로부터 점도 성분비와 완화 시간의 변화가 최대 퍼짐 거동에 함께 영향을 미치는 경향을 확인할 수 있다. 𝛽가 작을수록 전체 점도에서 점탄성 기여 점도의 비중이 증가하고, 본 계산에서는 이러한 조건에서 Weissenberg 수 변화에 따른 최대 퍼짐 직경의 차이가 상대적으로 크게 나타났다. 다만 이러한 변화의 구체적인 물리적 원인을 규명하기 위해서는 점탄성 응력 및 에너지 거동에 대한 추가적인 분석이 필요하다.

결과적으로 Fig. 6은 점탄성 액적의 최대 퍼짐 거동이 점도 성분비뿐만 아니라 응력의 시간적 완화 특성에도 영향을 받는다는 점을 보여준다. 따라서 점탄성 액적의 최대 퍼짐 거동을 분석할 때 점도 성분의 상대적 구성비 𝛽와 완화 특성을 함께 고려할 필요가 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-03/N0500310310/images/jkscfe_2026_313_163_F6.jpg
Fig. 6.

Effect of Weissenberg number on maximum spreading diameter at fixed total viscosity

4. 결 론

본 연구에서는 고체 표면에 충돌하는 점탄성 액적의 최대 퍼짐 거동을 분석하기 위해 LCRM 기반 액적 충돌 해석 체계에 Oldroyd-B 구성방정식을 적용하였다. 이를 통해 Reynolds 수, 점도 성분비 및 Weissenberg 수를 체계적으로 변화시키며 여러 점탄성 조건에서의 최대 퍼짐 직경을 비교하고, 각 물성 조건이 최대 퍼짐 거동에 미치는 영향을 종합적으로 분석하였다.

연구 결과, 점탄성 액적의 최대 퍼짐 직경은 Reynolds 수뿐만 아니라 점도 성분의 상대적 구성과 완화 특성에 따라서도 변화하는 것으로 나타났다. 𝛽=1인 조건에서는 점탄성 기여 점도가 0이 되어 뉴턴 유체와 유사한 퍼짐 경향을 보였으며, Reynolds 수가 증가함에 따라 최대 퍼짐 직경이 증가하였다. 또한 𝛽가 감소할수록 동일한 Reynolds 수에서 최대 퍼짐 직경이 증가하는 경향이 나타났다. 전체 점도를 일정하게 유지한 조건에서도 뉴턴 점도 성분과 점탄성 기여 점도의 상대적 구성 변화에 따라 최대 퍼짐 직경이 달라져, 점탄성 액적의 퍼짐 거동을 전체 점도만으로 설명하기 어렵다는 것을 확인하였다. Weissenberg 수의 경우 그 값이 증가함에 따라 최대 퍼짐 직경이 증가하였으며, 이러한 변화는 점탄성 기여 점도의 비중이 큰 조건에서 상대적으로 뚜렷하게 나타났다.

결과적으로 본 연구에서는 여러 점탄성 물성 조건에 따른 최대 퍼짐 직경을 직접 비교함으로써 Reynolds 수, 점도 성분비 및 Weissenberg 수를 함께 고려할 필요가 있음을 확인하였다. 다만 현재의 최대 퍼짐 직경 결과만으로 각 조건에서 나타난 변화의 세부적인 물리적 메커니즘을 단정하기에는 한계가 있다. 향후에는 점탄성 응력의 시간적·공간적 분포와 점성 소산 및 탄성에너지 거동을 추가적으로 분석하고, 초기 액적 직경, 충돌 속도 및 젖음성 등 다양한 조건으로 연구 범위를 확장하여 점탄성 액적의 충돌 거동을 보다 구체적으로 규명하고자 한다.

Acknowledgements

본 연구는 정부(과학기술정보통신부)의 재원으로 한국연구재단(No. RS-2025-02302984)의 지원을 받아 수행된 연구입니다.

References

1

2024, Shah, P. and Driscoll, M.M., “Drop Impact Dynamics of Complex Fluids: A Review,” Soft Matter, Vol.20, pp.4839-4858.

10.1039/D4SM00145A
2

2004, Clanet, C., Béguin, C., Richard, D. and Quéré, D., “Maximal Deformation of an Impacting Drop,” J. Fluid Mech., Vol.517, pp.199-208.

10.1017/S0022112004000904
3

2000, Bergeron, V., Bonn, D., Martin, J.Y. and Vovelle, L., “Controlling Droplet Deposition with Polymer Additives,” Nature, Vol.405, pp.772-775.

10.1038/35015525
4

2007, Bartolo, D., Boudaoud, A., Narcy, G. and Bonn, D., “Dynamics of Non-Newtonian Droplets,” Phys. Rev. Lett., Vol.99, 174502.

10.1103/PhysRevLett.99.174502
5

2001, Pillapakkam, S.B. and Singh, P., “A Level-Set Method for Computing Solutions to Viscoelastic Two-Phase Flow,” J. Comput. Phys., Vol.174, pp.552-578.

10.1006/jcph.2001.6927
6

2005, Chinyoka, T., Renardy, Y.Y., Renardy, M. and Khismatullin, D.B., “Two-Dimensional Study of Drop Deformation under Simple Shear for Oldroyd-B Liquids,” J. Non-Newton. Fluid Mech., Vol.130, pp.45-56.

10.1016/j.jnnfm.2005.07.005
7

2012, Oishi, C.M., Martins, F.P., Tomé, M.F. and Alves, M.A., “Numerical Simulation of Drop Impact and Jet Buckling Problems Using the eXtended Pom-Pom Model,” J. Non-Newton. Fluid Mech., Vol.169-170, pp.91-103.

10.1016/j.jnnfm.2011.12.001
8

2016, Figueiredo, R.A., Oishi, C.M., Afonso, A.M., Tasso, I.V.M. and Cuminato, J.A., “A Two-Phase Solver for Complex Fluids: Studies of the Weissenberg Effect,” Int. J. Multiph. Flow, Vol.84, pp.98-115.

10.1016/j.ijmultiphaseflow.2016.04.014
9

2012, Xu, X., Ouyang, J., Jiang, T. and Li, Q., “Numerical Simulation of 3D-Unsteady Viscoelastic Free Surface Flows by Improved Smoothed Particle Hydrodynamics Method,” J. Non-Newton. Fluid Mech., Vol.177-178, pp.109-120.

10.1016/j.jnnfm.2012.04.006
10

2018, Viezel, C., Tomé, M.F., Pinho, F.T. and McKee, S., “A Numerical Technique for Solving the Oldroyd-B Model for the Whole Range of Viscosity Ratios,” Proceedings of the 6th European Conference on Computational Mechanics (ECCM 6) and 7th European Conference on Computational Fluid Dynamics (ECFD 7), Glasgow, UK, pp.2913-2926.

11

2014, Figueiredo, R.A., Oishi, C.M., Cuminato, J.A., Azevedo, J.C., Afonso, A.M. and Alves, M.A., “Numerical Investigation of Three Dimensional Viscoelastic Free Surface Flows: Impacting Drop Problem,” Proceedings of the 6th European Conference on Computational Fluid Dynamics (ECFD VI), pp.5368-5380.

12

2009, Shin, S. and Juric, D., “A Hybrid Interface Method for Three-Dimensional Multiphase Flows Based on Front Tracking and Level Set Techniques,” Int. J. Numer. Methods Fluids, Vol.60, pp.753-778.

10.1002/fld.1912
13

1950, Oldroyd, J.G., “On the Formulation of Rheological Equations of State,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, Vol.200, pp.523-541.

10.1098/rspa.1950.0035
14

2018, Shin, S., Chergui, J. and Juric, D., “Direct Simulation of Multiphase Flows with Modeling of Dynamic Interface Contact Angle,” Theor. Comput. Fluid Dyn., Vol.32, pp.655-687.

10.1007/s00162-018-0470-4
15

1968, Chorin, A.J., “Numerical Solution of the Navier-Stokes Equations,” Math. Comput., Vol.22, pp.745-762.

10.1090/S0025-5718-1968-0242392-2
페이지 상단으로 이동하기