朱奎鑫，蓋 孟，賴舜男

大時(shí)間步長下體積守恒的實(shí)時(shí)水體模擬

朱奎鑫1，蓋 孟2,3，賴舜男2

(1. 北京大學(xué)深圳研究生院，廣東 深圳 518055；2.北京大學(xué)信息科學(xué)技術(shù)學(xué)院，北京 100871；3. 北京大學(xué)北京市虛擬仿真與可視化工程中心，北京 100871)

基于物理的流體仿真可以真實(shí)地捕捉水的運(yùn)動(dòng)，一種常用的實(shí)時(shí)模擬方法是基于二維淺水方程進(jìn)行水體仿真。首先介紹了二維淺水方程，并提出一種新的求解方法，通過一種體積守恒的隱式半拉格朗日方法進(jìn)行計(jì)算，在保持穩(wěn)定的同時(shí)允許較大的時(shí)間步長，然后額外增加修正步以保證計(jì)算過程中體積守恒。此外提出了一種基于動(dòng)量守恒的流固耦合方法，可以實(shí)時(shí)地模擬出較為真實(shí)的水體，并且保證了水體體積始終守恒，效果良好。

淺水方程；半拉格朗日；體積守恒；大時(shí)間步長；流固耦合

流體運(yùn)動(dòng)是一個(gè)復(fù)雜的現(xiàn)象，很難實(shí)時(shí)在計(jì)算機(jī)上進(jìn)行有效模擬，折中的方法是通過模擬水面來近似模擬水體，只需要使用二維區(qū)域上的高度場即可進(jìn)行實(shí)時(shí)模擬。常用的方法有2種，一種是基于譜的模擬[1]，此方法非常穩(wěn)定且細(xì)節(jié)豐富，但缺點(diǎn)是難以與移動(dòng)的障礙物相互作用；另一種是使用基于物理的淺水方程進(jìn)行模擬。本研究以二維淺水方程組為基礎(chǔ)，提出了一種高效、穩(wěn)定、體積守恒的流體仿真算法。通過半拉格朗日方法的使用，實(shí)現(xiàn)了較大的時(shí)間步長下的穩(wěn)定模擬，這種基于物理的流體模型能夠產(chǎn)生真實(shí)的波運(yùn)動(dòng)，對(duì)流體做到實(shí)時(shí)合理的仿真。

本文首先對(duì)淺水方程進(jìn)行了描述，以此為基礎(chǔ)建立了流體仿真模型。淺水方程是由納維-斯托克斯方程推導(dǎo)而來，雖然無法模擬三維流體或帶有黏度的流體，但是由于在高度項(xiàng)上進(jìn)行了簡化，計(jì)算量大大降低，使得求解效率提高，可以實(shí)現(xiàn)對(duì)流體的實(shí)時(shí)模擬。

然后，詳述了一種質(zhì)量守恒的半拉格朗日方法。人們知道，常規(guī)的半拉格朗日方法無法做到質(zhì)量守恒，無法滿足物理仿真與實(shí)際相符，通過在常規(guī)方法之后進(jìn)行的擴(kuò)散操作的改進(jìn)，保證了平流前后質(zhì)量守恒，對(duì)于流體仿真具有重大意義。

之后，對(duì)基于二維淺水方程的流體仿真進(jìn)行實(shí)現(xiàn)。在該模型中，采用改進(jìn)的隱式半拉格朗日積分法求解，再使用修正步彌補(bǔ)求解時(shí)的體積損耗，在保證穩(wěn)定性的同時(shí)，允許較大的時(shí)間步長。

在此模型上實(shí)現(xiàn)了流固耦合的技術(shù)，借助有向距離場來判斷交互區(qū)域，然后基于動(dòng)量守恒進(jìn)行單向交互，取得了很好的效果。

最后，通過演示該模型展示流體模擬的效果，并與其他模型進(jìn)行了比較，證明了本文方法體積守恒、算法穩(wěn)定、計(jì)算效率高。

1 淺水方程

淺水方程是描述淺水流動(dòng)的數(shù)學(xué)模型，常用于實(shí)時(shí)流體運(yùn)動(dòng)的仿真計(jì)算。在流體動(dòng)力學(xué)中，通常水平長度比垂直長度大得多，在此條件下，可以得到沿水深方向的靜水壓力分布，并根據(jù)納維-斯托克斯方程導(dǎo)出淺水方程。有時(shí)也會(huì)引入水深方向的平均化，導(dǎo)出淺水方程的另一種形式。本文使用的淺水方程式為

其中，為水平面高度；為水的深度，即水平面高度減去地形高度；為速度；為重力加速度；為時(shí)間。

淺水方程沒有黏性項(xiàng)，水的深度方向速度可忽略，所以模擬時(shí)僅限于二維無黏流體，對(duì)于高黏度流體，飛濺或破裂的流體無法模擬。然而，其可充分地模擬一大類流體的運(yùn)動(dòng)，且計(jì)算效率極高，是進(jìn)行實(shí)時(shí)流體模擬的最好選擇。

求解淺水方程有多種方法，中外研究者大都采用有限體積法以保證體積守恒[2-3]，不過該方法在較大時(shí)間步長下表現(xiàn)不穩(wěn)定；有些隱式差分方法，可以保證無條件穩(wěn)定[4]，卻不能保證體積守恒，本文方法希望可以同時(shí)具有這2個(gè)優(yōu)點(diǎn)。

2 半拉格朗日方法

半拉格朗日方法[5]可以看作是歐拉方法和拉格朗日方法的混合體。歐拉方法保持網(wǎng)格的規(guī)則性，但為了保持穩(wěn)定，需要較小的時(shí)間步長。拉格朗日方案不受穩(wěn)定性要求的限制，允許更大的時(shí)間步長。半拉格朗日方法試圖將2種方案的優(yōu)點(diǎn)，即歐拉方法的規(guī)則性和拉格朗日方法的穩(wěn)定性結(jié)合起來，在每個(gè)時(shí)間點(diǎn)，進(jìn)行網(wǎng)格點(diǎn)值計(jì)算時(shí)，沿著粒子軌跡進(jìn)行求解。

以式(1)為例，等式左邊展開為

公式右邊的第1項(xiàng)稱為局部項(xiàng)，第2項(xiàng)為平流項(xiàng)。其實(shí)，半拉格朗日法就是用拉格朗日的思想求解平流項(xiàng)，然后在固定的網(wǎng)格上求解局部項(xiàng)。

類似的，有

其中，為當(dāng)前時(shí)刻網(wǎng)格的速度。

圖1 半拉格朗日方法

圖2中的()和偽代碼中的()代表的是α，為每個(gè)網(wǎng)格的流入量權(quán)重，而()代表的是網(wǎng)格存儲(chǔ)的物理量。在循環(huán)過程中，α逐漸向1靠攏。

上述擴(kuò)散權(quán)重的方法，也存在一個(gè)問題，即迭代次數(shù)過多，對(duì)實(shí)時(shí)仿真造成影響。一種加快該過程的方法是繼承上一次的誤差，在每幀權(quán)重?cái)U(kuò)散之前，加上上一幀最后計(jì)算的結(jié)果與1的差，然后再進(jìn)行權(quán)重?cái)U(kuò)散。在實(shí)際工作中，這樣可以減少迭代次數(shù)，提高仿真效率。

圖2 擴(kuò)散權(quán)重方法

3 方程求解

進(jìn)行方程求解模擬流體運(yùn)動(dòng)的流程如下：

(1) 使用本文改進(jìn)的半拉格朗日方法進(jìn)行平流。

(2) 將方程離散化，隱式求解。

(3) 考慮求解中損失的體積，在修正時(shí)進(jìn)行體積補(bǔ)償。

3.1 半拉格朗日平流

半拉格朗日平流，使用本文改進(jìn)的半拉格朗日方法，滿足任意β和α均等于1的條件，這樣既保證了體積守恒，也實(shí)現(xiàn)了較合理地仿真效果。使用了繼承誤差的方法后，在仿真過程中，迭代次數(shù)不超過10次即可收斂，可既達(dá)到目的，也滿足了的實(shí)時(shí)仿真的要求。

3.2 隱式求解

在進(jìn)行平流之后，部分研究者采用半隱式的方法[8]，本文采用文獻(xiàn)[9]隱式求解的方法，1階精度，但對(duì)于流體模擬已經(jīng)足夠，式(1)和式(2)可以寫成如下形式

式(7)~(9)中，為速度方向的分量；為速度方向的分量。由于采用隱式求解，所以等號(hào)右邊均使用了下一時(shí)刻的狀態(tài)，式(7)等號(hào)右側(cè)的h，本應(yīng)為h+1，但這將變成非線性方程組求解問題，求解困難，因此本文使用h來近似h+1。

現(xiàn)在開始討論網(wǎng)格如何離散化。在整個(gè)模擬過程中，必須在空間的各個(gè)點(diǎn)存儲(chǔ)許多不同的量(速度、高度等)。很明顯，需要采用一種合理的布局。最簡單的方法是，將所有信息存儲(chǔ)在同一個(gè)網(wǎng)格中心上，但這并不是最佳的方法。很早以前，MAC網(wǎng)格就被提出[10]，目前在流體模擬中被廣泛使用。該網(wǎng)格也被稱為交錯(cuò)網(wǎng)格，因其在不同的位置存儲(chǔ)不同的物理量，如在網(wǎng)格中心存儲(chǔ)高度，網(wǎng)格邊上存儲(chǔ)速度。在二維情況，MAC網(wǎng)格中的單元格結(jié)構(gòu)如圖3所示。

根據(jù)這個(gè)網(wǎng)格結(jié)構(gòu)，可將高度、速度相關(guān)的公式進(jìn)行空間離散化，形式如下

將速度公式帶入到高度公式中即可得到一個(gè)線性方程組，使用雅克比迭代求解該大型稀疏線性方程組，時(shí)間步長較小時(shí)，迭代10次以內(nèi)就可以得到滿意的結(jié)果。

3.3 體積修正

在實(shí)現(xiàn)的流體模擬中，每一步都要保證體積守恒。平流時(shí)采用了改進(jìn)的半拉格朗日方法，保證了平流前后體積守恒，但是在進(jìn)行隱式求解時(shí)，無法保證體積守恒，為此，需要找到求解過程中的體積損失。

根據(jù)高度公式，不難發(fā)現(xiàn)，如果要保證最后體積守恒，線性方程組的所有方程等號(hào)右側(cè)相加應(yīng)該等于0，隱式求解時(shí)的方程組并未滿足該條件，所以體積不守恒?，F(xiàn)在要對(duì)體積進(jìn)行修正，在更新高度時(shí)，等號(hào)右側(cè)統(tǒng)一地加入一個(gè)相同形式的量，使得所有網(wǎng)格等號(hào)右側(cè)相加等于0。

通過上述描述已可以實(shí)現(xiàn)在開闊區(qū)域的水面模擬，如圖4所示。

圖4 水面上的水波

3.4 邊界處理

干濕邊界處理在淺水方程中也很重要，有很多方法得到了廣泛應(yīng)用[11]，將邊界進(jìn)行標(biāo)記，然后特殊處理。對(duì)于本文的模擬來說，在水流開始運(yùn)動(dòng)之前，有水的網(wǎng)格水平面應(yīng)該高于附近沒有水的網(wǎng)格地形。否則，兩者之間的面就被標(biāo)為障礙。對(duì)于在MAC網(wǎng)格中標(biāo)記為障礙的面，速度在每個(gè)時(shí)間步結(jié)束時(shí)設(shè)置為0，而且在平流步中不更新。滿足以下條件的面視為障礙，即

其中，為一個(gè)大于0的小常數(shù)，均定義為0.000 1。類似的情況也適用于交錯(cuò)網(wǎng)格的其他面。按照本文的處理方法，圖4所示的水體碰到邊界后如圖5所示。

圖5 水波碰到邊界

3.5 流固耦合

流固耦合在水體模擬中是一個(gè)重要的話題，為了使剛體與水體之間能有效地相互作用，本文提出了一種模擬流固相互作用的單向耦合方法，即水體的運(yùn)動(dòng)只能受到剛體的影響，而不對(duì)剛體的運(yùn)動(dòng)產(chǎn)生影響。

流固耦合有2個(gè)關(guān)鍵技術(shù)，一是如何確定交互區(qū)域，二是確定交互區(qū)域之后的相互作用。

使用有向距離場來確定交互區(qū)域。有向距離場(signed distance field, SDF)，“有向”、“距離”、“場”這3個(gè)詞非常精確的描述了SDF究竟是什么。SDF是到(多邊形模型)物體表面最近距離的采樣網(wǎng)格。距離的正負(fù)號(hào)表示在物體的外部或內(nèi)部，距離值為正表示點(diǎn)在物體外部，反之在內(nèi)部，若距離值為零表示點(diǎn)在物體表面?？蓪⒔换サ墓腆w模型文件轉(zhuǎn)化為有向距離場。對(duì)于淺水方程的每個(gè)網(wǎng)格，可對(duì)高度進(jìn)行離散，依次就此網(wǎng)格的當(dāng)前高度判斷距離值是否為負(fù)，并判斷是否在交互物體內(nèi)部。如果在交互物體內(nèi)部，那么此網(wǎng)格的當(dāng)前高度就是交互區(qū)域，這樣，在二維的水面上得到了一個(gè)三維的交互區(qū)域。

確定交互區(qū)域之后，根據(jù)動(dòng)量守恒進(jìn)行交互。對(duì)于當(dāng)前網(wǎng)格，如果處于交互區(qū)域，認(rèn)為交互區(qū)域的速度等于固體的速度。然而，由于交互區(qū)域是三維的，可能某些高度處于交互區(qū)域，或不屬于交互區(qū)域，但是對(duì)于當(dāng)前網(wǎng)格，只有一個(gè)速度，即

其中，1為固體速度；2為交互之前的水體速度；3為交互之后的水體速度；1為當(dāng)前網(wǎng)格根據(jù)深度離散化后，處于交互區(qū)域的高度和；2為不處于交互區(qū)域的高度和；3為水深，1+2=3。

交互算法流程如下：

步驟1.對(duì)于交互的固體模型轉(zhuǎn)化為有向距離場。

步驟2. 對(duì)于每個(gè)時(shí)間步長：

(1) 確定有向距離場在整個(gè)計(jì)算網(wǎng)格上的大致投影區(qū)域。

(2) 對(duì)于投影區(qū)域的每個(gè)網(wǎng)格：①將網(wǎng)格上的高度離散化，從上到下遍歷每一個(gè)離散坐標(biāo)，判斷坐標(biāo)點(diǎn)是否在物體內(nèi)部；②對(duì)每個(gè)網(wǎng)格，按照動(dòng)量守恒，更新當(dāng)前網(wǎng)格速度。

(3) 繼續(xù)進(jìn)行淺水方程計(jì)算。

按照此方法，水面上劃過的木塊尾跡，如圖6所示。

圖6 水面上的劃過的木塊

4 結(jié)果和分析

4.1 實(shí)驗(yàn)結(jié)果

本文模型可以模擬許多的流體運(yùn)動(dòng)。為了滿足實(shí)時(shí)要求，還借助了GUDA進(jìn)行GPU并行加速。所有的流體模擬都是在CPU為inter(R) Core(TM) i5-4570 3.20 GHz，顯卡為NVIDIAGeForce GTX970，操作系統(tǒng)為Windows10，使用CUDA9.1進(jìn)行并行加速計(jì)算。

圖7使用固定的邊界，對(duì)帶有障礙物的矩形水池中的水波進(jìn)行了實(shí)時(shí)模擬。此實(shí)驗(yàn)中模擬的網(wǎng)格大小為256×256，每幀間隔為0.1 s。

圖7 矩形水池中的水波

本文的模型是基于物理的，因此只需要為其提供一組初始條件，也就是模擬開始時(shí)水面的形狀，水面就會(huì)根據(jù)淺水方程自然演變。圖8是水從柱子間流出，較為真實(shí)地模擬了水體效果。

圖8 柱子間的水流

由于本文采用隱式半拉格朗日解法，可以使用較大時(shí)間步長進(jìn)行模擬，與較小時(shí)間步長相比，效果近似，但速度提升較快。圖9是時(shí)間步長分別為0.1 s和1.0 s的2個(gè)實(shí)驗(yàn)結(jié)果對(duì)比。在時(shí)間步長較大的情況下，隱式求解迭代次數(shù)增多，在保證體積守恒的情況下，依然可以保持實(shí)時(shí)模擬。時(shí)間步長為0.1 s和1.0 s時(shí)，模擬100幀的時(shí)間，分別為26 s和17 s。

圖9 不同時(shí)間步長 (左側(cè)為0.1 s，右側(cè)為1.0 s)

關(guān)于固體與水體之間的交互，從效果上看還是很好的，較好地模擬了固體對(duì)水體的作用力，如圖10所示。

圖10 水面上行駛的小船

4.2 比較與分析

顯式方法由于其簡單性，被用于流體模型的許多實(shí)現(xiàn)中。但其主要缺點(diǎn)是對(duì)時(shí)間步長有嚴(yán)格的限制，數(shù)值解可能與真實(shí)解呈指數(shù)發(fā)散。通過選擇隱式積分方法，保證了穩(wěn)定性不受梯度項(xiàng)大小的限制。通過采用半拉格朗日方法，還確保了時(shí)間步長不受CFL條件的限制。只要有足夠的精度估計(jì)出發(fā)點(diǎn)，可確保算法是穩(wěn)定的。比較該算法與采用半拉格朗日方法的算法的穩(wěn)定性，以2010年CHENTANEZ和MüLLER[12]的算法作為比較，其通過有限體積法同樣實(shí)現(xiàn)了體積守恒。兩者之間的效果對(duì)比如圖11所示，初始場景是3個(gè)水柱，2種方法的效果基本相同。在本文的測試案例中，隱式半拉格朗日方法得到的最大時(shí)間步長比文獻(xiàn)[12]方法大10倍以上，本文的CFL數(shù)在超過10時(shí)依然保持了模擬的穩(wěn)定性。

圖11 本文方法與文獻(xiàn)[12]方法的效果比較

進(jìn)一步研究了本文算法的體積守恒特性。在不增加或去除水的情況下，通過對(duì)水深進(jìn)行求和得出的總體積，在整個(gè)模擬過程中應(yīng)該保持不變。在本文所有的實(shí)驗(yàn)中，水的體積變化小于0.01%，由于變化極小，可認(rèn)為實(shí)現(xiàn)了體積守恒。與文獻(xiàn)[9]的方法相比，如圖12所示，使用2種方法模擬圖7的場景，其效果基本相同。不過同樣是半拉格朗日隱式求解，如圖13所示，本文方法保證了體積守恒。

圖12 本文方法與文獻(xiàn)[9]方法的效果對(duì)比

圖13 2種方法體積變化對(duì)比

4.3 進(jìn)一步討論

到目前為止，通過本文方法已經(jīng)實(shí)現(xiàn)了較大時(shí)間步長下的實(shí)時(shí)水體模擬，模擬過程中體積守恒，邊界處理令人滿意，單向流固耦合效果尚佳。但是水體模擬還有一些更有意思的事情，比如想要在水體的任意位置施加一個(gè)力、或者改變?nèi)我馑w的速度。對(duì)于上述用戶指定的約束條件，本文方法可以實(shí)現(xiàn)。根據(jù)淺水方程來看，如果要給水體加上一定的約束條件，那么在每一步模擬的最后，改變水體的速度是最好的選擇。外力等其他條件，可以先將其轉(zhuǎn)化為動(dòng)量通量，再轉(zhuǎn)化為速度的變化。通過該方法，可以更加靈活的模擬水體。圖14為本文方法模擬的風(fēng)對(duì)水面的作用力[13]。

圖14 風(fēng)對(duì)水面作用

5 結(jié) 論

本文提出了一個(gè)基于物理的水體模擬模型，該模型能夠模擬真實(shí)的水面，并且可以進(jìn)行有效的單向流固耦合，即使在較大的時(shí)間步長下也能保持穩(wěn)定，同時(shí)水體在模擬過程始終保持體積守恒。此外，因?yàn)榛跍\水方程，所以效率很高，可以做到實(shí)時(shí)模擬。不過，模型的局限性也在于此，只能模擬淺水的無黏流動(dòng)，對(duì)于更復(fù)雜的現(xiàn)象無法仿真。盡管如此，本模型在模擬溫和的海浪，對(duì)體積變化敏感的流體以及流固耦合方面還是可以起到很大的作用。下一步工作是實(shí)現(xiàn)雙向流固耦合，以及借助粒子的方法使得表面更富有細(xì)節(jié)且可以模擬一些破碎波浪等現(xiàn)象。

[1] JESCHKE S, WOJTAN C. Water wave animation via wavefrontparameter interpolation [J]. ACM Transactions on Graphics, 2015, 34(3): 1-14.

[2] 房克照, 尹晶, 孫家文, 等. 基于二維淺水方程的滑坡體興波數(shù)值模型[J]. 水科學(xué)進(jìn)展, 2017, 28(1): 96-105.

[3] MICHEL-DANSACV, BERTHON C, CLAINS, et al. A well-balanced scheme for the shallow-water equations with topography or Manning friction [J]. Journal of Computational Physics, 2017, 335: 115-154.

[4] 張迪, 繆小平, 彭福勝, 等. 一種求解對(duì)流擴(kuò)散方程的無條件穩(wěn)定算法[J]. 水動(dòng)力學(xué)研究與進(jìn)展: A輯, 2017, 32(2): 158-164.

[5] ROBERT A. A semi-Lagrangian and semi-implicit numerical integration scheme for the primitive meteorological equations [J]. Journalof the Meteorological Society of Japan Ser II, 1982, 60(1): 319-325.

[6] LENTINEM, GRéTARSSON J T, FEDKIW R. An unconditionally stable fully conservative semi-Lagrangian method [J]. Journal of Computational Physics, 2011, 230(8): 2857-2879.

[7] LENTINE M, AANJANEYA M, FEDKIW R. Mass and momentum conservation for fluid simulation [C]// Proceedings of the 2011 ACM SIGGRAPH/ Eurographics Symposium on Computer Animation - SCA '11. New York: ACM Press, 2011: 91-100.

[8] LAYTON A T, LAYTON H E. A semi-Lagrangiansemi-implicit numerical method for models of the urine concentrating mechanism [J]. SIAM Journal on Scientific Computing, 2002, 23(5): 1526-1548.

[9] LAYTON A T, VAN DE PANNE M. A numerically efficient and stable algorithm for animating water waves [J]. The Visual Computer, 2002, 18(1): 41-53.

[10] HARLOW F H, WELCH J E. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface [J]. Physics of Fluids, 1965, 8(12): 2182.

[11] QIAN S G, LI G, SHAO F J, et al. Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water flows in open channels [J]. Advances in Water Resources, 2018, 115: 172-184.

[12] CHENTANEZ N, MüLLER M. Real-time simulation of large bodies of water with small scale details [C]// Proceedings of the 2010 ACM SIGGRAPH/ Eurographics Symposium on Computer Animation.Goslar: Eurographics Association, 2010: 197-206.

[13] 鄒仲水, 趙棟梁, 黃健, 等. 海-氣界面動(dòng)量通量的估計(jì)方法分析與應(yīng)用[J]. 海洋學(xué)報(bào):中文版, 2014, 36(9): 75-83.

Real-Time Simulation of Water with Volume Conservation in Large Time Step

ZHU Kui-xin1, GAI Meng2,3, LAI Shun-nan2

(1. Peking University Shenzhen Graduate School, Shenzhen Guangdong 518055, China; 2. School of Electronics Engineering and Computer Science, Peking University, Beijing 100871, China; 3. Beijing Engineering Technology Research Center of Virtual Simulation and Visualization, Peking University, Beijing 100871, China)

Physics-based fluid simulation can capture the movement of water. A common real-time simulation method is based on two-dimensional shallow water equation. Firstly, we introduce the two-dimensional shallow water equation, and then propose a new solution method, which is calculated by an implicit semi-Lagrangian method with volume conservation. This method allows a large time step while maintaining stability. Finally, we add additional correction steps to ensure volume conservation in the calculation process. In addition, a fluid solid coupling method based on momentum conservation is proposed, which works well. The model can simulate water more realistically in real time and ensure that the water volume is always conserved.

shallow water equation; semi-Lagrangian; volume conservation; large time step; fluid solid coupling

TP 391

10.11996/JG.j.2095-302X.2019040725

A

2095-302X(2019)04-0725-08

2019-03-26；

定稿日期：2019-04-26

國家重點(diǎn)研發(fā)計(jì)劃課題(2017YFB1002705)；國家自然基金面上項(xiàng)目(61872398，61632003)；裝備預(yù)研基金項(xiàng)目(315050501)

朱奎鑫(1997-)，男，河南商丘人，碩士研究生。主要研究方向?yàn)榱黧w仿真。E-mail：zkx@pku.edu.cn

賴舜男(1965-)，女，重慶人，工程師，碩士。主要研究方向?yàn)樘摂M仿真。E-mail：snlai@pku.edu.cn