王 飛，李兆敏，李松巖，杜慶軍（.中國(guó)石油大學(xué)（華東）石油工程學(xué)院，山東青島 66580；.中國(guó)石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，山東青島 66580）

文章編號(hào)：1001?246X（2015）01?0058?07

泡沫混排解堵數(shù)學(xué)模型及數(shù)值模擬

王 飛1，李兆敏1，李松巖1，杜慶軍2
（1.中國(guó)石油大學(xué)（華東）石油工程學(xué)院，山東青島 266580；2.中國(guó)石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，山東青島 266580）

建立泡沫井筒流動(dòng)模型與泡沫地層滲流模型，得到完整的泡沫混排解堵數(shù)學(xué)模型；利用數(shù)值方法對(duì)模型進(jìn)行耦合求解，得到泡沫壓力，質(zhì)量，密度在井筒中的分布及井口井底壓力的變化規(guī)律.在固定井口回壓情況下，得到井底壓差的變化規(guī)律.結(jié)果表明：隨著井深的增加，泡沫壓力增加，質(zhì)量減小，密度增加；固定井口回壓，井底壓差則逐漸變小.

泡沫混排；地層滲流；井筒流動(dòng)；數(shù)學(xué)模型；井底壓差

0 引言

在鉆井，完井，采油等過(guò)程中，油藏近井地帶常受到不同程度的傷害.造成近井地層污染和堵塞，油井產(chǎn)量下降.目前，根據(jù)泡沫流體性質(zhì)提出的泡沫混排解堵技術(shù)能有效地解決此類問(wèn)題.通過(guò)井筒將泡沫注入近井地層，再利用泡沫的可壓縮特性將其迅速放噴，從而把近井地帶堵塞物攜帶排出，達(dá)到解除堵塞、改善近井地帶滲流狀況的目的.

此項(xiàng)技術(shù)在現(xiàn)場(chǎng)生產(chǎn)中得到了成功應(yīng)用［1－2］，實(shí)驗(yàn)研究［3］也有了一定進(jìn)展，但完善的理論還未建立，現(xiàn)有的模型對(duì)于地層滲流未考慮［4－6］或與井筒流動(dòng)過(guò)程分開(kāi)考慮［7］，而數(shù)值模擬常采用相關(guān)軟件［8－9］進(jìn)行，缺乏準(zhǔn)確的結(jié)論及認(rèn)識(shí).本文將注入返排中的井筒管流過(guò)程與泡沫在地層的滲流［10－11］過(guò)程結(jié)合，建立了耦合模型，完善了泡沫混排理論，數(shù)值模擬結(jié)果可以指導(dǎo)現(xiàn)場(chǎng)應(yīng)用.

1 物理模型

混排過(guò)程分為兩個(gè)部分：①泡沫注入部分，通過(guò)油管將泡沫注入近井地層；②泡沫返排部分，將泡沫從地層沿油管放噴（見(jiàn)圖1）.

圖1 泡沫混排示意圖Fig.1 Schematic of foam plug removal

2 數(shù)學(xué)模型

根據(jù)物理模型，泡沫混排數(shù)學(xué)模型分為兩個(gè)部分，泡沫注入部分和泡沫返排部分.每個(gè)部分又分為兩個(gè)過(guò)程：①是泡沫在井筒中的流動(dòng)過(guò)程；②是泡沫在地層中的滲流過(guò)程.建立數(shù)學(xué)模型就是將這兩個(gè)過(guò)程耦合，然后將兩個(gè)部分結(jié)合，得到整個(gè)施工流程上的模型.

2.1 泡沫注入模型建立

2.1.1 井筒管流模型

泡沫在井筒內(nèi)的流動(dòng)基于質(zhì)量守恒方程和動(dòng)量守恒方程［12－13］，假設(shè)

1）泡沫流體處于穩(wěn)定的流動(dòng)狀態(tài)；

2）泡沫流體的可壓縮性完全取決于泡沫內(nèi)氣體的可壓縮性；

3）井筒的溫度按地溫梯度計(jì)算.

根據(jù)以上假設(shè)，可以列出注入時(shí)泡沫流體在井筒中流動(dòng)的數(shù)學(xué)模型.

質(zhì)量守恒方程

動(dòng)量守恒方程

式中，ρg，ρl為泡沫中氣體和液體的密度，kg·m－3；qg，ql為 泡沫中氣體和液體的體積流量，m3·s－1；（d p/d H）f＝/D為泡沫流體的摩阻壓降；（d p/d H）f＝－ρfgsinθ為泡沫流體的重力壓降；（d p/d H）ac為泡沫流體的加速壓降，一般很小，可以忽略不計(jì)；ρf為泡沫密度；Vf為泡沫流體的流速，m·s－1；D為井筒直徑；θ為管路傾角度.

2.1.2 泡沫滲流模型

泡沫地層中滲流過(guò)程可以近似為單相徑向流［14－15］.泡沫具有暫堵分流特性［16］，對(duì)于同一小層，根據(jù)泡沫圈閉氣體和泡沫質(zhì)量守恒可以得到以下方程.

地層內(nèi)泡沫前沿半徑

氣相圈閉擬表皮系數(shù)

地層的吸液能力

泡沫體積流量

氣體體積流量

泡沫質(zhì)量

式中，Rf為注入泡沫前沿半徑，m；Rw為井筒半徑，m；qf為注入地層的泡沫流量，m3·s－1；Γ為注入泡沫質(zhì)量，小數(shù)；?為地層孔隙度，小數(shù)；SP為注入氣相圈閉擬表皮系數(shù)，無(wú)量綱；K，Kf分別為地層原始滲透率和地層注入泡沫后滲透率，D；K/Kf反映了泡沫的穩(wěn)定性，可由實(shí)驗(yàn)確定；pw為注入井底壓力，MPa；pe為油藏邊界壓力，MPa；μf為注入泡沫流體黏度，Pa·s；Re為泄油邊緣半徑，m；S為地層初始表皮系數(shù)，無(wú)量綱；Swc為束縛水飽和度，無(wú)量綱；ql，qg為注入泡沫中液體和氣體體積流量，m3·s－1；Te為返排地層溫度，K；zg為注入氣體壓縮因子，無(wú)量綱；qgsc為注入泡沫中氣體在標(biāo)準(zhǔn)狀態(tài)下的體積流量，m3·s－1；pgsc為標(biāo)準(zhǔn)狀態(tài)下的壓力，0.1 MPa；Tgsc為標(biāo)準(zhǔn)狀態(tài)下的溫度，273.15 K.

2.1.3 邊界條件

在泡沫混排中，注入井口泡沫流體的標(biāo)方流量保持為定值，返排井口回壓保持為定值.

qgsc＝const； ql井口＝const； p井底（t） ＝pw（t）（耦合邊界）．

2.1.4 初始條件

p井底（0） ＝pe．

將井筒流動(dòng)模型與泡沫在地層中的滲流模型耦合得到泡沫注入模型.

2.2 泡沫返排模型建立

1）井筒管流模型

與注入部分相同假設(shè)，同樣可以列出返排時(shí)泡沫流體在井筒中流動(dòng)的數(shù)學(xué)模型.

質(zhì)量守恒方程

動(dòng)量守恒方程

2）泡沫滲流模型

返排時(shí)地層內(nèi)泡沫前沿半徑逐漸減小，根據(jù)注入后半徑得

其中：Qt為注入地層的泡沫量，m3；′為返排地層的泡沫流量，m3·s－1；

與注入部分類似可以列出：

返排氣相圈閉擬表皮系數(shù)

返排地層的吸液能力

返排泡沫體積流量

返排氣體體積流量

返排泡沫質(zhì)量

3）邊界條件

在泡沫混排中，注入井口泡沫流體的標(biāo)方流量保持為定值，返排井口回壓保持為定值.

4）初始條件

將井筒流動(dòng)模型與泡沫在地層中的滲流模型耦合得到泡沫返排模型.

2.3 輔助方程

1）油管溫度按地溫梯度計(jì)算，傾角按井身結(jié)構(gòu)參數(shù)計(jì)算.

T＝T（H）；θ＝θ（H）．

2）假設(shè)注入返排氣液質(zhì)量比不變，即注入標(biāo)方下氣液比與返排標(biāo)方下氣液比相等.

約束條件：井底壓力pw（t）小于地層破裂壓力pfrac；泡沫基液返排量小于泡沫基液注入量，即∑q′l≤∑ql；返排泡沫前沿半徑大于井筒半徑，即′≥Rw．

3 模型離散

將注入模型與返排模型在時(shí)間和空間上離散，沿井筒分為M段，時(shí)間上分為N份，將上述數(shù)學(xué)模型進(jìn)行離散.

3.1 泡沫注入模型離散

1）井筒流動(dòng)模型離散

2）泡沫滲流模型離散

3.2 泡沫返排模型離散

1）井筒流動(dòng)模型離散

2）泡沫滲流模型離散

4 模型求解

1）根據(jù)泡沫注入地層滲流模型，結(jié)合初始條件，邊界條件，得到井筒井底壓力，并作為下一時(shí)刻泡沫滲流模型初始條件.

2）由井底壓力對(duì)井筒進(jìn)行迭代求解，得到沿井筒壓力分布.

3）重復(fù)上述步驟，直到注入結(jié)束.

4）輸出所有時(shí)刻沿井筒壓力分布，泡沫質(zhì)量，泡沫密度及地層滲流狀況等.

5）將泡沫注入結(jié)束參數(shù)作為泡沫返排初始條件，根據(jù)返排泡沫地層滲流模型，得到井底壓力.根據(jù)井筒流動(dòng)模型，計(jì)算沿井筒壓力分布.結(jié)合井口回壓，試算泡沫流量.

6）試算結(jié)束得到井底壓力，作為下一時(shí)刻泡沫滲流模型初始條件.

7）重復(fù)4）到6），直到返排量等于注入量.

8）輸出返排時(shí)間，返排井底壓差，井筒泡沫流量，泡沫質(zhì)量，泡沫密度等.

5 算例

根據(jù)計(jì)算步驟，編制了計(jì)算程序.該程序可用于泡沫混排參數(shù)設(shè)計(jì)，確定井口注入壓力、返排井底壓差、返排時(shí)間、氣體和液體的流量等關(guān)鍵參數(shù)，指導(dǎo)現(xiàn)場(chǎng)的返排作業(yè).下面以一口直井為例，對(duì)泡沫流體混排的過(guò)程進(jìn)行了模擬，給出相應(yīng)的計(jì)算結(jié)果，并畫出了各參數(shù)之間相互影響的關(guān)系曲線.井地層及油井施工參數(shù)如表1所示.

5.1 注入部分計(jì)算結(jié)果分析

圖2可以看出，隨著井深的增加，井筒壓力逐漸增加.反應(yīng)了注入中泡沫沿井筒受重力及摩阻壓降的綜合影響.

圖3可以看出，泡沫質(zhì)量隨著井深增加，逐漸減小，這是由于隨井深增加壓力增加，泡沫流體中氣體流量減小，體現(xiàn)了泡沫的可壓縮性.

圖4可以看出，隨著井深增加，泡沫密度逐漸增加，這是由于壓力隨井深的增加，泡沫中氣體逐漸被壓縮，液相所占體積分?jǐn)?shù)增加.

圖5可以看出注入泡沫時(shí)井口及井底壓力逐漸增加，開(kāi)始階段增加迅速，最終逐漸趨于平緩.這是泡沫對(duì)于地層的封堵作用及泡沫破滅共同導(dǎo)致的.

5.2 返排部分計(jì)算結(jié)果分析

圖6可以看出，固定井口回壓，井底返排壓差逐漸減小.在返排開(kāi)始階段，迅速放噴造成井底壓差處于很高的狀態(tài)，隨著返排的進(jìn)行，壓力逐漸趨于平衡.

圖7可以看出，固定井口回壓，井口和井底泡沫返排量變化不大，呈略微減小趨勢(shì).泡沫流量主要包含氣體和液體兩部分流量，在返排過(guò)程中，井底壓差的變化導(dǎo)致泡沫流量的變化，而在固定井口回壓下，考慮到泡沫摩阻壓降的變化范圍不大，泡沫流量基本也穩(wěn)定在一定的范圍.

表1 地層及油井施工參數(shù)Table 1 Formation and operation parameters of a well

圖2 注入泡沫壓力沿井筒變化Fig.2 Evolution of foam injection pressure along wellbore

圖3 注入泡沫質(zhì)量沿井筒變化Fig.3 Evolution of foam quality along wellbore

圖4 注入泡沫密度沿井筒變化Fig.4 Evolution of foam density along wellbore

圖5 注入泡沫壓力隨時(shí)間變化Fig.5 Evolution of foam injection pressure

圖8和9可以看出，泡沫流量和質(zhì)量隨著井深增加，逐漸減小，這是由于隨井深增加壓力增加，泡沫流體中氣體流量減小，體現(xiàn)了泡沫的可壓縮性.

圖10可以看出，隨著井深增加，泡沫密度逐漸增加，這是由于壓力隨井深的增加，泡沫中氣體逐漸被壓縮，液相所占體積分?jǐn)?shù)增加.整個(gè)井筒內(nèi)泡沫的密度要遠(yuǎn)遠(yuǎn)小于水的密度，這充分體現(xiàn)了泡沫流體低密度的特點(diǎn).

圖6 返排井底壓差隨時(shí)間變化Fig.6 Evolution of bottomhole differential pressure during foam flowback

圖7 泡沫返排量隨時(shí)間變化Fig.7 Evolution of foam rate during foam flowback

圖8 返排泡沫流量沿井筒變化Fig.8 Evolution of flow rate during foam flowback along wellbore

圖9 返排泡沫質(zhì)量沿井筒變化Fig.9 Evolution of foam quality during foam flowback along wellbore

圖10 返排泡沫密度沿井筒變化Fig.10 Evolution of foam density during foam flowback along wellbore

6 結(jié)論

1）建立了泡沫混排數(shù)學(xué)模型，計(jì)算得到了泡沫混排過(guò)程注入及返排中，井筒內(nèi)壓力，泡沫流量，泡沫質(zhì)量，泡沫密度的分布，并得到了隨著注入時(shí)間的增加，井口及井底壓力的變化，特別在固定返排井口回壓的情況下，得到了返排井底壓差及返排時(shí)間.可以對(duì)混排過(guò)程進(jìn)行設(shè)計(jì).

2）泡沫混排中，隨著井深的增加，注入及返排中沿井筒壓力逐漸增加；沿井筒泡沫質(zhì)量及流量逐漸減小；沿井筒泡沫密度逐漸增加.隨著注入過(guò)程的進(jìn)行，井口及井底壓力增加，最終趨于平緩，施工時(shí)可以據(jù)此控制井口及井底壓力.返排時(shí)，固定井口回壓，井底壓差開(kāi)始一段時(shí)間穩(wěn)定在較高值，之后迅速降低，現(xiàn)場(chǎng)施工中，返排放噴開(kāi)始階段需要嚴(yán)格監(jiān)控.

［1］ Price B P，Gothard B.Foam?assisted lift?importance of selection and application［C］.SPE 106465，2007.

［2］ Yang J，Siddiqui S.The use of foam to improve liquid lifting from low?pressure gas wells［P］.PETSOC?99?126?P.

［3］ Chen Z，Ahmed R M，Miska S Z，et al.Experimental study on cuttings transport with foam under simulated horizontal downhole conditions［C］.SPE 99201，2006.

［4］ 楊肖曦，李松巖，林日億，等.泡沫流體攜砂能力的數(shù)值模擬［J］.中國(guó)石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2006，30（3）：89－96.

［5］ Duan Mingqin，Miska Stefan，Yu Mengjiao，et al.Experimental study and modeling of cuttings transport using foam with drillpipe rotation［C］.SPE 116300，2008.

［6］ 李兆敏，王冠華，曹小朋.泡沫混排解堵技術(shù)研究及應(yīng)用［J］.西南石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2010，32（2）：178－181.

［7］ 曹小朋.泡沫混排解堵的理論研究及應(yīng)用［D］.青島：中國(guó)石油大學(xué)（華東），2008.

［8］ 萬(wàn)里平，何保生，唐洪明，等.沖砂洗井泡沫攜砂規(guī)律數(shù)值模擬研究［J］.西南石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2013，35 （4）：101－106.

［9］ 袁照永.水平井環(huán)空泡沫攜巖流動(dòng)規(guī)律研究［D］.青島：中國(guó)石油大學(xué)（華東），2009.

［10］ 王玉斗，李茂輝，溫科揚(yáng)，等.泡沫滲流激勵(lì)及滲流模型研究［J］.石油鉆探技術(shù)，2010，38（4）：104－107.

［11］ Sim joo M，Iran，Zitha P L J.Effects of oil on foam generation and propagation in porousmedia［C］.SPE 165271，2013.

［12］ 張小寧，李根生，黃中偉，等.泡沫鉆井液在井筒中的流動(dòng)與傳熱［J］.石油學(xué)報(bào)，2010，31（1）：134－138.

［13］ 李松巖，李兆敏，林日億.泡沫舉升排酸過(guò)程井筒壓力溫度數(shù)學(xué)模型研究［J］.西南石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2009，31（2）：59－63.

［14］ 李松巖.氮?dú)馀菽至魉峄に嚰夹g(shù)研究［D］.青島：中國(guó)石油大學(xué)（華東），2009.

［15］ 周生田，李兆敏，王飛.水平井泡沫酸化分流數(shù)值模擬［J］.石油學(xué)報(bào)，2012，33（3）：448－452.

［16］ 李兆敏，孫茂盛，林日億，等.泡沫封堵及選擇性分流實(shí)驗(yàn)研究［J］.石油學(xué)報(bào)，2007，28（4）：115－118.

Received date： 2014－01－14；Revised date： 2014－06－10

Article ID：1001?246X（2015）01?0065?10

Influence of Fem tosecond Laser Chirp on Optical Limiting

and Dynam ical Two?photon Absorption of 4，4'?bis（di?n?butylam ino）stilbene Com pounds

Abstract： Optical limiting（OL）and dynamical two?photon absorption（TPA）of 4，4'?bis（di?n?butylamino）stilbene（BDBAS）molecules in chirped femtosecond laser pulses are studied by solving Maxwell?Bloch equationswith an iterative predictor?corrector finite?difference time?domain（FDTD）method.It shows that both sign and magnitude of chirp rate influence greatly spectrum evolution and OL behavior.Spectra exhibit obvious carrier frequency shifts depending on sign of chirp rate，blue shift for positive chirp rate and red shift for negative chirp rate.As absolute chirp rate increases，shift becomesmore obvious，OL window gets narrower and saturation of output intensity becomes greater.Interestingly，self?induced transparency（SIT）appears as a negative chirp rate reduces to a certain value（－0.025 fs－2）.Dynamical TPA cross section is reduced as chirp effect is considered.It provides a method for controlling nonlinear optical absorptions.

Key words： two?photon absorption；optical limiting；chirped pulse；organicmolecule

0 Introduction

With remarkable progress in ultrashort and ultraintense laser technology，it is possible to generate precisely defined laser pulses［1］which brings a vast variety of investigations in the light?matter interaction［2－5］and leads to great potential of applications［6］.Among all of the applications，optical limiter，whosemajormechanism is dynamical two?photon absorption（TPA），hasmotivated an extensive research in order to protect delicate sensors，especially human eyes.On the basis of quantum approaches at ab initio level，the TPA cross section of a molecule was theoretically calculated［7－8］.This TPA cross section is normally noted as static since it is only referred to the molecule itself.When the TPA cross section of a molecular system is measured，there exists interaction between the molecule and laser.Thus，one should solve the combined Maxwell?Bloch equations［6，9－10］to simulate correctly experimental results and explore nonlinear optical processes. As a result，the TPA cross section performs as a dynamical value because it is related to laser parameters.

Considering that a laser pulse usually has a large frequency chirp that may change characteristics of the pulse significantly［11－12］，one expect that the laser chirp would have influence on nonlinear optical properties of molecules.In the last decades，chirped laser pulses have been widely applied to generate high?order harmonic（HHG）［13］，transfer population［14］，and transform multiphoton between Rydberg states［15］，etc.Desaix et al［16］studied chirped soliton pulses propagating in optical nonlinear Kerr media，indicating obvious influences on properties ofasymptotically emerging solitons by the initial chirp.Moreover，Centinietal［17］measured signal and group velocities of chirped pulses propagating through a GaAs cavity，which showing substantial modification of the group velocity of pulse by the chirp.Song et al［18］investigated theoretically coherent control of spectra with chirped femtosecond laser pulses propagating in a two?level?atom medium.They found that the sign and magnitude of the chirp rate have obvious effects on the spectral feature.

To the best of our knowledge，there are few works to study influence of chirp on OL properties and dynamical TPA cross sections of organic compounds.In this paper，taking the BDBAS compound as the medium reported by Ehrlich［19］，we investigate spectra of chirped femtosecond laser pulses and OL behaviors of the compounds by solving Maxwell?Bloch equations using an iterative predictor?corrector FDTDmethod［20］.Chirp effecton TPA cross section is analyzed further.

1 Theoreticalmethods

1.1 M axwell?Bloch equations of a three?level system

Based on a semiclassical theory for interaction between laser pulses and molecules，electromagnetic radiation is described classically by Maxwell equations，and themolecular system is treated quantum mechanically with Bloch equations.The density matrix equation with relaxation effect is written as

here，Γnmgives the rate per atom atwhich population decays from level m to level n，andγnmgives relaxation rate of the density matrix elementρnm．is Hamiltonian of the system，which can be expressed as the sum of Hamiltonian of the field?free moleculeand the interaction Hamiltonian H＾′．Within the dipole approximation，it can be written as

Settingρ01＝（u0＋i v0）/2，ρ12＝（u1＋i v1）/2，ρ02＝（u2＋i v2）/2，we get Bloch equations for the three?level system ?u

where dmn＝dnm（m，n＝0，1，2）are the permanent electric dipolemoments（m＝n）or the transition electric dipolemoments（m≠n）of the molecule，and?ωmnis the excitation energy between the states m and n．

We assume that the incident electromagnetic field polarizes along the x axis and propagates along the z axis to an input interface at z＝0．Then the Maxwell equations take the form

whereμ0andε0are the permeability and permittivity of free space，respectively.

The Maxwell and the Bloch equations are coupled to each other by themacroscopic polarization Px，

where N is themolecular density.

1.2 Dynam ical TPA cross section

In the absence of significant recombination，diffusion and thermal，the differentialequation that describes propagation of pulses in the presence of linear and TPA can be written as［21］

where I is the field intensity，αdenotes the linear absorption coefficientandβis the TPA coefficient. The analytical solution of Eq.（6）is

As one knows，the TPA coefficient depends on the input intensity I0［22］.If I0is smaller than the saturation absorption intensity，the TPA coefficientβcan be assumed as a linear function of I0［23］，

whereβ0is the steady?state TPA coefficient.Input Eq.（8）into（7），one can see that the inverse transmission 1/T is a quadratic function of the input field intensity I0，

αandβ0can thus be estimated from Eq.（9）.Themolecular TPA cross sectionσtpis related toβ0by

where hνis the input photon energy.

2 Results and discussions

As shown in Fig.1，BDBAS is a typical one?dimensional symmetrical conjugated organic molecule.The time?dependent DFT/B3LYP method implemented in DALTON package［24］is employed to obtain the excitation energy Enand the dipolemoment dmnof themolecule.The 6?31G basis set is chosen for all calculations.It is shown in Table 1 that thismolecule has only one charge transfer state in low energy region，which is the first excited state S1withδop＝1.630.The transition dipolemoment between S1and the fourth excited state S2ismuch larger（d14＝4.154×10－29C·m）than the others.ab initio calculations show that the fourth excited state S2has the maximum TPA cross section in the low energy region.The excitation energies of S1and S2are 3.408 eV and 4.206 eV，respectively.The permanent dipole moments are approximately equal to zero because of symmetry of themolecule.When the interaction between themolecule and ultrashortpulse in the low energy region is considered，a cascade three?level system including the ground state S0，the first excited state S1and the fourth excited state S2ismodeled as shown in Fig.1.

Fig.1 Molecular structure of BDBAS（left）and transitions（right）

Table 1 Excitation energy En，oscillator strengthδopand dipolemoment dmnof the first five excited states

We choose the initial chirped electric field with a hyperbolic?secant functional for Ex（z＝0，t）＝F0sech［1．76（t－z0/c）/ ］cos［ω0t＋χt2/2］，where F0is the peak amplitude of the inputelectric field， is the fullwidth at halfmaximum（FWHM）of the pulse intensity envelope，which is set as a few femtoseconds，andχis the chirp rate.The choice of z0ensures that the pulse penetrates negligibly into themedium at t＝0．Themolecule is assumed in its ground state before the light is turned on，i.e.，ρ00（t＝0）＝1，ρ11（t＝0）＝ρ22（t＝0）＝0．Themolecular density is taken as N＝7.0×1025m－3.The decay rates of the densitymatrix elementsγnmcan be chosen as 1.0× 1013s－1，while the decay rates of excited statesΓ01，Γ12andΓ02are assumed equal to 1.0×109s－1，1.0×1012s－1and 0［25］，respectively.The time and space incrementsΔt andΔz are chosen to ensure cΔt≤Δz［26］.

2.1 Two?photon resonant propagation of fem tosecond pulse

Frequency of the incident ultrashort pulseω0is taken as half of the frequency between S0and S2states，namely，ω0＝ω20/2．Taking FWHM of the incident pulse as 5 fs，we simulated temporal evolution of pulses with different input peak amplitudes F0＝4.0×109V·m－1，6.0×109V·m－1and 8.0×109V·m－1for chirp ratesχ＝0，0.002 5 fs－2and－0.002 5 fs－2，respectively（Fig. 2）.As one can see in Fig.2，the pulses increase area by splitting，not by pulse broadening as in the long pulse cases，which is in agreementwith a previouswork［27］.Thismeans that the traditional area theorem based on the slowly?varying envelope approximation and rotatingwave approximation［28］is invalid in describing an ultrashort pulse propagation accurately.As the input peak amplitudes increases，the splitbecomes obvious.Compared with the chirp free case，the splits of chirped pulses are weaker，which is consistentwith the two?level case［18］.When a chirp rate exists，themain pulse has little variation，and the subpulse intensity becomes weaker.The time?delay between the main pulse and subpulse turns longer compared with the case without a chirp rate.

Fig.2 Evolution of electric fields with different input peak amplitudes（a）χ＝0，（b）χ＝0.002 5 fs－2，and（c）χ＝－0.002 5 fs－2with pulse width 5 fs at propagation distance of 7.0μm

Corresponding spectra of Fig.2 are shown in Fig.3.One can see that except for the base frequencyω0，lower and higher frequencies appear mainly due to self?phase modulation（SPM）during the propagation of pulse，such asω21，ω10，ω10＋2ω21，3ω0and 5ω0.The odd harmonic components（3ω0and 5ω0）appear for the symmetry compounds with a strong TPA process，while the even harmonic components are restrained.With the enhancement of peak amplitudes，the spectra display more components.When chirp rate is considered，the spectra of chirped pulses exhibit obvious shift of carrier frequency during propagation in medium depending on the sign of chirp rate，namely，blue shifts for positive chirp rates and red shifts for negative cases，which is consistentwith previous results［18］.

Fig.3 Corresponding spectra in Fig.2（a）χ＝0，（b）χ＝0.002 5 fs－2，and（c）χ＝－0.002 5 fs－2.

In order to further elucidate the influence of chirp on ultrashort pulse propagation，we consider different chirp rates.Fig.4（left）shows a chirped pulse of F0＝4.0×109V·m－1withχ＝±0.005 fs－2andχ＝±0.025 fs－2，propagating through BDBAS at respective distances of0 and 14.0μm.One can see that as absolute value of the chirp rate increases，the intensity of main pulse increases obviously，while the intensity of subpulse components decreases.Corresponding spectra are shown in Fig.4（right）.One can see that the spectra shift ismore evident for pulse with a larger absolute value of chirp rate.Moreover，higher and lower frequency components（ω21andω10）disappear for a chirped pulse withχ＝±0.025 fs－2.The center frequency deviates further from TPA resonant frequencyω20/2 as the absolute value of chirp rate increases，leading to weaker TPA from S0to S2and thus smaller population of S2.Specially，when the chirp rate takes－0.025 fs－2，themedium becomes approximately transparent to pulse.

Fig.4 Evolution of electric fields（left）at distances of 0（solid line）and 14.0μm（dash line）and corresponding spectra（right）at 14.0μm with different chirp rates（F0＝4.0×109V·m－1， ＝5 fs）

In order to further observe this phenomenon，F(xiàn)ig.5 shows a chirped pulse of F0＝4.0×109V· m－1withχ＝－0.025 fs－2propagating through BDBAS at respective distances of 0，3.5μm，7.0μm，10.5μm，14.0μm and 17.5μm（Fig.5（a））along with the corresponding spectra（Fig.5 （b））.One can see that the shape of pulse is quite stable and no additional frequency appears.This can be explained by the resonant frequency shiftandmuch weak absorption during pulse propagating process.

Fig.5 （a）Chirped pulse ofχ＝－0.025 fs－2propagating through BDBAS at respective distances of 0μm （solid line），3.5μm（dash line），7.0μm（dot line），10.5μm（dash dot line），14.0μm（dash dot dot line）and 17.5μm（short dash line），and（b）corresponding spectra（F0＝4.0×109V·m－1， ＝5 fs）

2.2 OL and dynam ical TPA cross section

Fig.6 Output fluence Soutat z＝7.0μm versus input fluence Sinof a 5 fs pulses

To illustrate influence of chirp rate on OL ability，we present output fluence Soutat z＝7.0μm as a function of input fluence Sinof 5 fs pulses with different chirp ratesχ＝0， ±0.002 5 fs－2， ± 0.005 fs－2and－0.025 fs－2in Fig.6.As shown in Fig.6，one can see that the chirp rate has an obvious influence on OL in both OL window and higher output saturation intensity.When the absolute value of chirp rate increases，the medium shows a narrow OL window and higher output saturation intensity，showing weaker OL ability.Specially，when the chirp rate reaches－0.025 fs－2，the medium turns transparent，meaning breakdown of OL behavior.This phenomenon is due to the resonant frequency shift caused by chirp rate，which results in much weak absorption from ground state S0to the maximum TPA state S2.

Dependence of dynamical TPA cross sectionσtpon pulse width in femtosecond time domain with different chirp rates ata distance of7.0μm is shown in Fig.7.One can see that the TPA cross section increases as thewidth of the input pulse is broadened.The TPA cross sections is obtained as 1 145.3 GM（1 GM＝1×10－50cm4s/photon）for chirp free pulse width of 5 fs at a distance of7.0μm，which is in the same order of magnitude with the static TPA cross sectionσtp＝1 430 GM.The calculated value is smaller than the experiment result of σtp＝9 300 GM［19］.The discrepancy mainly results from the experimental condition with pulse duration of 5 ns andλ＝600 nm under which a two?step TPA process takes place.Furthermore，the TPA cross section is reduced by the chirp，and theχ＝0.002 5 fs－2pulse shows a larger TPA cross section than theχ＝－0.002 5 fs－2pulse.Therefore，both sign and magnitude of pulse chirp rates should be considered when calculating TPA cross sections of different pulses.

Fig.7 Dependence of TPA cross sectionσtpon pulse widths and chirp rates at a distance of 7.0μm

3 Conclusions

In conclusion，we investigate influence of chirp rates on OL and dynamical TPA in a cascade three?level medium BDBAS theoretically by solving Maxwell?Bloch equations using an iterative predictor?corrector FDTDmethod.It is demonstrated that both sign and magnitude of the chirp rate influence spectra and nonlinear properties of the medium crucially.The OL ability breaks down as magnitude of the negative chirp rate is larger than 0.025 fs－2.The OL behavior ismore obvious and the dynamical TPA cross section ismuch larger for a positive chirp rate laser compared with negative one.It shows thatboth sign andmagnitude should be taken into consideration when nonlinear optical properties of medium are measured and simulated.Our investigation predicts a way to control nonlinear optical absorption by changing chirp rate.

［1］ Kong L J，Zhao L M，Lefrancois S，et al.Generation ofmegawatt peak power picosecond pulses from a divided?pulse fiber amplifier［J］.Opt Lett，2012，37（2）：253－255.

［2］ Liu Z B，F(xiàn)an X J.Phase correlation in propagation effect of a doppler broadening quasi?Λtype four?level system［J］.Chinese JComput Phys，2012，29（6）：881－890.

［3］ Tan X，Wang L，Wang Z D，et al.Population of dense V?type three?level medium in few?cycle laser pulses［J］.Chinese JComput Phys，2009，26（5）：773－780.

［4］ Zhou Y，Miao Q，Wang C K.Effect of Lorentz local field correction on propagation of ultrashort laser pulse in one?dimensional para?nitroaniline（PNA）molecules［J］.Chinese Phys B，2011，20（4）：044205（1－5）.

［5］ Wang L，Wang Z D，Lian B，et al.Spatial distribution of few?cycle laser pulses propagating in a ladder?type three?levelmedium［J］.Chinese JComput Phys，2011，28（1）：81－86.

［6］ Liu JC，Wang C K，Gel'mukhanov1 F.Optical limiting of short laser pulses［J］.Phys Rev A，2007，76 （5）：053804（1－6）.

［7］ Sa?ek P，Vahtras O，Helgaker T，et al.Density?functional theory of linear and nonlinear time?dependent molecular properties［J］.JChem Phys，2002，117（21）：9630－9645.

［8］ Ding H J，Sun J，Zhang Y J，et al.Solvent effects on one?and two?photon optical properties of a novel2，2'：6'，2”?terpyridine?based chromophore［J］.JAt Mol Sci，2013，4（4）：349－356.

［9］ Wang C K，Zhao P，Miao Q，et al.Optical limiting and dynamical two?photon absorption of organic compounds for a nanosecond pulse［J］.JPhys B：At Mol Opt Phys，2010，43（10）：105601（1－5）.

［10］ Zhou Y，Miao Q，Sun Y P，et al.Solvent effect on dynamical TPA and optical limiting of BDMAS molecularmedia for nanosecond and femtosecond laser pulses［J］.JPhys B：AtMol Opt Phys，2011，44 （3）：035103（1－7）.

［11］ Zheng H J，Liu S L，Xu JP.Research on nonlinear propagation of chirped exponential pulses［J］.Stud On Opt Commun，2007，2（140）：11－41.

［12］ Zheng H J，Liu S L，Lin Xin，et al.Effect of initial frequency chirp on the linear propagation characteristics of the hyperbolic secant optical pulse［J］.Acta Phys Sin，2007，56（4）：2286－2292.

［13］ Lee D G，Kim JH，Hong K H，et al.Coherent control of high?order harmonicswith chirped femtosecond laser pulses［J］.Phys Rev Lett，2001，87（24）：243902（1－4）.

［14］ Netz R，Nazarkin A，Sauerbrey R.Observation of selectivity of coherent population transfer induced by optical interference［J］.Phys Rev Lett，2003，90（6）：063001（1－4）.

［15］ Conover CW S，Doogue M C，Struwe F J.Chirped?pulse multiphoton transitions between Rydberg states ［J］.Phys Rev A，2002，65（3）：033414（1－8）.

［16］ Desaix M，Helczynski L，Anderson D，et al.Propagation properties of chirped soliton pulses in optical nonlinear Kerrmedia［J］.Phys Rev E，2002，65（5）：056602（1－7）.

［17］ Centini M，Bloemer M，Myneni K，etal.Signal velocity and group velocity for an optical pulse propagating through a GaAs cavity［J］.Phys Rev E，2003，68（1）：016602（1－4）.

［18］ Song X H，Gong SQ，YangW F，et al.Coherent control of spectra effectswith chirped femtosecond laser pulse［J］.Opt Commun，2004，236（1）：151－157.

［19］ Ehrlich JE，Wu X L，Lee I?Y S，et al.Two?photon absorption and broadband optical limiting with bis?donor stilbenes［J］.Opt Lett，1997，22（24）：1843－1845.

［20］ Ziolkowski RW，Arnold JM，Gobny D M.Ultrafast pulse interactionswith two?level atoms［J］.Phys Rev A，1995，52（4）：3082－3094.

［21］ Bhawalkar JD，He G S，Prasad P N.Nonlinearmultiphoton processes in organic and polymeric materials ［J］.Rep Prog Phys，1996，59（9）：1041－1070.

［22］ He G S，Zheng Q D，Baev A，et al.Saturation ofmultiphoton absorption upon strong and ultrafast infrared laser excitation［J］.JAppl Phys，2007，101（8）：083108（1－6）.

［23］ Sun Y P，Liu J C，Wang C K.Effect of time?dependent ionization on dynamical two?photon absorption cross sections ofmolecularmedia［J］.Acta Opt Sin，2009，29（6）：1621－1625.

［24］ DALTON.A molecular electronic structure program［CP/OL］.Dalton 2013.http：//dalton?program.org/. ［25］ Gavrilyuk S，Polyutov S，Jha PC，et al.Many?photon dynamics of photobleaching［J］.JPhys Chem A，2007，111（47）：11961－11975.

［26］ Taflove A，Brodwin M E.Numerical solution of steady?state electromagnetic scattering problems using the time?dependent Maxwell's equations［J］.Micro Theo Tech，1975，23（8）：623－630.

［27］ Xiao J，Wang Z Y，Xu Z Z.Area evolution of a few?cycle pulse laser in a two?level?atom medium［J］. Phys Rev A，2002，65（3）：031402（1－4）.

［28］ Mccall S L，Hahn E L.Self?induced transparency［J］.Phys Rev，1969，183（2）：457－485.

M athematical M odel and Numerical Simulation of Foam Plug Removal

WANG Fei1，LIZhaomin1，LISongyan1，DU Qingjun2
（1.College ofPetroleum Engineering，China University ofPetroleum，Qingdao 266580，China；
2.School ofGeosciences，China University of Petroleum，Qingdao 266580，China）

With amodel of foam wellbore flow and a model of foam seepage flow，a mathematicalmodel of foam plug removal is given.Themodel is solved with numerical method.Distributions of foam pressure，foam quality，foam density along wellbore and wellhead and bottom pressure are discussed.Furthermore，variation of bottom hole differential pressure is given as wellhead back pressure is fixed.It shows that as foam pressure and density increases，foam quality decreases with increase of well depth and bottom hole differential pressure declines.

foam plug removal；seepage flow；wellbore flow；mathematicalmodel；bottomhole differential pressure

ZHANG Yujin，ZHANG Qiuyue，SONG Yuzhi，WANG Chuankui

（College ofPhysics and Electronics，Shandong Normal University，Jinan 250014，China）

O437；O561.1 Document code：A

TE319

A

2014－01－14；

2014－06－10

國(guó)家高科技研究發(fā)展計(jì)劃（863計(jì)劃）（2013AA064801）和國(guó)家自然科學(xué)基金青年基金（11102236）資助項(xiàng)目

王飛（1988－），男，博士生，從事泡沫流體模擬計(jì)算及在油田中的應(yīng)用等研究，E?mail：goodboygreatman＠163.com

Received date：2014－01－08；Revised date：2014－03－30

Foundation item s：Supported by 973 Program（2011CB808100）

Biography：Zhang Yujin（1989－），female，PhD，major in molecular photonics，E?mail：zhangyujin312＠163.com