비베 함수

편집 역사 토론
계산화학 시뮬레이션 수치해석 마지막 수정: 2026-09-03 04:27:36

1. 개요[편집]

비베 함수
Wiebe (Vibe) Function
제안이반 이바노비치 비베 (И. И. Вибе), 1956
정체기연질량분율 xb(θ) 의 경험적 S 곡선
계수효율계수 a, 형상계수 m, 연소시작 θs, 지속기간 Δθ
수학적 정체와이불 누적분포함수 (형상 m+1)
쓰는 곳1D 엔진 사이클 코드, 단일영역 열역학 해석, 제어 지향 모델
못 하는 것난류·화염 전파·화학의 예측

연소를 물리로 풀 수 없으면, 연소를 곡선으로 그려 놓고 나머지를 물리로 푼다.

비베 함수(Wiebe function, 유럽 문헌에서는 Vibe function)는 실린더 안에서 연료가 타 들어간 비율을 크랭크각의 함수로 표현하는 경험식으로, 연소의 상세 물리를 계수 네 개로 압축해 열역학 사이클 계산에 꽂아 넣기 위한 도구다.

xb(θ)=1exp ⁣[a(θθsΔθ)m+1]x_b(\theta) = 1 - \exp\!\left[-a\left(\frac{\theta - \theta_s}{\Delta\theta}\right)^{m+1}\right]

xbx_b기연질량분율(mass fraction burned, MFB), θs\theta_s 는 연소 시작 크랭크각, Δθ\Delta\theta 는 연소 지속기간, aa 는 효율계수, mm 은 형상계수다. θ<θs\theta < \theta_s 에서는 xb=0x_b = 0, θ>θs+Δθ\theta > \theta_s + \Delta\theta 에서는 xb1x_b \approx 1 로 잘라 쓴다.

소련의 내연기관 학자 이반 이바노비치 비베(1902~1969)가 1956년에 제시했고, 이후 70년간 엔진 1D 사이클 해석의 사실상 표준 열방출 모델로 살아남았다. 이유는 단순하다. 계수 네 개면 실측 압력선도를 놀랄 만큼 잘 재현하고, 미분이 해석적으로 떨어지며, 계산 비용이 사실상 0 이다. 대신 대가도 명확하다 — 이 함수 안에는 난류도 화염도 화학도 없다.

2. 계수의 뜻[편집]

효율계수 aa 는 “지속기간이 끝났을 때 몇 %까지 태울 것인가”만 결정한다. θ=θs+Δθ\theta = \theta_s + \Delta\theta 를 넣으면 xb=1eax_b = 1 - e^{-a} 이므로 거꾸로 풀면

a=ln(1xb,end)a = -\ln\left(1 - x_{b,\mathrm{end}}\right)
  • xb,end=99.9%x_{b,\mathrm{end}} = 99.9\%a=ln1000=6.908a = \ln 1000 = 6.908
  • xb,end=99%x_{b,\mathrm{end}} = 99\%a=ln100=4.605a = \ln 100 = 4.605
  • xb,end=99.3%x_{b,\mathrm{end}} = 99.3\%a=5a = 5 (하이우드의 교과서가 쓰는 값)

널리 인용되는 6.908 은 99.9% 에 대응하는 값이고, 하이우드 계열 문헌의 a=5a = 5 는 99.3% 기준이다. 즉 aa 는 물리량이 아니라 “Δθ 를 무엇의 지속기간으로 정의했는가”라는 약속에 불과하다. 서로 다른 문헌의 Δθ 를 비교하려면 aa 를 먼저 맞춰야 한다.1

형상계수 mm 이 곡선의 성격을 결정한다. 미분형(= 정규화된 열방출률)을 보면 바로 드러난다.

dxbdθ=a(m+1)Δθ(θθsΔθ)mexp ⁣[a(θθsΔθ)m+1]\frac{dx_b}{d\theta} = \frac{a(m+1)}{\Delta\theta}\left(\frac{\theta-\theta_s}{\Delta\theta}\right)^{m} \exp\!\left[-a\left(\frac{\theta-\theta_s}{\Delta\theta}\right)^{m+1}\right]

앞의 멱함수가 초기 상승을, 지수함수가 후기 소멸을 담당한다. 정규화 시간 τ=(θθs)/Δθ\tau = (\theta-\theta_s)/\Delta\theta 로 두고 최댓값 조건을 풀면 열방출률의 피크 위치가 닫힌 형태로 나온다.

τ=[ma(m+1)]1/(m+1)\tau^* = \left[\frac{m}{a(m+1)}\right]^{1/(m+1)}

피크에서의 기연분율은 더 예쁘게 정리된다. aτm+1=m/(m+1)a\tau^{*\,m+1} = m/(m+1) 이므로

xb(τ)=1exp ⁣(mm+1)x_b(\tau^*) = 1 - \exp\!\left(-\frac{m}{m+1}\right)

열방출률이 최대가 되는 시점의 연소 진행도는 aa 와 무관하게 mm 만으로 결정된다. a=6.908a = 6.908 기준으로 정리하면 이렇다.

m열방출률 피크 위치 τ*그때 기연분율성격
00 (단조 감소)0초기 집중형, 비물리적
10.270.39전반부 집중
20.460.49거의 대칭 — SI 엔진의 국룰
40.650.55후반부 집중

mm 이 커지면 피크가 뒤로 밀리고 곡선이 대칭에 가까워진다. SI 엔진 예혼합 연소는 화염이 커지다가 벽에 닿아 잦아드는 거동이라 대략 대칭이고, 그래서 m2m \approx 2, Δθ4060\Delta\theta \approx 40 \sim 60 °CA 가 출발점으로 통한다. 디젤은 초반에 예혼합 연소가 확 터지므로 이 한 개의 곡선으로는 표현이 안 되고, 뒤에 나오는 이중 비베로 간다.

수학적으로 보면 xbx_b형상모수 m+1m+1, 척도모수 a1/(m+1)a^{-1/(m+1)} 의 와이불 누적분포함수와 정확히 같다. 비베가 연쇄반응 논리로 유도했다고 알려져 있지만, 실용적으로는 “0에서 1로 올라가는 S 곡선 중 계수 두 개로 비대칭성까지 조절되는 가장 값싼 형태”라고 이해하는 게 정직하다. 물리적 유도의 권위를 빌릴 필요가 없는 함수다.

3. 실린더 압력을 적분하는 절차[편집]

비베 함수 자체는 압력을 주지 않는다. 압력은 단일영역(single-zone) 열역학 1법칙에서 나온다. 실린더 내용물을 균일한 이상기체 한 덩어리로 보고 비열비 γ\gamma 를 상수로 두면

dpdθ=γ1V(θ)(dQdθdQhtdθ)γpVdVdθ\frac{dp}{d\theta} = \frac{\gamma - 1}{V(\theta)}\left(\frac{dQ}{d\theta} - \frac{dQ_{ht}}{d\theta}\right) - \gamma\,\frac{p}{V}\,\frac{dV}{d\theta}

가 된다. 오른쪽 첫 항이 연소에 의한 압력 상승, 둘째 항이 피스톤 운동에 의한 압축·팽창이다. 여기에 비베가 들어가는 자리는 딱 하나다.

dQdθ=Qindxbdθ,Qin=mfLHVηcomb\frac{dQ}{d\theta} = Q_{\mathrm{in}}\,\frac{dx_b}{d\theta}, \qquad Q_{\mathrm{in}} = m_f\,\mathrm{LHV}\,\eta_{\mathrm{comb}}

체적은 슬라이더-크랭크 기구학에서 정확하게 온다. 압축비 rcr_c, 커넥팅로드/크랭크반경 비 R=l/aR = l/a 로 쓰면

V(θ)Vc=1+rc12(R+1cosθR2sin2θ)\frac{V(\theta)}{V_c} = 1 + \frac{r_c - 1}{2}\left(R + 1 - \cos\theta - \sqrt{R^2 - \sin^2\theta}\right)

이고 이것을 미분해 dV/dθdV/d\theta 를 넣는다. dQht/dθdQ_{ht}/d\theta 는 벽면 열손실로, 보통 보슈니(Woschni) 상관식으로 순간 열전달계수를 주고 hA(TTw)h A (T - T_w) 를 쓴다. 열손실을 빼먹으면 압력이 눈에 띄게 과대예측된다 — 비베로 맞춘 것처럼 보이던 계수가 열손실 모델을 켜자마자 전부 달라지는 이유이고, 계수 피팅과 열전달 모델은 항상 짝으로 검증해야 한다.

이렇게 하면 p(θ)p(\theta) 를 크랭크각에 대한 상미분방정식 하나로 룽게-쿠타법으로 적분할 수 있다. 초기조건은 흡기밸브 닫힘(IVC) 시점의 압력·온도·조성이고, 그 조성에 잔류가스와 배기가스 재순환 분율이 들어간다.

비베 연소율 x_b(θ) 하나를 단일영역 에너지식 dp/dθ = (γ−1)/V·dQ/dθ − γ(p/V)dV/dθ 에 넣고 RK4 로 0.25° 씩 −180°…+180° 를 적분한다(r=10, γ=1.35, Q_tot=1800 J, 열손실·블로바이 없이 압축·팽창만). 점화시기 θ_s 를 옮기면 IMEP 가 최댓값을 갖는 MBT 가 나타나며, Δθ=40° 에서 실측 MBT θ_s = −18.8°, IMEP 19.68 bar, η = 0.547 이다. 아래 스윕 곡선이 θ_s → IMEP 를 21점 재적분해 그 꼭대기를 표시한다.

4. IMEP — 선도를 면적으로 바꾸기[편집]

압력 이력이 나오면 일은 적분이다. 도시평균유효압력(IMEP)은 사이클당 도시일을 행정체적으로 나눈 값이다.

IMEP=1VdpdV\mathrm{IMEP} = \frac{1}{V_d}\oint p\,dV

여기서 적분 구간을 어디로 잡느냐가 실무에서 사고를 부른다.

  • 총 IMEP(gross): 압축·팽창 행정만, 즉 360° 구간. 연소 품질만 보고 싶을 때.
  • 정미 IMEP(net): 720° 전체. 흡배기 행정의 펌핑 루프까지 포함한다.
  • 둘의 차이가 PMEP(펌핑 평균유효압력)이고, 정미 IMEP 에서 마찰 FMEP 를 빼면 축토크에 해당하는 BMEP 가 된다.

pp-VV 선도에서 압축·팽창 루프는 시계 방향(양의 일), 부분부하 흡배기 루프는 반시계 방향(음의 일)으로 그려진다. 가변 밸브 타이밍이나 밀러 사이클이 효율을 올린다고 할 때 실제로 줄이는 것이 이 반시계 루프의 면적이며, 비베 기반 사이클 모델이 그 효과를 정량화하는 최소 도구다.

도시열효율은 ηi=pdV/(mfLHV)\eta_i = \oint p\,dV / (m_f \mathrm{LHV}) 로 바로 나온다. 이 단일영역 모델로도 점화시기 스윕을 돌려 최적 진각(MBT)을 찾는 것이 가능한데, 지각하면 팽창 일이 줄고 진각하면 압축 중 열방출이 늘어 손실이 커지는 두 효과의 균형이 자동으로 잡히기 때문이다.

5. 이중 비베 — 디젤의 두 얼굴[편집]

디젤 연소는 성격이 다른 두 단계로 나뉜다. 착화 지연 동안 섞여 있던 연료가 한꺼번에 타는 예혼합 연소(급격한 스파이크)와, 그 뒤 분사 속도가 열방출을 지배하는 확산 연소(넓고 완만한 꼬리)다. 단일 비베로는 이 둘을 동시에 못 그린다. 그래서 두 개를 선형결합한다.

xb(θ)=βxb,1(θ; θs,1,Δθ1,m1)+(1β)xb,2(θ; θs,2,Δθ2,m2)x_b(\theta) = \beta\, x_{b,1}(\theta;\ \theta_{s,1}, \Delta\theta_1, m_1) + (1-\beta)\, x_{b,2}(\theta;\ \theta_{s,2}, \Delta\theta_2, m_2)

β\beta 는 예혼합 연소가 태우는 연료 분율이다. 예혼합 항은 mm 을 작게(0.5~1) 잡아 앞으로 쏠린 좁은 봉우리를, 확산 항은 mm 을 크게 잡아 완만한 봉우리를 만든다. 계수가 7 개로 늘어나므로 자유도가 넘치기 시작하고, 서로 다른 계수 조합이 거의 같은 압력선도를 주는 식별 불가능성이 본격적으로 문제가 된다(역문제에서 익숙한 그 구조다).

다단 분사(파일럿 + 메인 + 포스트)를 쓰는 현대 디젤은 분사 펄스마다 비베를 하나씩 붙여 3~5 중 비베까지 가고, 이 지점부터 “모델”이라기보다 압력선도를 재생하는 압축 포맷에 가까워진다. 예측력이 필요하면 결국 분무·착화 물리가 들어간 모델(연소 시뮬레이션 계열)로 갈아타야 한다.

6. 역산 — 측정 압력에서 계수를 뽑기[편집]

실무에서 비베 계수는 손으로 고르는 것이 아니라 실측 압력선도에서 역산한다. 절차는 정형화돼 있다.

1단계, 겉보기 열방출률(ROHR). 위의 1법칙을 압력 쪽에서 뒤집는다.

dQndθ=γγ1pdVdθ+1γ1Vdpdθ\frac{dQ_n}{d\theta} = \frac{\gamma}{\gamma-1}\,p\,\frac{dV}{d\theta} + \frac{1}{\gamma-1}\,V\,\frac{dp}{d\theta}

QnQ_n 을 “정미(net)” 열방출이라 부르는 것은 벽면 열손실이 이미 빠진 값이기 때문이다. 열손실 모델을 따로 더하면 총 열방출이 되지만, 계수 피팅에는 정미 쪽이 오히려 일관적이다.

2단계, 기연분율로 정규화. xb(θ)=θsθdQn/totaldQnx_b(\theta) = \int_{\theta_s}^{\theta} dQ_n / \int_{\mathrm{total}} dQ_n. 압력 데이터에서 폴리트로픽 지수만으로 MFB 를 뽑는 라스바일러-위드로(1938) 방법도 여전히 널리 쓰인다 — 계산이 가볍고 γ\gamma 가정에 덜 민감해 연소해석 장비의 실시간 표시에 들어가 있다.

3단계, 계수 맞추기. xbx_b 곡선에 비베를 비선형 최소제곱으로 피팅한다. 그런데 닫힌 형태로 푸는 지름길이 있다. 식을 두 번 로그 취하면

ln[ln(1xb)]=lna+(m+1)lnθθsΔθ\ln\left[-\ln(1-x_b)\right] = \ln a + (m+1)\ln\frac{\theta-\theta_s}{\Delta\theta}

이므로 ln[ln(1xb)]\ln[-\ln(1-x_b)]ln(θθs)\ln(\theta-\theta_s)직선이고, 기울기가 m+1m+1, 절편에서 Δθ\Delta\theta 가 나온다. 이 좌표계에서 최소자승법 한 번이면 끝이다. 더 거칠게는 10% 와 90% 연소각 두 점만 써도 된다.

m+1=ln ⁣[ln0.9/ln0.1]ln ⁣[(θ10θs)/(θ90θs)]m + 1 = \frac{\ln\!\left[\ln 0.9 / \ln 0.1\right]}{\ln\!\left[(\theta_{10}-\theta_s)/(\theta_{90}-\theta_s)\right]}

Δθ\Delta\theta 는 두 점 중 하나를 다시 대입해 얻는다. 현장에서 “MFB10-90 구간이 몇 도”라는 지표를 그렇게 자주 말하는 이유가 이것이다 — 그 두 숫자가 사실상 비베 계수 그 자체다.

함정 목록. 이 역산이 조용히 틀리는 경로가 많다.

  • 압력 절대값 보정(pegging). 압전 센서는 상대압만 준다. 기준점을 IVC 흡기압으로 잡거나 압축 구간의 폴리트로픽 피팅으로 잡는데, 0.1 bar 만 어긋나도 초기 열방출률이 통째로 뒤집힌다.
  • 상사점 위상(TDC phasing). 압력 데이터의 크랭크각 0° 가 실제 상사점과 어긋나면 IMEP 가 몇 %씩 흔들린다. 1 °CA 오차가 10% 수준까지 갈 수 있어, 모티어링(motored) 압력의 피크 위치로 열손실 보정을 거쳐 상사점을 따로 결정하는 것이 표준 절차다.
  • γ\gamma 를 상수로 두는 것. 온도 700 K 와 2500 K 에서 γ\gamma 는 1.35 와 1.25 쯤으로 다르다. 상수 γ\gamma 는 총 열방출량을 수 % 틀리게 하고, 그 오차가 ηcomb\eta_{\mathrm{comb}} 로 흡수되어 보이지 않게 숨는다.
  • 평균의 함정. 사이클 변동이 있으니 100~300 사이클을 평균하는데, 평균 압력선도에서 뽑은 열방출률은 열방출률의 평균이 아니다. 변동이 큰 조건(대량 EGR·희박)에서는 사이클별로 역산해 계수의 분포를 보는 것이 맞다.

7. 1D 사이클 코드에서의 위치와 한계[편집]

GT-Power·Ricardo WAVE·AVL BOOST 같은 1D 엔진 코드에서 실린더는 배관망 한가운데 붙은 가변체적 용기다. 흡배기관은 1차원 압축성 유동으로, 밸브는 유출계수를 갖는 오리피스로 풀고, 실린더 안의 연소는 비베가 담당한다. 여기서 나온 p(θ)p(\theta) 가 다시 다음 것들의 입력이 된다.

  • 노킹 판정 — 말단가스의 온도·압력 이력으로 착화 적분(Livengood-Wu)을 돌린다.
  • 질소산화물 예측 — 기연영역 온도에 젤도비치 기구를 얹어 NO 를 적분한다. 이때는 단일영역을 미연·기연 2구역으로 쪼개야 한다(단일영역 평균온도로는 NOx 가 자릿수로 틀린다).
  • 과급기 정합 — 배기 엔탈피와 맥동이 터빈 일을 결정한다.

그리고 여기서 한계가 정직하게 드러난다. 비베 계수 a,m,Δθ,θsa, m, \Delta\theta, \theta_s 안에는 회전속도·부하·점화시기·희석률·텀블 강도·연료 종류가 전부 뭉개져 들어가 있다. 그러므로

  • 맞춘 조건 밖으로 외삽할 수 없다. EGR 을 10% 에서 25% 로 올리면 실제로는 층류 화염 속도가 떨어져 Δθ 가 길어지는데, 비베는 그것을 스스로 알지 못한다. 계수를 운전점 맵으로 만들어 보간하는 것이 관행이고, 그 맵 바깥은 그냥 미지다.
  • 사이클 변동을 표현할 수 없다. 계수가 결정론적이니 매 사이클 같은 압력이 나온다. COV 를 보려면 계수에 확률분포를 얹는 편법을 쓴다.
  • 화학이 지배하는 연소는 원리적으로 불가능하다. HCCI 처럼 열방출 시점이 반응속도로 정해지는 연소에서는 비베로 맞춘 계수가 조건을 조금만 바꿔도 무너진다. 이 영역은 다구역 모델이나 상세 화학 직접 결합으로 가야 한다.

그래서 예측이 필요한 자리에는 난류 연소 모델을 쓴다. 블리자드-켁 계열의 2구역 와동 유입(eddy entrainment) 모델은 화염면이 난류에 의해 삼켜 들이는 질량과 층류 화염 속도로 타 들어가는 질량을 따로 적분해 연소율을 계산한다. 입력이 계수 네 개가 아니라 난류강도·화염 반경·SLS_L 상관식이므로, 회전속도와 희석을 바꿔도 스스로 반응한다. GT-Power 의 SITurb 계열이 이 부류다.

정리하면 비베 함수의 위치는 명확하다. 연소를 설명하는 모델이 아니라, 연소를 통과시켜 그 뒤의 열역학·가스동역학·후처리 물리를 풀게 해 주는 통로다. 이 함수에 연소 물리를 기대하지 않는 한 여전히 가장 효율 좋은 도구고, 기대하는 순간 가장 위험한 도구가 된다.23

8. 관련 문서[편집]

9. Footnotes[편집]

  1. “Δθ = 45 °CA” 라는 문장 하나로 두 논문을 비교하다가, 한쪽은 a=5a=5(99.3%)이고 다른 쪽은 a=6.908a=6.908(99.9%)이라 실제 연소 지속기간이 10% 이상 다르다는 걸 뒤늦게 깨닫는 일이 흔하다. 표 각주를 읽는 습관이 여기서 밥값을 한다.

  2. 비베 계수 세 개로 실측 압력선도를 0.1 bar 이내로 재현해 놓고 “우리 연소 모델이 검증됐다”고 쓴 보고서를 본다면, 그건 검증이 아니라 곡선 맞추기가 성공했다는 보고다. V&V 정신에서 이건 validation 이 아니라 calibration 이다.

  3. 그래도 이 함수가 70년을 버틴 이유는 분명하다. 노트북에서 전체 운전맵 수천 점을 몇 분에 돌릴 수 있는 모델은 지금도 이것뿐이고, 제어기 개발과 초기 설계 탐색은 정확도보다 그 속도를 필요로 한다. “일단 돌려”가 정당한 몇 안 되는 자리다.